sfe_py的应力云图计算与显示step by step

子正

已于 2025-07-15 15:57:43 修改

阅读量88

点赞数

CC 4.0 BY-SA版权

分类专栏： # 预防性维护工程实践红宝书 - 结构篇问题建模 # 技术追踪文章标签：数学建模应力分析笔记 PHM

于 2025-07-14 17:40:53 首次发布

本文链接：https://blog.csdn.net/twicave/article/details/149335291

问题建模同时被 3 个专栏收录

166 篇文章

订阅专栏

技术追踪

77 篇文章

订阅专栏

预防性维护工程实践红宝书 - 结构篇

5 篇文章

订阅专栏

1.依赖环境

首先是相关环境的安装，如果你使用wsl模式的linux环境，你还是趁此机会给ubuntu装个虚拟显示器。相关的一系列安装动作：

pip install sfepy
1. 如果你安装的是2025.2版，你需要：
2. pip install jax
  1. 这是sfe_py矩阵运算的数学库依赖
  2. 大概率是numpy的速度不够快
sudo apt update && sudo apt install libgl1-mesa-glx x11-apps -y
1. 如果你使用了ubuntu 22.4以上的版本：
  1. sudo apt update && sudo apt install libglx-mesa0 x11-apps -y

echo "export DISPLAY=$(grep nameserver /etc/resolv.conf | awk '{print $2}'):0.0" >> ~/.bashrc
source ~/.bashrc
注意这里只是一次性地得到了宿主机的ip，如果需要长期使用x11，需要修改~/.bashrc:

# enable programmable completion features (you don't need to enable
# this, if it's already enabled in /etc/bash.bashrc and /etc/profile
# sources /etc/bash.bashrc).
#if [ -f /etc/bash_completion ] && ! shopt -oq posix; then
#    . /etc/bash_completion
#fi
#export DISPLAY=172.26.208.1:0.0
# 方法1：从路由表提取默认网关（宿主机IP）
host_ip_via_route=$(ip route | awk '/default/ {print $3}' | head -n1)

# 方法2：从 resolv.conf 提取DNS服务器（宿主机IP）
host_ip_via_dns=$(grep -oP 'nameserver \K[\d.]+' /etc/resolv.conf | head -n1)

# 优先使用路由表结果，失败则回退DNS解析
host_ip="${host_ip_via_route:-$host_ip_via_dns}"
export DISPLAY="${host_ip}:0.0"

windows环境安装：
1. vcxsrv-64.1.20.14.0.installer.exe

2. 例程

代码不建议上来徒手做.mesh文件建模，可以先跑一下系统自带的例程：

比如：

sfepy-run ./mixed_mesh.py

...................... .mixed_mesh.py 在sfepy的安装目录，如果你的python启用了.venv，它位置很好找：

find .venv -name "mixed_mesh.py"

sfepy-view beam_h* -f u:s0:wu:e:p0 u:s1:wu:e:p0

............................显示生成的.vtk文件

例程里的beam：

这个还没有仔细校对参数，应该是匀质梁的单端固定的效果。
接下来可以先从改变外力，叠加外部载荷开始。

2.1 还有一个可以徒手跟的例程

sfepy的手册，Page12~18页还列出了一个教程。是一个二维面受力引发的形变，建议手敲一遍，仔细理解每一行代码的具体含义。参见：

什么叫有限元分析？有限元分析软件有哪些？ - 知乎https://www.zhihu.com/question/49990690/answer/1928457972356461459

2.2 注释过的最简有限元解算源码

#you can run this .py by python3 ./regions_test.py

import numpy as nm
from sfepy.discrete.fem import Mesh, FEDomain, Field

#.mesh是原始模型，看到部分，由点，三角形，三角锥，6面体这类元素组成
mesh = Mesh.from_file('/home/git/xjtu_sy/.venv/lib/python3.12/site-packages/sfepy/meshes/2d/rectangle_tri.mesh')
#域可以认为是mesh文件的一个用于有限元结算的内部结构，它可能会有一些额外的数据单元，比如受力，比如形变等关联的数据域
domain = FEDomain('domain', mesh)
#[:,0]是这个点阵的x坐标，这里统计出了x坐标的最大最小值
min_x, max_x = domain.get_mesh_bounding_box()[:, 0]
#这是一个误差上限
eps = 1e-8 * (max_x - min_x)
#region似乎是一些点阵的子集，分区，。这里吧所有的点阵都加入了这个叫"Omega"的域
omega = domain.create_region('Omega', 'all')
#这里又拆分了两个分区，这似乎定义了x轴方向的一个非常薄的表面层，这是在模拟热处理过程吗？
gamma1 = domain.create_region('Gamma1','vertices in x < %.10f' % (min_x + eps),'facet')
gamma2 = domain.create_region('Gamma2','vertices in x > %.10f' % (max_x - eps),'facet')
#这个似乎是应力场，可能对应着下面的uvmf这四元
field = Field.from_args('fu', nm.float64, 'vector', omega, approx_order=2)

