弹性力学基本方程及有限元编程
-
概述
在工程领域，无论是设计一座桥梁、制造一台精密机械，还是研发航空航天部件，都需要弄清楚材料在受力时会发生怎样的变化。弹性力学就是专门研究这个问题的学科，它通过微分平衡方程、几何方程和物理方程，为我们揭示材料内部应力、应变和位移的规律。而有限元方法，则是帮助我们把这些理论应用到实际工程计算中的强大工具。
提醒：个别地方公式识别有问题，感兴趣的读者可以到本人公众号阅读 冬生亦冬生（全网首发，欢迎一起交流学习）
一、弹性力学三大基本方程
微分平衡方程（力的平衡）
描述弹性体微元的力学平衡，基于小变形、连续均匀各向同性假设：
应用：求解结构应力场（如机械零件、建筑承重结构强度分析）。
几何方程（变形几何关系）
建立位移与应变的线性关系（小变形假设）：
线应变： 切应变：（其余切应变类似）
应用举例：在分析建筑物在地震中的晃动时，通过测量建筑物各部分的位移，利用几何方程就能计算出建筑物的应变，判断是否会出现裂缝等损坏。
物理方程（本构关系）
各向同性弹性材料的胡克定律（矩阵形式）：
应用：连接应力与应变，用于材料选型和实验应力计算。
二、边界条件
位移边界条件 指定边界位移值：在边界（如固定端。
应力边界条件 指定边界面力：在边界 上，
（(l,m,n) 为边界外法线方向余弦，为面力分量）。 3. 混合边界条件 部分边界给定位移，部分边界给定应力（如简支梁两端约束）。
三、方程联合求解
三大方程与边界条件联立，形成定解问题： 几何方程将位移 转换为应变； 物理方程将应变转换为应力； 微分平衡方程结合体积力建立平衡关系； 边界条件约束解的唯一性。典型问题：简支梁受均布载荷，通过假设位移函数、代入方程并满足边界条件求解应力与挠度。
四、简单弹性问题的有限元编程方法
有限元方法是一种将复杂问题简单化的数值计算方法，它把连续的弹性体划分成许多小的单元，通过计算每个单元的特性，再把它们组合起来，得到整个弹性体的结果。下面以平面应力三角形单元为例，介绍有限元编程的核心步骤：
问题离散化 将弹性体划分为个单元（如三角形、四边形单元），节点总数为，节点位移向量（以平面问题为例）。
单元刚度矩阵 以平面应力三角形单元（3 节点，坐标 为例： 形函数为单元面积，为几何系数）；
应变 - 位移关系 为应变矩阵）；
单元刚度矩阵
整体刚度矩阵组装 将所有单元刚度矩阵按节点编号叠加，形成整体刚度矩阵 ，尺寸为。
载荷与边界条件处理 节点载荷向量 包含集中力、分布力等效节点力； 位移边界条件：直接约束中对应自由度（如固定节点设为 0）； 应力边界条件：转化为节点载荷加入。
求解系统方程 通过稀疏矩阵求解器（如高斯消元法）解得节点位移，再计算应变和应力。 简单 Python 代码框架（平面应力三角形单元）
import numpy as np
import scipy.sparse as sp
import scipy.sparse.linalg as spla
import matplotlib.pyplot as plt
import matplotlib.tri as tri
# ============================
# 1. 网格生成模块
# ============================
def generate_mesh(L=10, H=2, nx=10, ny=4):
"""
生成矩形区域的三角形网格
参数：
L : 长度 (mm)
H : 高度 (mm)
nx : x方向单元数
ny : y方向单元数
返回：
nodes : 节点坐标数组 (N,2)
elements : 单元节点索引数组 (M,3)
"""
# 生成节点
x = np.linspace(0, L, nx+1)
y = np.linspace(0, H, ny+1)
xx, yy = np.meshgrid(x, y)
nodes = np.column_stack([xx.ravel(), yy.ravel()])
# 生成三角形单元
elements = []
for i in range(ny):
for j in range(nx):
# 单个矩形分为两个三角形
n0 = i*(nx+1) + j
n1 = n0 + 1
n2 = n0 + (nx+1) + 1
n3 = n0 + (nx+1)
# 添加两个三角形单元
elements.append([n0, n1, n2])
elements.append([n0, n2, n3])
return np.array(nodes), np.array(elements)
# ============================
# 2. 材料参数和单元刚度矩阵计算
# ============================
def linear_triangle_stiffness(nodes, E=210e3, nu=0.3, t=1.0):
"""
计算线性三角形单元的刚度矩阵（平面应力）
参数：
nodes : 单元节点坐标 (3,2)
E : 弹性模量 (MPa)
nu : 泊松比
t : 厚度 (mm)
返回：
ke : 单元刚度矩阵 (6,6)
"""
# 节点坐标
x1, y1 = nodes[0]
x2, y2 = nodes[1]
x3, y3 = nodes[2]
# 计算单元面积
A = 0.5 * abs((x2-x1)*(y3-y1) - (x3-x1)*(y2-y1))
# 弹性矩阵D
D = E/(1-nu**2) * np.array([
[1, nu, 0 ],
[nu, 1, 0 ],
[0, 0, (1-nu)/2]
])
# 应变位移矩阵B
b = np.array([y2-y3, y3-y1, y1-y2])
c = np.array([x3-x2, x1-x3, x2-x1])
B = 1/(2*A) * np.array([
[b[0], 0, b[1], 0, b[2], 0 ],
[0, c[0], 0, c[1], 0, c[2]],
[c[0], b[0], c[1], b[1], c[2], b[2]]
])
# 单元刚度矩阵
ke = (B.T @ D @ B) * A * t
return ke
# ============================
# 3. 总刚度矩阵组装
# ============================
def assemble_global_stiffness(nodes, elements, E, nu, t):
"""
组装全局刚度矩阵
参数：
nodes : 节点坐标数组
elements : 单元连接数组
E, nu, t : 材料参数
返回：
K : 全局刚度矩阵 (2N,2N) 稀疏格式
"""
n_nodes = nodes.shape[0]
n_dof = 2 * n_nodes
K = sp.lil_matrix((n_dof, n_dof))
for elem in elements:
# 获取单元节点索引
n1, n2, n3 = elem
# 计算单元刚度矩阵
elem_nodes = nodes[[n1, n2, n3]]
ke = linear_triangle_stiffness(elem_nodes, E, nu, t)
# 组装到全局矩阵
dofs = [2*n1, 2*n1+1, 2*n2, 2*n2+1, 2*n3, 2*n3+1]
for i in range(6):
for j in range(6):
K[dofs[i], dofs[j]] += ke[i,j]
return K.tocsr()
# ============================
# 4. 边界条件和载荷施加
# ============================
def apply_boundary_conditions(K, F, nodes, fix_left=True):
"""
应用边界条件（左侧固定，右侧施加集中力）
返回：
K_modified : 处理后的刚度矩阵
F_modified : 处理后的载荷向量
"""
# 找到左侧节点（x=0）
left_nodes = np.where(nodes[:,0] == 0)[0]
fixed_dofs = []
for n in left_nodes:
fixed_dofs.extend([2*n, 2*n+1]) # 固定x和y方向
# 找到右侧节点（x=max）
right_nodes = np.where(nodes[:,0] == nodes[:,0].max())[0]
for n in right_nodes:
F[2*n+1] = -100.0 # 施加y方向集中力
# 处理固定自由度
dofs = np.arange(K.shape[0])
free_dofs = np.setdiff1d(dofs, fixed_dofs)
K_modified = K[free_dofs, :][:, free_dofs]
F_modified = F[free_dofs]
return K_modified, F_modified, free_dofs, fixed_dofs
# ============================
# 5. 求解和后处理
# ============================
def post_process(nodes, elements, U, scale=100):
"""
后处理：显示变形前后的网格和应力分布
参数：
scale : 变形放大系数
"""
# 创建三角剖分对象
triang = tri.Triangulation(nodes[:,0], nodes[:,1], elements)
# 绘制原始网格
plt.figure(figsize=(12,4))
plt.subplot(121)
plt.triplot(triang, color='gray', lw=0.5)
plt.title("Original Mesh")
# 绘制变形后网格
plt.subplot(122)
deformed_nodes = nodes + U[::2].reshape(-1,1) * scale # 仅显示x方向位移
deformed_triang = tri.Triangulation(deformed_nodes[:,0],
deformed_nodes[:,1],
elements)
plt.triplot(deformed_triang, color='red', lw=0.5)
plt.title(f"Deformed Mesh (Scale: {scale}X)")
plt.show()
# ============================
# 主程序
# ============================
if __name__ == "__main__":
# 参数设置
E = 210000 # 弹性模量 (MPa)
nu = 0.3 # 泊松比
t = 10.0 # 厚度 (mm)
L = 100 # 长度 (mm)
H = 20 # 高度 (mm)
nx, ny = 20, 4 # 网格划分
# 1. 生成网格
nodes, elements = generate_mesh(L, H, nx, ny)
print(f"Generated mesh with {nodes.shape[0]} nodes and {elements.shape[0]} elements")
# 2. 组装总刚度矩阵
K = assemble_global_stiffness(nodes, elements, E, nu, t)
print("Global stiffness matrix assembled")
# 3. 初始化载荷向量
F = np.zeros(2 * nodes.shape[0])
# 4. 应用边界条件
K_mod, F_mod, free_dofs, fixed_dofs = apply_boundary_conditions(K, F, nodes)
print(f"Applied boundary conditions: fixed {len(fixed_dofs)} DOFs")
# 5. 求解方程组
U = np.zeros(2 * nodes.shape[0])
U[free_dofs] = spla.spsolve(K_mod, F_mod)
print("Displacement solution completed")
# 6. 后处理
post_process(nodes, elements, U, scale=500)
# 显示最大位移
max_ux = np.max(np.abs(U[::2]))
max_uy = np.max(np.abs(U[1::2]))
print(f"Maximum displacements: Ux={max_ux:.4f} mm, Uy={max_uy:.4f} mm")
五、总结
弹性力学三大方程与边界条件构成理论基础，有限元方法通过离散化实现工程问题的数值求解。实际应用中需根据结构复杂度选择单元类型（如梁、壳、三维实体单元），结合商业软件（ABAQUS、ANSYS）或自主编程实现精确分析，广泛应用于航空航天、机械、土木等领域的结构设计与优化。由于变分法是弹性力学的重要内容也是有限元的重要基础，有利于理解有限元的构造及编程方法，后续将介绍这部分内容。提醒：个别地方公式识别有问题，感兴趣的读者可以到本人公众号阅读 冬生亦冬生（全网首发，欢迎一起交流学习）