🎯要点

1. 数学​方法​：🎯一维线性边界值问题：🖊高斯求积法则 | 🖊洛巴托求积法则 | 🖊矩阵插值和微分计算 | 🖊在细化网格上生成值。🎯二维边界值问题：构建二维网格：🖊生成几何边界 | 🖊生成扩展转置连接矩阵 | 🖊生成互补网格 | 🖊绘制标记网格图和互补标记网格 | 🖊从定义几何尺寸生成网格求解：使用逐元素策略、部分（或全部）矢量化方法或稀疏矩阵的坐标形式解刚度矩阵，使用全部矢量化方法解全局力矢量，计算方程系数，计算狄利克雷边界条件，解有限元矩阵和偏微分方程。
2. 构件分析方法​：🎯矩阵分析：🖊平面桁架位移和形变状态：提取结构的刚度矩阵，绘制变形形状，节点位移后，计算杆件的内力和桁架 | 🖊 空间桁架位移和形变状态：提取结构的刚度矩阵，计算节点载荷下 | 🖊平面框架位移和形变状态：提取局部刚度矩阵和旋转矩阵，计算出内力，计算钢框架承载下 | 🖊空间框架位移和形变状态：提取结构刚度矩阵，绘制变形形状，计算钢框架承载下，计算承受横向载荷和力矩的空间钢框架 | 🖊网格结构位移和形变：提取结构刚度矩阵，绘制形变形状，计算承受集中力矩和点载荷，计算钢具有相同的圆形横截面网架结构承受两个集中力矩 | 🖊构件载荷位移和形变：承受集中和均匀荷载的钢网架结构 | 🖊支撑沉降平面框架位移：外力影响下形变 | 🖊温度因素产生的构件位移：构件的拉伸或压缩的应力 | 🖊非标件框架位移：提取结构的刚度矩阵，计算给定载荷下 | 🖊释放构件位移：节点载荷情况下 | 🖊弹性支撑构件位移：均匀分布载荷的节点。🎯线性弹性分析：🖊双节点和三节点桁架：计算刚度矩阵 | 🖊弹簧和不同横截面的一维结构：计算外力条件下的位移，计算节点载荷下位移和法向应力 | 🖊平面桁架：计算外力下的形变 | 🖊网格二维结构：计算平面应力、平面应变和轴对称单元的刚度矩阵，计算节点力矢量，绘制形变 | 🖊平面应力计算情景：受拉伸表面载荷作用下，均匀弯曲和表面牵引载荷作用下，平面应变结构，轴对称结构。🎯弹塑性分析 | 🎯有限变形和超弹性 | 🎯有限应变。

🍇Python平面应力下变形位移求解器

考虑长度为 50 毫米、宽度为 10 毫米、厚度为 1 毫米的条带。使用 200 N 的力，该条带在垂直的平面应变张力下变形。条带的底部保持固定。带材的弹性特性为E=100MPa，ν=0.48。

在此，我使用线性弹性方程：
$${\sigma }_{ij}=2\mu {\epsilon }_{ij}+\lambda {\epsilon }_{kk}{\delta }_{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_nxconn.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).TdN = 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.0K[2*i+1,:]   = 0.0K[2*i,2*i]   = 1.0K[2*i+1,2*i+1] = 1.0if nodes[i,1] == mesh_ly:x = nodes[i,0]f[2*i+1] = 20.0if x == 0.0 or x == mesh_lx:f[2*i+1] *= 0.5print('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，与上面的近似计算结果一致。

参阅一：计算思维

参阅二：亚图跨际