from sfepy.discrete import (FieldVariable, Material, Integral, Function,Equation, Equations, Problem)
u = FieldVariable('u', 'unknown', field)
v = FieldVariable('v', 'test', field, primary_var_name='u')

from sfepy.mechanics.matcoefs import stiffness_from_lame
m = Material('m', D=stiffness_from_lame(dim=2, lam=1.0, mu=1.0))
f = Material('f', val=[[0.02], [0.01]])

#三阶似乎是匀质切应力与最终挠度曲线的积分方程阶数
integral = Integral('i', order=3)
#下面似乎是内部约束
from sfepy.terms import Term
t1 = Term.new('dw_lin_elastic(m.D, v, u)',integral, omega, m=m, v=v, u=u)
t2 = Term.new('dw_volume_lvf(f.val, v)',integral, omega, f=f, v=v)
eq = Equation('balance', t1 + t2)
eqs = Equations([eq])

#这里是边界条件
from sfepy.discrete.conditions import Conditions, EssentialBC
fix_u = EssentialBC('fix_u', gamma1, {'u.all' : 0.0})
def shift_u_fun(ts, coors, bc=None, problem=None, shift=0.0):
      val = shift * coors[:,1]**2
      return val
bc_fun = Function('shift_u_fun', shift_u_fun,extra_args={'shift' : 0.01})
shift_u = EssentialBC('shift_u', gamma2, {'u.0' : bc_fun})

#最后是方程的求解时的一些配置，求解阶数，存至regions_test.vtk文件
#这里可以看到region的一些信息。
from sfepy.base.base import IndexedStruct
from sfepy.solvers.ls import ScipyDirect
from sfepy.solvers.nls import Newton
ls = ScipyDirect({})
nls_status = IndexedStruct()
nls = Newton({}, lin_solver=ls, status=nls_status)
pb = Problem('elasticity', equations=eqs)
pb.save_regions_as_groups('regions_test')
#! sfepy-view regions.vtk -2 --grid-vector1 "[2, 0, 0]"
#%history -f regions_test.py #save ipython input scripts.

#folded by force.
#代入边界条件后，后面才是求解的过程。这个原始的点阵会被拉伸
pb.set_bcs(ebcs=Conditions([fix_u, shift_u]))
pb.set_solver(nls)
status = IndexedStruct()
variables = pb.solve(status=status)
print('Nonlinear solver status:\n', nls_status)
print('Stationary solver status:\n', status)
pb.save_state('linear_elasticity.vtk', variables)
#sfepy-view linear_elasticity.vtk -2
#sfepy-view linear_elasticity.vtk -2 -f u:wu:p0 1:vw:p0

3.受力变化

例程里的受力是一种匀质梁一样的受力：

def get_force(ts, coors, mode=None, **kwargs):
    if mode == 'qp':
        F = 1e3

        val = nm.zeros_like(coors)[..., None]
        val[:, 2, 0] = -coors[:, 0] / 0.7 * F

        return {'val': val}

这里的F是受力
coors是所有点的坐标
然后coors[:,0]取的是所有的点的坐标。
val[:,2, 0]，第二个2，实际是z坐标，那个0大概是指压力，他可能会受到拉压扭力。

3.1 修改后的

def get_force(ts, coors, mode=None, **kwargs):
    if mode == 'qp':
        F = 8e4
        tol = 1e-5  # 容差范围，避免浮点误差
        
        # 计算中点坐标
        val = nm.zeros_like(coors)[..., None]
        x_mid = (nm.max(coors[:, 0]) + nm.min(coors[:, 0])) / 2.0
        #y_mid = (nm.max(coors[:, 1]) + nm.min(coors[:, 1])) / 2.0  # 若为三维模型
        #z_mid = (nm.max(coors[:, 2]) + nm.min(coors[:, 2])) / 2.0

        #仅有三点受力
        val[:, 2, 0] = 0

        # 中点受力
        dist = nm.abs(coors[:, 0] - x_mid)
        idx_center = nm.argmin(dist)  # 找到最近节点的索引
        #print(f'{x_mid} - {idx} - {nm.max(coors[:, 0])}, {nm.min(coors[:, 0])}')
        #os._exit(0)
        val[idx_center, 2, 0] = -F
        
        #左端点
        x_min = nm.min(coors[:,0])
        dist = nm.abs(coors[:, 0] - x_min)
        idx_left = nm.argmin(dist)  # 找到最近节点的索引
        val[idx_left, 2, 0] = F*0.5
        
        #右端点
        x_max = nm.max(coors[:,0])
        dist = nm.abs(coors[:, 0] - x_max)
        idx_right = nm.argmin(dist)  # 找到最近节点的索引
        val[idx_right, 2, 0] = F*0.5
        print(f'>>>>>>>>>>>>>>>>>>>>>>>>pos= {idx_left} {idx_center} {idx_right}')        
        #os._exit(0)
        return {'val': val}

