python 悬臂梁的有限元分析

依赖包

fenics是一种用于有限元计算的动态面向对象库，它提供了一种专用的数学语言UFL来表述变分形式，并自动生成底层C++代码。

fenics 名称释义：

• fe：finite element的简写
• cs：computational software的简写
• ni：有了fe和cs后，由于最初fenics软件是在芝加哥大学(简称为phoenix）编译的，故而在其间加入ni就很自然而然了

fenics 库实现

fenics 匹配的python 2.0，反正3.0以上的无法调包成功

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from fenics import *

# 定义梁的参数
L = 5.0  # 梁长度 (m)
b = 0.1  # 梁宽度 (m)
h = 0.2  # 梁高度 (m)
E = 2e11  # 弹性模量 (Pa)
F = 1000  # 施加的集中力 (N)
  
# 创建网格
nx, ny = 20, 5
mesh = RectangleMesh(Point(0, 0), Point(L, h), nx, ny)

# 定义函数空间
V = FiniteElement('Lagrange', mesh.ufl_cell(), 1)
u = TrialFunction(V)
v = TestFunction(V)

# 定义变分问题
a = dot(grad(u), grad(v)) * dx
L = Constant(F) * v * ds(2)

# 求解变分问题
u = Function(V)
solve(a == L, u)

# 计算位移和应力
x = mesh.coordinates()[:, 0]
y = mesh.coordinates()[:, 1]
u_values = u.vector().get_local()
sigma_xx = E * diff(u, 'x')

# 绘制梁的变形和应力分布
fig = plt.figure(figsize=(12, 6))


# 绘制变形图
ax1 = fig.add_subplot(121, projection='3d')
ax1.plot_trisurf(x, y, u_values, cmap='viridis')
ax1.set_xlabel('x (m)')
ax1.set_ylabel('y (m)')
ax1.set_zlabel('Displacement (m)')
ax1.set_title('Deformation of the Cantilever Beam')

# 绘制应力分布图
ax2 = fig.add_subplot(122)
ax2.plot(x, sigma_xx)
ax2.set_xlabel('x (m)')
ax2.set_ylabel('Stress (Pa)')
ax2.set_title('Stress Distribution in the Cantilever Beam')
plt.tight_layout()
plt.show()

有限元分析代码实现

# @Time:2024/10/17
# @Autho:Jiang Gang
# @Version:1.0
# @Contact:sivxiu@foxmail.com

import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import solve

# 定义悬臂梁的参数
L = 1.0  # 梁长度 (m)
b = 0.2  # 梁宽度 (m)
h = 0.2  # 梁高度 (m)
E = 250e6  # 杨氏模量 (Pa)，泊松比0.3，可以不需要密度
F = 10000  # 末端集中力 (N)
n = 100  # 节点数

# 计算截面属性
A = b * h  # 刚体横截面积 (m^2)
I = b * h ** 3 / 12  # 截面惯性矩 (m^4)


# 离散化梁
x = np.linspace(0, L, n) # 沿梁长度方向的节点坐标
y = np.linspace(-h/2, h/2, n)
X, Y = np.meshgrid(x, y) # 二维网格


# 建立刚度矩阵
K = np.zeros((n * n, n * n))
for i in range(n):
    for j in range(n):
        K[i * n + j, i * n + j] = 12 * E * I / (b * h ** 3)
        if i > 0:
            K[i * n + j, (i - 1) * n + j] = 6 * E * I / (b * h ** 2)
            K[(i - 1) * n + j, i * n + j] = 6 * E * I / (b * h ** 2)
K[:n, :n] = 1.0  # 固定端约束

# 计算位移 u，载荷向量 F_vec
F_vec = np.zeros(n * n)
F_vec[-n:] = -F
u = solve(K, F_vec).reshape((n, n)) # 求解线性方程组 K * u = F_vec

# 计算应力 simga
M = -F * (L - X)  # 弯矩
sigma = M * Y / I  # 应力 (Pa)
epsilon = sigma / E  # 应变

# 绘制应力云图
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 6))

im = ax1.imshow(np.abs (sigma) / 1e6, cmap='coolwarm', extent=(0, L, -h/2, h/2)) # imshow() 函数绘制二维应力分布,sigma为应力
ax1.set_xlabel('Position (m)')
ax1.set_ylabel('Height (m)')
ax1.set_title('Stress in Cantilever Beam')
cbar = fig.colorbar(im, ax=ax1) # 添加颜色条和标签
cbar.set_label('Stress (MPa)')

