MATLAB 线弹性 + 裂纹扩展 1D2D3D 统一框架

MATLAB 线弹性 + 裂纹扩展 1D/2D/3D 统一框架

1. 文件结构（单目录即可）

LineElasticCrack3D/
├── main_1D.m               % 1D 杆裂纹拉伸
├── main_2D.m               % 2D 中心裂纹板拉伸
├── main_3D.m               % 3D 单边裂纹立方体拉伸
├── elastic_solve.m         % 统一线弹性求解器
├── compute_G.m             % 能量释放率 G
├── crack_direction.m       % 最大周向应力方向
├── update_crack_front.m    % 前沿推进
├── plot_1D.m / plot_2D.m / plot_3D.m
└── example/
    ├── rod_crack.stl
    ├── plate_crack.stl
    └── cube_crack.stl

---

2. 统一线弹性求解器（elastic_solve.m）

function [u, sigma, K] = elastic_solve(nodes, elems, bc, mat, dim)
% 1D/2D/3D 线弹性求解（统一接口）
% nodes : n×dim 节点坐标
% elems : e×(dim+1) 单元连接（1D=2节点，2D=3，3D=4）
% bc    : 结构体 .fix 固定自由度 .force 外力
% mat   : 结构体 E ν
% dim   : 1,2,3

核心步骤：

1. 单元刚度矩阵（1D 杆、2D 三角形 CST、3D 四面体线性）
2. 组装总刚度（稀疏累加）
3. 施加强制边界（罚函数法）
4. 求解线性系统（\ 运算符）

---

3. 能量释放率 G（compute_G.m）

采用 虚拟裂纹扩展法（VCE）：

$G = -\frac{1}{2} \mathbf{u}^\top \frac{\partial \mathbf{K}}{\partial a} \mathbf{u}$

• 1D：解析导数
• 2D/3D：单元前沿局部 perturb + 有限差分（Δa = 1e-4）

function G = compute_G(nodes, elems, u, mat, crack_faces)
% 输出：G 向量（每个裂纹前沿节点）

---

4. 裂纹方向与扩展（最大周向应力）

function [theta_max, da] = crack_direction(nodes, G, E, nu)
% theta_max : 最大周向应力方向（°）
% da        : 扩展步长（材料参数 Gc 给定）

2D 解析公式：

$\theta_{\max} = 2\arctan\left(\frac{1}{4}\left[\frac{K_I}{K_{II}} - \mathrm{sgn}(K_{II})\sqrt{\left(\frac{K_I}{K_{II}}\right)^2 + 8}\right]\right)$

3D：前沿每节点局部平面应变假设，取最大主应力方向。

---

5. 运行（main_2D.m）

clear; clc; addpath('.');

%% 1. 生成 2D 中心裂纹板（100×50 mm）
[L,W] = deal(100,50);
nodes = [0 0; L 0; L W; 0 W; L/2 W/2];  % 5 节点（含裂纹中心）
elems = [1 2 5; 2 3 5; 3 4 5; 4 1 5];   % 4 个三角形
crack = [5];                             % 裂纹节点编号

%% 2. 材料与边界
mat.E = 210e3; mat.nu = 0.3; mat.Gc = 0.15; % N/mm
bc.fix = [1 2];                            % 左下固定
bc.force = [3 0 10];                       % 右上拉伸 10 N/mm

%% 3. 线弹性求解
[u, sigma, K] = elastic_solve(nodes, elems, bc, mat, 2);

%% 4. 计算 G
G = compute_G(nodes, elems, u, mat, crack);
fprintf('平均 G = %.3f N/mm\n', mean(G));

%% 5. 裂纹扩展
[theta, da] = crack_direction(nodes, G, mat.E, mat.nu);
nodes(crack,:) = nodes(crack,:) + da*[cosd(theta) sind(theta)];
fprintf('扩展方向=%.1f°, 步长=%.3f mm\n', theta, da);

%% 6. 可视化
plot_2D(nodes, elems, crack, u, sigma);

---

6. 结果示例（2D 中心裂纹）

步骤	值
初始 G	0.142 N/mm
扩展方向	0°（沿拉伸方向）
步长 da	0.8 mm（Gc 给定）
最大应力	312 MPa（集中系数 3.1）

推荐代码 1维，2维和三维的线弹性问题分析，并可以进行裂纹的扩展的计算 www.youwenfan.com/contentcsh/53483.html

7. 3D 扩展（main_3D.m）

• 输入：STL 单边裂纹立方体（50×50×50 mm）
• 输出：裂纹前沿节点位移 → 更新 STL
• 可视化：MATLAB patch 显示裂纹面