当前位置：首页 > news >正文

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

news 2025/9/17 5:42:50

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

核心步骤：

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

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

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

$-\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 解析公式：

$θmax⁡=2arctan⁡(14[KIKII−sgn(KII)(KIKII)2+8])\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);