im2 = ax2.imshow(np.abs(epsilon), cmap='coolwarm', extent=(0, L, -h/2, h/2))
ax2.set_xlabel('Position (m)')
ax2.set_ylabel('Height (m)')
ax2.set_title('Strain in Cantilever Beam')
cbar2 = fig.colorbar(im2, ax=ax2)
cbar2.set_label('Strain (%)')

plt.tight_layout()
plt.show()

代码结果，没有显示变形

应力

应变

使用ansys模拟相同边界和条件，除去应力奇异，仿真结果也是相当的给力，等效应力和等效应变

等效应力

等效应变

应力集中

Von Mises应力

最大应力理论或者Rankine理论：最大或者最小应力达到材料的屈服或者极限强度时，材料失效

• 对于延展性强的材料，该理论失效
• 脆性材料的失效理论还有Coulomb-Mohr理论

对于延性材料 ductile material（韧性材料），静水应力hydrostatic stress不会导致材料失效

• 韧性材料屈服的机制：剪切变形，屈服是由引起形变的应力导致的
• 莫尔圆学习
• Tresca失效准则：最大剪切应力理论，当最大剪切应力等于拉伸试验中屈服时的剪应力时，发生屈服
  • Tresca失效准则的屈服面完全位于Von Mises的屈服面内，Tresca评价更加保守，二者差异约为15.5%
• Von Mises失效准则：最大畸变能理论，当材料的最大变形能等于屈服时的变形能时，材料屈服；Von Mises更常用
  • 单元的变形能
    
  • 当Von Mises应力大于材料的屈服强度时，预计发生屈服

等效应力

• 范式等效应力（Von Mises Stress）是一种屈服准则，习惯称为Mises等效应力，它遵循材料力学第四强度理论(形状改变比能理论)
• stress intensity（应力强度），是由第三强度理论得到的当量应力，其值为第一主应力减去第三主应力
  • 三个方向的应力排列，最小应力为0
• 一般脆性材料，如铸铁、石料、混凝土，多用第一强度理论。考察绝对值最大的主应力。
• 一般材料在外力作用下产生塑性变形，以流动形式破坏时，应该采用第三或第四强度理论。压力容器上用第三强度理论（安全第一），其它多用第四强度理论

屈服准则

• A.受力物体内质点处于单向应力状态时，只要单向应力大到材料的屈服点时，则该质点开始由弹性状态进入塑性状态，即处于屈服。
• B.受力物体内质点处于多向应力状态时，必须同时考虑所有的应力分量。在一定的变形条件（变形温度、变形速度等） 下，只有当各应力分量之间符合一定关系时，质点才开始进入塑性状态，这种关系称为屈服准则，也称塑性条件。它是描述受力物体中不同应力状态下的质点进入塑性状态并使塑性变形继续进行所必须遵守的力学条件

第三强度理论认为最大剪应力是引起流动破坏的主要原因,如低碳钢拉伸时在与轴线成45度的截面上发生最大剪应力,材料沿着这个平面发生滑移,出现滑移线。这一理论比较好的解释了塑性材料出现塑性变形的现象。

第四强度理论认为形状改变比能是引起材料流动破坏的主要原因,钢材等塑性材料遵循第四强度理论,结果更符合实际。

工程应力-应变曲线

金属或者弹性体的应力-应变曲线，通常分为四个阶段： 弹性阶段、屈服阶段、应变硬化阶段和颈缩断裂阶段

塑件材料在拉伸试验进入屈服阶段以后，开始产生显著的塑性变形，其数值远比弹性变形大。此外，试件横截面也渐渐变小。进入强化阶段后，试件伸长和横截面收缩就更加明显。特别是在局部变形阶段，试件颈缩部分的拉伸应变比其余各处大，截面面积也与其余各处明显不同。

真应变和真应力比工程应变和工程应力更准确，但是我们却很少需要真实应变和真实应力。由于负荷值的变化随时可以读出，但瞬间截面积很难直接读出，因此，工程应力应变曲线较容易得到。

塑性形变过程中，假设体积不变，真应变${ϵ}_t$与工程应变$ϵ$：$${ϵ}_t=\int d{ϵ}_t={\int }_{l_0}^l\frac{dl}{l}=\ln \frac{l}{l_0}=\ln (1+ϵ)$$
其中，$l_0$为初始标距值，$l$为标距内瞬时长度值，$A_0$为初始横截面积，$A$为瞬时横截面积

由此，真应力${\sigma }_t$与工程应力$\sigma$：$${\sigma }_t=\frac{F}{A}=\frac{Fl}{A_0{l}_0}=\frac{F}{A_0}(1+ϵ)=\sigma (1+ϵ)$$
真应变或真应力均比对应的工程值要大