概述
在工程领域，无论是设计一座桥梁、制造一台精密机械，还是研发航空航天部件，都需要弄清楚材料在受力时会发生怎样的变化。弹性力学就是专门研究这个问题的学科，它通过微分平衡方程、几何方程和物理方程，为我们揭示材料内部应力、应变和位移的规律。而有限元方法，则是帮助我们把这些理论应用到实际工程计算中的强大工具。
提醒：个别地方公式识别有问题，感兴趣的读者可以到本人公众号阅读 冬生亦冬生（全网首发，欢迎一起交流学习）
一、弹性力学三大基本方程
-
微分平衡方程（力的平衡）
描述弹性体微元的力学平衡，基于小变形、连续均匀各向同性假设：
应用：求解结构应力场（如机械零件、建筑承重结构强度分析）。
-
几何方程（变形几何关系）
建立位移与应变的线性关系（小变形假设）：
线应变： 切应变：（其余切应变类似）
应用举例：在分析建筑物在地震中的晃动时，通过测量建筑物各部分的位移，利用几何方程就能计算出建筑物的应变，判断是否会出现裂缝等损坏。
-
物理方程（本构关系）
各向同性弹性材料的胡克定律（矩阵形式）：
应用：连接应力与应变，用于材料选型和实验应力计算。
二、边界条件
-
位移边界条件 指定边界位移值：在边界（如固定端。
-
应力边界条件 指定边界面力：在边界 上，
（(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）或自主编程实现精确分析，广泛应用于航空航天、机械、土木等领域的结构设计与优化。由于变分法是弹性力学的重要内容也是有限元的重要基础，有利于理解有限元的构造及编程方法，后续将介绍这部分内容。提醒：个别地方公式识别有问题，感兴趣的读者可以到本人公众号阅读 冬生亦冬生（全网首发，欢迎一起交流学习）
扫一扫
391
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