我相当于做了一个简支梁，仅考虑中点处的外力，但是应力云图不对：
sfepy-view beam_h*.vtk -f u:s1:wu:e:p0 u:s0:wu:e:p

附录A 原始文档，例程位置

SfePy: Simple Finite Elements in Python — SfePy version: 2025.2+git.e2f7d954 documentation

有一篇关于mesh的文可能会用到：

mesh算法进化如何实现？这里是成功实例......_几何_三角形_数量

它给出了更少点数，更好的效果的一种实施例。

附录B 关于测点位置与剪力滞

实际埋点时，剪力滞的问题可能会影响误差，目前已知的解决方案是，把传感器放置在腹板，而不是上下翼板，避开这个问题。

附录C .mesh文件的格式

Dimension 3 定义了这是一个3维结构
Vertices是点，每个点因为是3维，所以有3个点，第4列是点类目标的标记。
Hexahedra是六面体，由8个点组成，最后一列是类型标记，1.
Tetrahedra是三棱锥，也就是4面体，由4个点组成，最后一列是类型标记，2.
六面体，三棱锥中的坐标都是顶点坐标，Vertices.单独的点没有意义。只能在纳入一个结构中才有意义。

MeshVersionFormatted 2
Dimension 3

Vertices
32
0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0
0.0000000000000000e+00 1.0000000000000001e-01 0.0000000000000000e+00 0
0.0000000000000000e+00 1.0000000000000001e-01 1.0000000000000001e-01 0
0.0000000000000000e+00 0.0000000000000000e+00 1.0000000000000001e-01 0
1.0000000000000001e-01 0.0000000000000000e+00 0.0000000000000000e+00 0
1.0000000000000001e-01 1.0000000000000001e-01 0.0000000000000000e+00 0
1.0000000000000001e-01 1.0000000000000001e-01 1.0000000000000001e-01 0
1.0000000000000001e-01 0.0000000000000000e+00 1.0000000000000001e-01 0
2.0000000000000001e-01 0.0000000000000000e+00 0.0000000000000000e+00 0
2.0000000000000001e-01 1.0000000000000001e-01 0.0000000000000000e+00 0
2.0000000000000001e-01 1.0000000000000001e-01 1.0000000000000001e-01 0
2.0000000000000001e-01 0.0000000000000000e+00 1.0000000000000001e-01 0
2.9999999999999999e-01 0.0000000000000000e+00 0.0000000000000000e+00 0
2.9999999999999999e-01 1.0000000000000001e-01 0.0000000000000000e+00 0
2.9999999999999999e-01 1.0000000000000001e-01 1.0000000000000001e-01 0
2.9999999999999999e-01 0.0000000000000000e+00 1.0000000000000001e-01 0
4.0000000000000002e-01 0.0000000000000000e+00 0.0000000000000000e+00 0
4.0000000000000002e-01 1.0000000000000001e-01 0.0000000000000000e+00 0
4.0000000000000002e-01 1.0000000000000001e-01 1.0000000000000001e-01 0
4.0000000000000002e-01 0.0000000000000000e+00 1.0000000000000001e-01 0
5.0000000000000000e-01 0.0000000000000000e+00 0.0000000000000000e+00 0
5.0000000000000000e-01 1.0000000000000001e-01 0.0000000000000000e+00 0
5.0000000000000000e-01 1.0000000000000001e-01 1.0000000000000001e-01 0
5.0000000000000000e-01 0.0000000000000000e+00 1.0000000000000001e-01 0
5.9999999999999998e-01 0.0000000000000000e+00 0.0000000000000000e+00 0
5.9999999999999998e-01 1.0000000000000001e-01 0.0000000000000000e+00 0
5.9999999999999998e-01 1.0000000000000001e-01 1.0000000000000001e-01 0
5.9999999999999998e-01 0.0000000000000000e+00 1.0000000000000001e-01 0
6.9999999999999996e-01 0.0000000000000000e+00 0.0000000000000000e+00 0
6.9999999999999996e-01 1.0000000000000001e-01 0.0000000000000000e+00 0
6.9999999999999996e-01 1.0000000000000001e-01 1.0000000000000001e-01 0
6.9999999999999996e-01 0.0000000000000000e+00 1.0000000000000001e-01 0

Hexahedra
5
1 2 3 4 5 6 7 8 1
5 6 7 8 9 10 11 12 1
9 10 11 12 13 14 15 16 1
13 14 15 16 17 18 19 20 1
21 22 23 24 25 26 27 28 1

Tetrahedra
12
17 18 21 20 2
18 24 21 20 2
18 22 21 24 2
18 19 22 20 2
19 24 22 20 2
19 23 22 24 2
25 26 29 28 2
26 32 29 28 2
26 30 29 32 2
26 27 30 28 2
27 32 30 28 2
27 31 30 32 2

End