🎯要点
- 数学方法：🎯一维线性边界值问题：🖊高斯求积法则 | 🖊洛巴托求积法则 | 🖊矩阵插值和微分计算 | 🖊在细化网格上生成值。🎯二维边界值问题：构建二维网格：🖊生成几何边界 | 🖊生成扩展转置连接矩阵 | 🖊生成互补网格 | 🖊绘制标记网格图和互补标记网格 | 🖊从定义几何尺寸生成网格求解：使用逐元素策略、部分（或全部）矢量化方法或稀疏矩阵的坐标形式解刚度矩阵，使用全部矢量化方法解全局力矢量，计算方程系数，计算狄利克雷边界条件，解有限元矩阵和偏微分方程。
- 构件分析方法：🎯矩阵分析：🖊平面桁架位移和形变状态：提取结构的刚度矩阵，绘制变形形状，节点位移后，计算杆件的内力和桁架 | 🖊 空间桁架位移和形变状态：提取结构的刚度矩阵，计算节点载荷下 | 🖊平面框架位移和形变状态：提取局部刚度矩阵和旋转矩阵，计算出内力，计算钢框架承载下 | 🖊空间框架位移和形变状态：提取结构刚度矩阵，绘制变形形状，计算钢框架承载下，计算承受横向载荷和力矩的空间钢框架 | 🖊网格结构位移和形变：提取结构刚度矩阵，绘制形变形状，计算承受集中力矩和点载荷，计算钢具有相同的圆形横截面网架结构承受两个集中力矩 | 🖊构件载荷位移和形变：承受集中和均匀荷载的钢网架结构 | 🖊支撑沉降平面框架位移：外力影响下形变 | 🖊温度因素产生的构件位移：构件的拉伸或压缩的应力 | 🖊非标件框架位移：提取结构的刚度矩阵，计算给定载荷下 | 🖊释放构件位移：节点载荷情况下 | 🖊弹性支撑构件位移：均匀分布载荷的节点。🎯线性弹性分析：🖊双节点和三节点桁架：计算刚度矩阵 | 🖊弹簧和不同横截面的一维结构：计算外力条件下的位移，计算节点载荷下位移和法向应力 | 🖊平面桁架：计算外力下的形变 | 🖊网格二维结构：计算平面应力、平面应变和轴对称单元的刚度矩阵，计算节点力矢量，绘制形变 | 🖊平面应力计算情景：受拉伸表面载荷作用下，均匀弯曲和表面牵引载荷作用下，平面应变结构，轴对称结构。🎯弹塑性分析 | 🎯有限变形和超弹性 | 🎯有限应变。
🍇Python平面应力下变形位移求解器
考虑长度为 50 毫米、宽度为 10 毫米、厚度为 1 毫米的条带。使用 200 N 的力，该条带在垂直的平面应变张力下变形。条带的底部保持固定。带材的弹性特性为E=100MPa，ν=0.48。
在此，我使用线性弹性方程：
σ i j = 2 μ ε i j + λ ε k k δ i j \sigma_{i j}=2 \mu \varepsilon_{i j}+\lambda \varepsilon_{k k} \delta_{i j} σij=2μεij+λεkkδij
拉梅常数可以根据已知的弹性性质来确定。通过分配平面应变并假设 1 方向的自由位移，我可以执行快速计算，得出：最大顶部位移 = 8 毫米。
import numpy as np
import math
from matplotlib import pyplot as plt
def shape(xi):
x,y = tuple(xi)
N = [(1.0-x)*(1.0-y), (1.0+x)*(1.0-y), (1.0+x)*(1.0+y), (1.0-x)*(1.0+y)]
return 0.25*np.array(N)
print('create mesh')
mesh_ex = 9
mesh_ey = 49
mesh_lx = 10.0
mesh_ly = 50.0
# derived
mesh_nx = mesh_ex + 1
mesh_ny = mesh_ey + 1
num_nodes = mesh_nx * mesh_ny
num_elements = mesh_ex * mesh_ey
mesh_hx = mesh_lx / mesh_ex
mesh_hy = mesh_ly / mesh_ey
nodes = []
for y in np.linspace(0.0, mesh_ly, mesh_ny):
for x in np.linspace(0.0, mesh_lx, mesh_nx):
nodes.append([x,y])
nodes = np.array(nodes)
conn = []
for j in range(mesh_ey):
for i in range(mesh_ex):
n0 = i + j*mesh_nx
conn.append([n0, n0 + 1, n0 + 1 + mesh_nx, n0 + mesh_nx])
print ('material model - plane strain')
E = 100.0
v = 0.48
C = E/(1.0+v)/(1.0-2.0*v) * np.array([[1.0-v, v, 0.0],
[ v, 1.0-v, 0.0],
[ 0.0, 0.0, 0.5-v]])
print('create global stiffness matrix')
K = np.zeros((2*num_nodes, 2*num_nodes))
q4 = [[x/math.sqrt(3.0),y/math.sqrt(3.0)] for y in [-1.0,1.0] for x in [-1.0,1.0]]
B = np.zeros((3,8))
for c in conn:
xIe = nodes[c,:]
Ke = np.zeros((8,8))
for q in q4:
dN = gradshape(q)
J = np.dot(dN, xIe).T
dN = np.dot(np.linalg.inv(J), dN)
B[0,0::2] = dN[0,:]
B[1,1::2] = dN[1,:]
B[2,0::2] = dN[1,:]
B[2,1::2] = dN[0,:]
Ke += np.dot(np.dot(B.T,C),B) * np.linalg.det(J)
for i,I in enumerate(c):
for j,J in enumerate(c):
K[2*I,2*J] += Ke[2*i,2*j]
K[2*I+1,2*J] += Ke[2*i+1,2*j]
K[2*I+1,2*J+1] += Ke[2*i+1,2*j+1]
K[2*I,2*J+1] += Ke[2*i,2*j+1]
print('assign nodal forces and boundary conditions')
f = np.zeros((2*num_nodes))
for i in range(num_nodes):
if nodes[i,1] == 0.0:
K[2*i,:] = 0.0
K[2*i+1,:] = 0.0
K[2*i,2*i] = 1.0
K[2*i+1,2*i+1] = 1.0
if nodes[i,1] == mesh_ly:
x = nodes[i,0]
f[2*i+1] = 20.0
if x == 0.0 or x == mesh_lx:
f[2*i+1] *= 0.5
print('solving linear system')
u = np.linalg.solve(K, f)
print('max u=', max(u))
print('plotting displacement')
ux = np.reshape(u[0::2], (mesh_ny,mesh_nx))
uy = np.reshape(u[1::2], (mesh_ny,mesh_nx))
xvec = []
yvec = []
res = []
for i in range(mesh_nx):
for j in range(mesh_ny):
xvec.append(i*mesh_hx + ux[j,i])
yvec.append(j*mesh_hy + uy[j,i])
res.append(uy[j,i])
t = plt.tricontourf(xvec, yvec, res, levels=14, cmap=plt.cm.jet)
plt.scatter(xvec, yvec, marker='o', c='b', s=2)
plt.grid()
plt.colorbar(t)
plt.axis('equal')
plt.show()
下图显示了运行 Python 代码的结果，最大位移为6.75 mm，与上面的近似计算结果一致。