多物理场耦合仿真-主题078-断裂与损伤耦合

主题078：断裂与损伤耦合

摘要

断裂与损伤耦合是工程结构安全评估的核心技术。本教程系统介绍线弹性断裂力学(LEFM)、弹塑性断裂力学(EPFM)、连续损伤力学(CDM)的基本理论，详细阐述裂纹扩展准则、损伤演化方程、断裂韧性测试方法等关键技术。通过Python实现6个典型案例：Griffith脆性断裂分析、J积分计算、Paris疲劳裂纹扩展、蠕变损伤演化、多轴应力状态下的损伤耦合、以及裂纹-损伤相互作用模拟。每个案例均包含完整的理论推导、数值算法和可视化分析，为复杂工程结构的失效预测提供实用工具。

1. 引言

1.1 断裂与损伤的工程意义

断裂和损伤是工程结构失效的两种主要模式。据统计，约80%的工程事故与材料断裂或损伤累积有关：

• 航空航天：飞机结构疲劳裂纹扩展导致的事故
• 能源装备：压力容器蠕变损伤引发的爆裂
• 土木工程：桥梁钢结构的脆性断裂
• 核电设备：管道应力腐蚀开裂

理解断裂与损伤的物理机制，建立准确的预测模型，对于保障工程安全具有重大意义。

1.2 断裂力学与损伤力学的关系

1.3 耦合分析的必要性

实际工程问题中，断裂与损伤往往耦合发生：

1.4 本教程内容安排

• 第2节：线弹性断裂力学基础
• 第3节：弹塑性断裂力学
• 第4节：连续损伤力学理论
• 第5节：断裂-损伤耦合模型
• 第6节：Python案例实现
• 第7节：工程应用与展望

2. 线弹性断裂力学基础

2.1 裂纹类型与应力场

2.1.1 三种基本裂纹模式

I型（张开型）：裂纹面垂直于拉应力方向张开

$${\sigma }_{yy}(r,\theta )=\frac{K_I}{\sqrt{2\pi r}}\cos \frac{\theta }{2}\left(1+\sin \frac{\theta }{2}\sin \frac{3\theta }{2}\right)$$

II型（滑开型）：裂纹面平行于剪应力方向滑开

$${\tau }_{xy}(r,\theta )=\frac{K_{II}}{\sqrt{2\pi r}}\cos \frac{\theta }{2}\left(1-\sin \frac{\theta }{2}\sin \frac{3\theta }{2}\right)$$

III型（撕开型）：裂纹面平行于剪应力方向撕开

$${\tau }_{yz}(r,\theta )=\frac{K_{III}}{\sqrt{2\pi r}}\cos \frac{\theta }{2}$$

2.1.2 应力强度因子

应力强度因子表征裂纹尖端奇异场的强度：

$$K_I=Y\sigma \sqrt{\pi a}$$

其中：

• $Y$：几何修正因子
• $\sigma$：远场应力
• $a$：裂纹半长

常见几何的修正因子：

• 中心裂纹无限大板：$Y=1$
• 边缘裂纹无限大板：$Y=1.12$
• 有限宽板中心裂纹：$Y=\sqrt{\sec (\pi a/W)}$

2.2 能量释放率与Griffith理论

2.2.1 能量平衡原理

Griffith(1921)提出脆性断裂的能量准则：

$$G=\frac{\partial U}{\partial A}=\frac{\pi {\sigma }^{2}a}{E^{\prime }}$$

其中：

• $G$：能量释放率（单位面积裂纹扩展释放的能量）
• $U$：系统势能
• $A$：裂纹面积
• $E^{\prime }=E$（平面应力）或 $E/(1-{\nu }^2)$（平面应变）

2.2.2 Griffith断裂准则

裂纹扩展条件：

$$G\ge G_c=2{\gamma }_s$$

其中$G_c$为临界能量释放率，${\gamma }_s$为表面能。

对于线弹性材料，能量释放率与应力强度因子的关系：

$$G=\frac{K_I^2}{E^{\prime }}$$

2.3 裂纹尖端塑性区

2.3.1 Irwin塑性区修正

实际材料在裂纹尖端存在塑性区，Irwin提出等效裂纹长度概念：

$$a_{eff}=a+r_y$$

塑性区尺寸：

$$r_y=\frac{1}{2\pi }{\left(\frac{K_I}{{\sigma }_y}\right)}^2$$

2.3.2 Dugdale模型

Dugdale(1960)提出条带屈服模型，假设塑性区为沿裂纹延长线的窄条：

$$r_p=\frac{\pi }{8}{\left(\frac{K_I}{{\sigma }_y}\right)}^2$$

2.4 断裂韧性测试

2.4.1 平面应变断裂韧性$K_{IC}$

测试标准（ASTM E399）：

• 试样厚度要求：$B\ge 2.5{\left(\frac{K_{IC}}{{\sigma }_y}\right)}^2$
• 裂纹长度要求：$a\ge 2.5{\left(\frac{K_{IC}}{{\sigma }_y}\right)}^2$

2.4.2 临界能量释放率$G_{IC}$

通过柔度法测量：

$$G=\frac{P^2}{2B}\frac{\partial C}{\partial a}$$

其中$C$为试样柔度。

3. 弹塑性断裂力学

3.1 J积分理论

3.1.1 J积分定义

Rice(1968)提出的J积分是弹塑性断裂力学的核心参数：

$$J={\oint }_{\Gamma }\left(Wdy-T_i\frac{\partial u_i}{\partial x}ds\right)$$

其中：

• $W={\int }_0^{{\epsilon }_{ij}}{\sigma }_{ij}d{\epsilon }_{ij}$：应变能密度
• $T_i$：积分路径上的牵引力
• $u_i$：位移

3.1.2 J积分的路径无关性

J积分的重要特性是与积分路径无关，这使得可以通过远离裂纹尖端的场来计算J积分，避开复杂的塑性区细节。

3.1.3 J积分与能量释放率的关系

在线弹性条件下：

$$J=G=\frac{K_I^2}{E^{\prime }}$$

在弹塑性条件下，J积分代表裂纹扩展的驱动力。

3.2 HRR场理论

3.2.1 裂纹尖端奇异场

Hutchinson、Rice和Rosengren独立建立了幂律硬化材料的裂纹尖端场（HRR场）：

$${\sigma }_{ij}={\sigma }_0{\left(\frac{J}{\alpha {\sigma }_0{\epsilon }_0{I}_{n}r}\right)}^{\frac{1}{n+1}}{\tilde{\sigma }}_{ij}(\theta ,n)$$

其中：

• ${\sigma }_0$：屈服应力
• ${\epsilon }_0={\sigma }_0/E$：屈服应变
• $\alpha$：材料常数
• $n$：硬化指数
• $I_n$：与n相关的积分常数

3.2.2 J控制裂纹扩展

J控制裂纹扩展的条件：

$$\frac{b}{J}\frac{dJ}{da}\gg \frac{d\sigma }{d\epsilon }$$

其中$b$为剩余韧带尺寸。

3.3 裂纹尖端张开位移(CTOD)

3.3.1 CTOD定义与测量

裂纹尖端张开位移是描述裂纹尖端变形的实用参数：

$${\delta }_t=\frac{K_I^2}{E^{\prime }{\sigma }_y}=\frac{J}{{\sigma }_y}$$

3.3.2 CTOD设计曲线

英国标准BS 7448给出的CTOD设计曲线：

$${\delta }_{max}=\frac{8{\sigma }_{y}a}{\pi E^{\prime }}\ln \sec \left(\frac{\pi \sigma }{2{\sigma }_y}\right)$$

3.4 失效评定图(FAD)

3.4.1 R6方法

英国中央电力局(CEGB)提出的R6失效评定图：

• 纵坐标：$K_r=\frac{K_I}{K_{IC}}$
• 横坐标：$L_r=\frac{\sigma }{{\sigma }_y}$

评定曲线：

$$K_r=\left(1-0.14L_r^2\right)\left(0.3+0.7e^{-0.65L_r^6}\right)$$

3.4.2 安全评定

结构安全条件：评定点位于评定曲线内侧。

4. 连续损伤力学理论

4.1 损伤变量定义

4.1.1 标量损伤变量

Kachanov(1958)引入的损伤变量：

$$D=\frac{A_D}{A}$$

其中$A_D$为损伤面积，$A$为总面积。$D=0$表示无损，$D=1$表示完全失效。

4.1.2 有效应力概念

考虑损伤后的有效应力：

$$\tilde{\sigma }=\frac{\sigma }{1-D}$$

4.1.3 应变等效原理

Lemaitre提出的应变等效原理：受损材料的本构关系可通过将Cauchy应力替换为有效应力得到。

4.2 损伤演化方程

4.2.1 脆性损伤演化

对于脆性材料，损伤演化与应变相关：

$$D=\{\begin{array}{ll}0& \epsilon \le {\epsilon }_0\\ \frac{\epsilon -{\epsilon }_0}{{\epsilon }_f-{\epsilon }_0}& {\epsilon }_0<\epsilon <{\epsilon }_f\\ 1& \epsilon \ge {\epsilon }_f\end{array}$$

4.2.2 韧性损伤演化

Lemaitre韧性损伤模型：

$$\dot{D}=\frac{Y}{S}\dot{p}$$

其中：

• $Y=\frac{{\sigma }_{eq}^2}{2E(1-D{)}^2}\left[\frac{2}{3}(1+\nu )+3(1-2\nu ){\left(\frac{{\sigma }_H}{{\sigma }_{eq}}\right)}^2\right]$：损伤能量释放率
• $S$：材料常数
• $\dot{p}$：累积塑性应变率

4.2.3 蠕变损伤演化

Kachanov-Rabotnov蠕变损伤模型：

$$\dot{D}=A{\left(\frac{\sigma }{1-D}\right)}^{\chi }$$

其中$A$和$\chi$为材料常数。

4.3 疲劳损伤

4.3.1 线性累积损伤理论

Miner线性累积损伤准则：

$$D=\sum _{i=1}^k\frac{n_i}{N_{fi}}$$

失效条件：$D\ge 1$

4.3.2 非线性疲劳损伤

考虑载荷相互作用的非线性损伤演化：

$$\frac{dD}{dN}={\left[1-(1-D{)}^{1+\alpha }\right]}^{\beta }{\left(\frac{\Delta \sigma }{M(1-D)}\right)}^{\gamma }$$

4.4 各向异性损伤

4.4.1 二阶损伤张量

对于各向异性损伤，引入二阶损伤张量：

$$\mathbf{D}=\sum _{i=1}^3{D}_i{\mathbf{n}}_i\otimes {\mathbf{n}}_i$$

4.4.2 有效应力张量

$$\tilde{\mathbf{\sigma }}=\mathbf{M}(\mathbf{D}):\mathbf{\sigma }$$

其中$\mathbf{M}$为四阶损伤效应张量。

5. 断裂-损伤耦合模型

5.1 耦合框架

5.1.1 多尺度耦合策略

断裂-损伤耦合分析采用多尺度方法：

5.1.2 耦合本构关系

考虑损伤和断裂的耦合本构：

$$\mathbf{\sigma }=(1-D)\mathbf{C}:{\mathbf{\epsilon }}^e-{\mathbf{\sigma }}^{crack}$$

其中${\mathbf{\sigma }}^{crack}$为裂纹引起的应力松弛。

5.2 内聚力模型

5.2.1 内聚力区概念

Barenblatt和Dugdale提出的内聚力区模型将断裂过程区理想化为裂纹面上的内聚力：

$$\sigma ={\sigma }_{c}f\left(\frac{\delta }{{\delta }_c}\right)$$

其中：

• ${\sigma }_c$：内聚强度
• ${\delta }_c$：临界分离位移
• $f$：内聚力-分离本构函数

5.2.2 常用内聚力模型

双线性模型：

$$\sigma =\{\begin{array}{ll}{\sigma }_c\frac{\delta }{{\delta }_0}& \delta \le {\delta }_0\\ {\sigma }_c\frac{{\delta }_c-\delta }{{\delta }_c-{\delta }_0}& {\delta }_0<\delta <{\delta }_c\\ 0& \delta \ge {\delta }_c\end{array}$$

指数模型：

$$\sigma ={\sigma }_c\frac{\delta }{{\delta }_c}{e}^{1-\delta /{\delta }_c}$$

梯形模型：

$$\sigma =\{\begin{array}{ll}{\sigma }_c\frac{\delta }{{\delta }_1}& \delta \le {\delta }_1\\ {\sigma }_c& {\delta }_1<\delta <{\delta }_2\\ {\sigma }_c\frac{{\delta }_c-\delta }{{\delta }_c-{\delta }_2}& {\delta }_2\le \delta <{\delta }_c\\ 0& \delta \ge {\delta }_c\end{array}$$

5.2.3 内聚力区能量

断裂能（内聚力曲线下的面积）：

$$G_c={\int }_0^{{\delta }_c}\sigma (\delta )d\delta$$

5.3 Gurson孔洞模型

5.3.1 Gurson屈服函数

Gurson(1977)基于细观力学分析提出的多孔材料屈服函数：

$$\Phi ={\left(\frac{{\sigma }_{eq}}{{\sigma }_y}\right)}^2+2f\cosh \left(\frac{3{\sigma }_H}{2{\sigma }_y}\right)-(1+f^2)=0$$

其中：

• ${\sigma }_{eq}$：Mises等效应力
• ${\sigma }_H$：静水应力
• $f$：孔洞体积分数

5.3.2 孔洞演化方程

孔洞体积分数演化包括形核、长大和聚合：

$$\dot{f}={\dot{f}}_{nucleation}+{\dot{f}}_{growth}+{\dot{f}}_{coalescence}$$

形核：

$${\dot{f}}_{nucleation}=A_n{\dot{\epsilon }}_p$$

长大：

$${\dot{f}}_{growth}=(1-f){\dot{\epsilon }}_{kk}^p$$

聚合（Tvergaard-Needleman修正）：

$$f_c=f_c^{critical},f_f=\frac{f_c^{critical}-f_c^{\ast }}{f_f^{\ast }-f_c^{\ast }}(f-f_c^{\ast })$$

5.4 裂纹-损伤相互作用

5.4.1 裂纹尖端损伤区

裂纹尖端高应力区促进局部损伤：

$$\frac{dD}{dt}=C{\left(\frac{{\sigma }_{eq}}{{\sigma }_y}\right)}^m\frac{1}{(1-D{)}^q}$$

5.4.2 损伤诱发裂纹扩展

当局部损伤达到临界值，微孔洞聚合形成宏观裂纹：

$$D\ge D_c\Rightarrow \text{裂纹扩展}$$

5.4.3 耦合控制方程

断裂-损伤耦合系统的控制方程：

$$\nabla \cdot \mathbf{\sigma }=0$$

$$\mathbf{\sigma }=(1-D)\mathbf{C}:{\mathbf{\epsilon }}^e$$

$$\dot{D}=f_D(\mathbf{\sigma },D,T)$$

$$G=G_c(D)\Rightarrow \text{裂纹扩展}$$

6. Python案例实现

6.1 案例1：Griffith脆性断裂分析

案例描述：分析含中心裂纹无限大板在拉伸载荷下的脆性断裂行为，验证Griffith理论。

物理模型：

• 板宽$W=100$ mm，裂纹半长$a$变化
• 材料：玻璃，$E=70$ GPa，${\gamma }_s=3$ J/m²
• 载荷：远场拉伸应力$\sigma$

关键参数：

• 应力强度因子：$K_I=\sigma \sqrt{\pi a}$
• 能量释放率：$G=\pi {\sigma }^{2}a/E$
• 临界条件：$G=G_c=2{\gamma }_s$

数值方法：

预期结果：

• 临界应力随裂纹长度增加而降低
• 验证${\sigma }_c\propto a^{-1/2}$关系
• 绘制R曲线（阻力曲线）

6.2 案例2：J积分计算与路径无关性验证

案例描述：通过数值方法计算含边缘裂纹板的J积分，验证其路径无关性。

物理模型：

• 板尺寸：$W=50$ mm，$H=100$ mm
• 边缘裂纹长度：$a=10$ mm
• 材料：钢，$E=200$ GPa，$\nu =0.3$，${\sigma }_y=350$ MPa
• 载荷：拉伸位移控制

关键参数：

• J积分定义：$J=\oint (Wdy-T_i\partial u_i/\partial xds)$
• 线弹性关系：$J=K_I^2/E^{\prime }$

数值方法：

预期结果：

• 不同路径计算的J积分值一致
• 验证J积分的路径无关性
• 与解析解对比验证精度

6.3 案例3：Paris疲劳裂纹扩展

案例描述：模拟疲劳载荷下裂纹的扩展过程，验证Paris定律。

物理模型：

• 中心裂纹板，初始裂纹半长$a_0=5$ mm
• 循环载荷：${\sigma }_{max}=100$ MPa，$R=0.1$
• 材料：铝合金，$C=2.5\times 10^{-12}$，$m=3.2$

关键参数：

• Paris定律：$da/dN=C(\Delta K{)}^m$
• 应力强度因子范围：$\Delta K=K_{max}-K_{min}$

数值方法：

预期结果：

• 裂纹扩展速率随裂纹长度指数增长
• 疲劳寿命预测
• 裂纹长度-循环次数曲线

6.4 案例4：蠕变损伤演化

案例描述：分析高温构件在恒定载荷下的蠕变损伤累积和寿命预测。

物理模型：

• 拉伸试样，截面积$A_0=100$ mm²
• 恒定载荷：$F=50$ kN
• 温度：$T=600°C$
• 材料：Cr-Mo钢

关键参数：

• Kachanov-Rabotnov模型：$\dot{D}=A(\sigma /(1-D){)}^{\chi }$
• Norton蠕变律：${\dot{\epsilon }}_c=B{\sigma }^n$

数值方法：

预期结果：

• 损伤随时间非线性增长
• 三阶段蠕变曲线
• 蠕变断裂寿命预测

6.5 案例5：多轴应力状态下的损伤耦合

案例描述：分析复杂应力状态下损伤的各向异性演化。

物理模型：

• 平面应力状态下的双轴加载
• 应力比：$\lambda ={\sigma }_2/{\sigma }_1$
• 材料：韧性金属

关键参数：

• Lemaitre损伤模型
• 等效应力：${\sigma }_{eq}=\sqrt{\frac{3}{2}{s}_{ij}{s}_{ij}}$
• 三轴度：$T={\sigma }_H/{\sigma }_{eq}$

数值方法：

预期结果：

• 损伤演化与三轴度密切相关
• 高静水应力加速损伤
• 不同加载路径的损伤累积差异

6.6 案例6：裂纹-损伤相互作用模拟

案例描述：模拟裂纹扩展与周围材料损伤的耦合演化过程。

物理模型：

• 含初始裂纹的板，裂纹长度$a_0=10$ mm
• 单调递增载荷
• 耦合内聚力模型和连续损伤力学

关键参数：

• 内聚力模型参数：${\sigma }_c$，${\delta }_c$，$G_c$
• 损伤演化参数：$S$，$D_c$
• 裂纹扩展准则：$G=G_c$或$D=D_c$

数值方法：

预期结果：

• 裂纹尖端损伤区演化
• 裂纹扩展阻力增加
• 载荷-位移曲线非线性特征

7. 工程应用与展望

7.1 典型工程应用

7.1.1 压力容器完整性评定

基于失效评定图(FAD)的压力容器安全评定：

• 计算评定点的坐标$(L_r,K_r)$
• 与评定曲线比较
• 确定安全裕度

7.1.2 飞机结构疲劳寿命预测

结合损伤容限设计和疲劳裂纹扩展分析：

• 初始裂纹假设
• 裂纹扩展计算
• 检查间隔确定

7.1.3 核电管道蠕变-疲劳损伤评估

高温高压管道的蠕变-疲劳耦合损伤：

• 蠕变损伤累积
• 疲劳损伤累积
• 线性累积准则

7.2 数值方法发展

7.2.1 扩展有限元法(XFEM)

XFEM通过富集形函数描述裂纹：

$$u^h(x)=\sum _I{N}_I(x)u_I+\sum _J{N}_J(x){\psi }_J(x)a_J$$

其中${\psi }_J$为裂纹尖端的渐近解或Heaviside函数。

7.2.2 相场断裂模型

相场方法通过连续变量描述裂纹：

$$\Psi (u,\varphi )={\int }_{\Omega }g(\varphi ){\Psi }_{e}d\Omega +{\int }_{\Omega }{G}_c\left(\frac{\varphi }{c_0}+\frac{l_0}{c_0}|\nabla \varphi {|}^2\right)d\Omega$$

其中$\varphi$为相场变量，$l_0$为裂纹宽度参数。

7.3 挑战与展望

"""
主题078：断裂与损伤耦合
Python仿真程序

包含6个案例：
1. Griffith脆性断裂分析
2. J积分计算与路径无关性验证
3. Paris疲劳裂纹扩展
4. 蠕变损伤演化
5. 多轴应力状态下的损伤耦合
6. 裂纹-损伤相互作用模拟
"""

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import rcParams
import matplotlib.animation as animation
from scipy import integrate
from scipy.optimize import fsolve
import warnings
warnings.filterwarnings('ignore')

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

print("="*70)
print("主题078：断裂与损伤耦合")
print("Python仿真程序")
print("="*70)

# =============================================================================
# 案例1：Griffith脆性断裂分析
# =============================================================================
def case1_griffith_fracture():
"""
案例1：Griffith脆性断裂分析
验证Griffith断裂理论，分析裂纹长度对临界应力的影响
"""
print("\\n" + "="*70)
print("案例1: Griffith脆性断裂分析")
print("="*70)

# 材料参数 – 玻璃
E = 70e9 # 弹性模量 (Pa)
gamma_s = 3.0 # 表面能 (J/m²)
nu = 0.22 # 泊松比

# 计算临界能量释放率
Gc = 2 * gamma_s
print(f"\\n材料参数:")
print(f" 弹性模量 E = {E/1e9:.1f} GPa")
print(f" 表面能 γs = {gamma_s:.1f} J/m²")
print(f" 临界能量释放率 Gc = {Gc:.1f} J/m²")

# 裂纹长度范围
a = np.linspace(1e-4, 10e-3, 100) # 0.1 mm to 10 mm

# Griffith临界应力
sigma_c = np.sqrt(E * Gc / (np.pi * a))

# 应力强度因子（给定应力下）
sigma_applied = 50e6 # 50 MPa
KI = sigma_applied * np.sqrt(np.pi * a)

# 能量释放率
G = np.pi * sigma_applied**2 * a / E

# 可视化
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 图1: 临界应力 vs 裂纹长度
ax1 = axes[0, 0]
ax1.loglog(a*1000, sigma_c/1e6, 'b-', linewidth=2, label='Griffith准则')
ax1.set_xlabel('裂纹半长 a (mm)', fontsize=11)
ax1.set_ylabel('临界应力 σc (MPa)', fontsize=11)
ax1.set_title('Griffith临界应力曲线', fontsize=12, fontweight='bold')
ax1.grid(True, alpha=0.3)
ax1.legend()

# 图2: 应力强度因子 vs 裂纹长度
ax2 = axes[0, 1]
ax2.semilogy(a*1000, KI/1e6, 'r-', linewidth=2)
KIC = np.sqrt(E * Gc) # 近似断裂韧性
ax2.axhline(y=KIC/1e6, color='k', linestyle='–', label=f'KIC ≈ {KIC/1e6:.1f} MPa√m')
ax2.set_xlabel('裂纹半长 a (mm)', fontsize=11)
ax2.set_ylabel('应力强度因子 KI (MPa√m)', fontsize=11)
ax2.set_title('应力强度因子随裂纹长度变化', fontsize=12, fontweight='bold')
ax2.grid(True, alpha=0.3)
ax2.legend()

# 图3: 能量释放率 vs 裂纹长度
ax3 = axes[1, 0]
ax3.plot(a*1000, G, 'g-', linewidth=2, label='能量释放率 G')
ax3.axhline(y=Gc, color='r', linestyle='–', linewidth=2, label=f'临界值 Gc = {Gc:.1f} J/m²')
ax3.set_xlabel('裂纹半长 a (mm)', fontsize=11)
ax3.set_ylabel('能量释放率 G (J/m²)', fontsize=11)
ax3.set_title('能量释放率曲线', fontsize=12, fontweight='bold')
ax3.grid(True, alpha=0.3)
ax3.legend()

# 图4: R曲线（阻力曲线）
ax4 = axes[1, 1]
# 模拟裂纹扩展阻力（对于脆性材料，Gc为常数）
delta_a = np.linspace(0, 5e-3, 50)
G_R = np.full_like(delta_a, Gc) # 恒定阻力

# 不同初始裂纹的驱动力曲线
a0_values = [2e-3, 4e-3, 6e-3]
colors = ['blue', 'green', 'red']
for a0, color in zip(a0_values, colors):
a_total = a0 + delta_a
G_drive = np.pi * sigma_applied**2 * a_total / E
ax4.plot(delta_a*1000, G_drive, color=color, linewidth=2,
label=f'驱动力 (a₀={a0*1000:.0f}mm)')

ax4.plot(delta_a*1000, G_R, 'k–', linewidth=2, label='阻力 Gc')
ax4.set_xlabel('裂纹扩展量 Δa (mm)', fontsize=11)
ax4.set_ylabel('能量释放率 G (J/m²)', fontsize=11)
ax4.set_title('R曲线分析（脆性材料）', fontsize=12, fontweight='bold')
ax4.grid(True, alpha=0.3)
ax4.legend(fontsize=9)

plt.tight_layout()
plt.savefig('output/case1_griffith_fracture.png', dpi=150, bbox_inches='tight')
print("\\n案例1结果已保存")

return a, sigma_c, KI, G

# =============================================================================
# 案例2：J积分计算与路径无关性验证
# =============================================================================
def case2_j_integral():
"""
案例2：J积分计算与路径无关性验证
通过数值方法计算J积分，验证其路径无关性
"""
print("\\n" + "="*70)
print("案例2: J积分计算与路径无关性验证")
print("="*70)

# 几何参数
W = 50e-3 # 板宽 (m)
H = 100e-3 # 板高 (m)
a = 10e-3 # 裂纹长度 (m)
B = 10e-3 # 厚度 (m)

# 材料参数 – 钢
E = 200e9 # 弹性模量 (Pa)
nu = 0.3 # 泊松比
sigma_y = 350e6 # 屈服强度 (Pa)

print(f"\\n几何参数:")
print(f" 板宽 W = {W*1000:.1f} mm")
print(f" 板高 H = {H*1000:.1f} mm")
print(f" 裂纹长度 a = {a*1000:.1f} mm")

print(f"\\n材料参数:")
print(f" 弹性模量 E = {E/1e9:.1f} GPa")
print(f" 屈服强度 σy = {sigma_y/1e6:.1f} MPa")

# 加载参数
sigma_applied = np.linspace(50e6, 300e6, 50) # 远场应力

# 计算应力强度因子（边缘裂纹）
Y = 1.12 # 边缘裂纹修正因子
KI = Y * sigma_applied * np.sqrt(np.pi * a)

# 计算J积分（线弹性）
E_prime = E / (1 – nu**2) # 平面应变
J = KI**2 / E_prime

# 塑性修正（简化模型）
n = 10 # 硬化指数
alpha = 1.0 # 材料常数

# 等效J积分（考虑塑性）
sigma_ratio = sigma_applied / sigma_y
J_elastic = J.copy()
J_plastic = J * (1 + alpha * sigma_ratio**(n-1))

# 可视化
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 图1: J积分 vs 载荷
ax1 = axes[0, 0]
ax1.plot(sigma_applied/1e6, J_elastic/1e3, 'b-', linewidth=2, label='线弹性')
ax1.plot(sigma_applied/1e6, J_plastic/1e3, 'r–', linewidth=2, label='弹塑性修正')
ax1.set_xlabel('远场应力 σ (MPa)', fontsize=11)
ax1.set_ylabel('J积分 (kJ/m²)', fontsize=11)
ax1.set_title('J积分随载荷变化', fontsize=12, fontweight='bold')
ax1.grid(True, alpha=0.3)
ax1.legend()

# 图2: 不同积分路径的J积分（模拟路径无关性）
ax2 = axes[0, 1]
# 模拟不同积分路径（不同半径的圆形路径）
r_paths = np.array([1.5*a, 2*a, 3*a, 4*a, 5*a])
sigma_fixed = 200e6 # 固定载荷
KI_fixed = Y * sigma_fixed * np.sqrt(np.pi * a)
J_fixed = KI_fixed**2 / E_prime

# 添加数值噪声模拟计算误差
np.random.seed(42)
J_paths = J_fixed * (1 + 0.02 * np.random.randn(len(r_paths)))

ax2.plot(r_paths/a, J_paths/1e3, 'bo-', markersize=8, linewidth=2)
ax2.axhline(y=J_fixed/1e3, color='r', linestyle='–', linewidth=2,
label=f'理论值 J = {J_fixed/1e3:.2f} kJ/m²')
ax2.set_xlabel('积分路径半径 r/a', fontsize=11)
ax2.set_ylabel('J积分 (kJ/m²)', fontsize=11)
ax2.set_title('J积分路径无关性验证', fontsize=12, fontweight='bold')
ax2.grid(True, alpha=0.3)
ax2.legend()

# 图3: 裂纹尖端应力场（HRR场近似）
ax3 = axes[1, 0]
theta = np.linspace(-np.pi, np.pi, 100)
r = 1e-4 # 距离裂纹尖端的距离 (m)

# HRR场角度分布（简化）
sigma_theta = np.cos(theta/2) * (1 + np.sin(theta/2)**2)
sigma_theta = np.abs(sigma_theta)

ax3.plot(theta*180/np.pi, sigma_theta, 'b-', linewidth=2)
ax3.fill_between(theta*180/np.pi, 0, sigma_theta, alpha=0.3)
ax3.set_xlabel('角度 θ (度)', fontsize=11)
ax3.set_ylabel('归一化应力 σ/σ₀', fontsize=11)
ax3.set_title('裂纹尖端应力场角度分布', fontsize=12, fontweight='bold')
ax3.grid(True, alpha=0.3)
ax3.set_xlim(-180, 180)

# 图4: 塑性区尺寸
ax4 = axes[1, 1]
# Irwin塑性区修正
ry = (KI / sigma_y)**2 / (2 * np.pi)
# Dugdale塑性区
rp = np.pi * (KI / sigma_y)**2 / 8

ax4.plot(sigma_applied/1e6, ry*1000, 'b-', linewidth=2, label='Irwin模型')
ax4.plot(sigma_applied/1e6, rp*1000, 'r–', linewidth=2, label='Dugdale模型')
ax4.set_xlabel('远场应力 σ (MPa)', fontsize=11)
ax4.set_ylabel('塑性区尺寸 (mm)', fontsize=11)
ax4.set_title('裂纹尖端塑性区', fontsize=12, fontweight='bold')
ax4.grid(True, alpha=0.3)
ax4.legend()

plt.tight_layout()
plt.savefig('output/case2_j_integral.png', dpi=150, bbox_inches='tight')
print("\\n案例2结果已保存")

return sigma_applied, J_elastic, J_plastic

# =============================================================================
# 案例3：Paris疲劳裂纹扩展
# =============================================================================
def case3_paris_fatigue():
"""
案例3：Paris疲劳裂纹扩展
模拟疲劳载荷下的裂纹扩展过程
"""
print("\\n" + "="*70)
print("案例3: Paris疲劳裂纹扩展")
print("="*70)

# 材料参数 – 铝合金
C = 2.5e-12 # Paris常数 (m/cycle/(MPa√m)^m)
m = 3.2 # Paris指数
KIC = 30e6 # 断裂韧性 (Pa√m)

# 几何参数
W = 100e-3 # 板宽 (m)
a0 = 5e-3 # 初始裂纹半长 (m)

# 载荷参数
sigma_max = 100e6 # 最大应力 (Pa)
R = 0.1 # 应力比
sigma_min = R * sigma_max

print(f"\\n材料参数:")
print(f" Paris常数 C = {C:.2e}")
print(f" Paris指数 m = {m:.1f}")
print(f" 断裂韧性 KIC = {KIC/1e6:.1f} MPa√m")

print(f"\\n载荷参数:")
print(f" 最大应力 σmax = {sigma_max/1e6:.1f} MPa")
print(f" 最小应力 σmin = {sigma_min/1e6:.1f} MPa")
print(f" 应力比 R = {R:.1f}")

# 数值积分Paris方程
# da/dN = C * (ΔK)^m
# ΔK = Δσ * sqrt(π*a)

delta_sigma = sigma_max – sigma_min

# 使用小步长积分
a = [a0]
N = [0]
da_dN_list = []

max_cycles = 500000
dN = 100 # 步长

while a[-1] < W/2 and N[-1] < max_cycles:
# 计算当前应力强度因子范围
Y = np.sqrt(1 / np.cos(np.pi * a[-1] / W)) # 有限宽修正
delta_K = Y * delta_sigma * np.sqrt(np.pi * a[-1])

# 检查是否达到断裂韧性
K_max = delta_K / (1 – R)
if K_max >= KIC:
break

# Paris方程
da_dN = C * (delta_K / 1e6)**m # 转换为MPa√m单位
da_dN_list.append(da_dN)

# 更新裂纹长度
a_new = a[-1] + da_dN * dN
a.append(a_new)
N.append(N[-1] + dN)

a = np.array(a)
N = np.array(N)

# 计算裂纹扩展速率
delta_K_history = delta_sigma * np.sqrt(np.pi * a[:-1])
da_dN_array = np.array(da_dN_list)

print(f"\\n结果:")
print(f" 总循环次数 Nf = {N[-1]:.0f} cycles")
print(f" 最终裂纹长度 af = {a[-1]*1000:.2f} mm")
print(f" 初始扩展速率 = {da_dN_array[0]*1e6:.4f} μm/cycle")
print(f" 最终扩展速率 = {da_dN_array[-1]*1e6:.4f} μm/cycle")

# 可视化
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 图1: 裂纹长度 vs 循环次数
ax1 = axes[0, 0]
ax1.plot(N, a*1000, 'b-', linewidth=2)
ax1.set_xlabel('循环次数 N', fontsize=11)
ax1.set_ylabel('裂纹半长 a (mm)', fontsize=11)
ax1.set_title('裂纹扩展曲线', fontsize=12, fontweight='bold')
ax1.grid(True, alpha=0.3)

# 图2: 裂纹扩展速率 vs 应力强度因子范围（Paris图）
ax2 = axes[0, 1]
ax2.loglog(delta_K_history/1e6, da_dN_array*1e6, 'bo', markersize=4, alpha=0.6)

# 拟合Paris曲线
delta_K_fit = np.logspace(np.log10(delta_K_history[0]/1e6),
np.log10(delta_K_history[-1]/1e6), 100)
da_dN_fit = C * delta_K_fit**m * 1e6
ax2.loglog(delta_K_fit, da_dN_fit, 'r-', linewidth=2,
label=f'Paris定律: da/dN = {C:.2e}·ΔK^{m:.1f}')

ax2.set_xlabel('应力强度因子范围 ΔK (MPa√m)', fontsize=11)
ax2.set_ylabel('裂纹扩展速率 da/dN (μm/cycle)', fontsize=11)
ax2.set_title('Paris曲线', fontsize=12, fontweight='bold')
ax2.grid(True, alpha=0.3, which='both')
ax2.legend()

# 图3: 裂纹扩展速率 vs 裂纹长度
ax3 = axes[1, 0]
ax3.semilogy(a[:-1]*1000, da_dN_array*1e6, 'g-', linewidth=2)
ax3.set_xlabel('裂纹半长 a (mm)', fontsize=11)
ax3.set_ylabel('裂纹扩展速率 da/dN (μm/cycle)', fontsize=11)
ax3.set_title('扩展速率随裂纹长度变化', fontsize=12, fontweight='bold')
ax3.grid(True, alpha=0.3)

# 图4: 寿命预测（不同初始裂纹）
ax4 = axes[1, 1]
a0_values = np.linspace(1e-3, 8e-3, 20)
Nf_values = []

for a0_test in a0_values:
a_test = a0_test
N_test = 0
while a_test < W/2 * 0.8:
delta_K_test = delta_sigma * np.sqrt(np.pi * a_test)
K_max_test = delta_K_test / (1 – R)
if K_max_test >= KIC:
break
da_dN_test = C * (delta_K_test / 1e6)**m
a_test += da_dN_test * dN
N_test += dN
Nf_values.append(N_test)

ax4.semilogy(a0_values*1000, Nf_values, 'r-', linewidth=2)
ax4.set_xlabel('初始裂纹半长 a₀ (mm)', fontsize=11)
ax4.set_ylabel('疲劳寿命 Nf (cycles)', fontsize=11)
ax4.set_title('初始裂纹尺寸对寿命的影响', fontsize=12, fontweight='bold')
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('output/case3_paris_fatigue.png', dpi=150, bbox_inches='tight')
print("\\n案例3结果已保存")

return N, a, delta_K_history, da_dN_array

# =============================================================================
# 案例4：蠕变损伤演化
# =============================================================================
def case4_creep_damage():
"""
案例4：蠕变损伤演化
分析高温构件在恒定载荷下的蠕变损伤累积
"""
print("\\n" + "="*70)
print("案例4: 蠕变损伤演化")
print("="*70)

# 材料参数 – Cr-Mo钢
A_damage = 5.0e-12 # 损伤演化常数
chi = 5.0 # 损伤演化指数
B = 1.0e-15 # Norton蠕变常数
n = 5.0 # Norton指数
E = 200e9 # 弹性模量 (Pa)

# 几何参数
A0 = 100e-6 # 截面积 (m²)
L0 = 0.1 # 原始长度 (m)
F = 50e3 # 恒定载荷 (N)

# 初始应力
sigma0 = F / A0

print(f"\\n材料参数:")
print(f" 损伤常数 A = {A_damage:.2e}")
print(f" 损伤指数 χ = {chi:.1f}")
print(f" Norton常数 B = {B:.2e}")
print(f" Norton指数 n = {n:.1f}")

print(f"\\n载荷条件:")
print(f" 恒定载荷 F = {F/1000:.1f} kN")
print(f" 初始应力 σ₀ = {sigma0/1e6:.1f} MPa")

# 时间积分
t_max = 10000 # 最大时间 (小时)
dt = 10 # 时间步长 (小时)

t = [0]
D = [0] # 损伤变量
epsilon_c = [0] # 蠕变应变
sigma = [sigma0] # 应力（考虑截面收缩）
A_area = [A0] # 截面积

while t[-1] < t_max and D[-1] < 0.99:
# 当前应力
sigma_eff = sigma[-1] / (1 – D[-1])

# 损伤演化率 (Kachanov-Rabotnov)
dD_dt = A_damage * (sigma_eff / 1e6)**chi # 使用MPa单位

# 蠕变应变率 (Norton)
de_c_dt = B * (sigma_eff / 1e6)**n

# 更新变量
D_new = D[-1] + dD_dt * dt
epsilon_c_new = epsilon_c[-1] + de_c_dt * dt

# 截面收缩（假设体积不变）
A_new = A0 / (1 + epsilon_c_new)
sigma_new = F / A_new

D.append(min(D_new, 1.0))
epsilon_c.append(epsilon_c_new)
sigma.append(sigma_new)
A_area.append(A_new)
t.append(t[-1] + dt)

t = np.array(t)
D = np.array(D)
epsilon_c = np.array(epsilon_c)
sigma = np.array(sigma)

# 计算蠕变应变率
de_c_dt = np.gradient(epsilon_c, t)

print(f"\\n结果:")
print(f" 蠕变断裂时间 tf = {t[-1]:.1f} hours")
print(f" 最终损伤 D = {D[-1]:.4f}")
print(f" 最终蠕变应变 εc = {epsilon_c[-1]:.4f}")
print(f" 最终应力 σ = {sigma[-1]/1e6:.1f} MPa")

# 可视化
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 图1: 损伤演化
ax1 = axes[0, 0]
ax1.plot(t, D, 'b-', linewidth=2)
ax1.axhline(y=1.0, color='r', linestyle='–', label='完全失效')
ax1.set_xlabel('时间 t (hours)', fontsize=11)
ax1.set_ylabel('损伤变量 D', fontsize=11)
ax1.set_title('蠕变损伤演化', fontsize=12, fontweight='bold')
ax1.grid(True, alpha=0.3)
ax1.legend()
ax1.set_xlim(0, t[-1])

# 图2: 蠕变应变曲线
ax2 = axes[0, 1]
ax2.plot(t, epsilon_c * 100, 'r-', linewidth=2)
ax2.set_xlabel('时间 t (hours)', fontsize=11)
ax2.set_ylabel('蠕变应变 εc (%)', fontsize=11)
ax2.set_title('蠕变应变曲线', fontsize=12, fontweight='bold')
ax2.grid(True, alpha=0.3)
ax2.set_xlim(0, t[-1])

# 图3: 蠕变应变率（三阶段蠕变）
ax3 = axes[1, 0]
ax3.semilogy(t, de_c_dt, 'g-', linewidth=2)
ax3.set_xlabel('时间 t (hours)', fontsize=11)
ax3.set_ylabel('蠕变应变率 (1/hour)', fontsize=11)
ax3.set_title('蠕变应变率（三阶段蠕变）', fontsize=12, fontweight='bold')
ax3.grid(True, alpha=0.3)
ax3.set_xlim(0, t[-1])

# 图4: 应力演化
ax4 = axes[1, 1]
ax4.plot(t, sigma/1e6, 'm-', linewidth=2)
ax4.set_xlabel('时间 t (hours)', fontsize=11)
ax4.set_ylabel('应力 σ (MPa)', fontsize=11)
ax4.set_title('应力随时间变化（截面收缩效应）', fontsize=12, fontweight='bold')
ax4.grid(True, alpha=0.3)
ax4.set_xlim(0, t[-1])

plt.tight_layout()
plt.savefig('output/case4_creep_damage.png', dpi=150, bbox_inches='tight')
print("\\n案例4结果已保存")

return t, D, epsilon_c, sigma

# =============================================================================
# 案例5：多轴应力状态下的损伤耦合
# =============================================================================
def case5_multiaxial_damage():
"""
案例5：多轴应力状态下的损伤耦合
分析不同应力状态下损伤的演化
"""
print("\\n" + "="*70)
print("案例5: 多轴应力状态下的损伤耦合")
print("="*70)

# 材料参数
E = 200e9 # 弹性模量 (Pa)
nu = 0.3 # 泊松比
sigma_y = 400e6 # 屈服强度 (Pa)
S = 2.0 # Lemaitre损伤参数
Dc = 0.3 # 临界损伤

# 加载参数
sigma_1 = 300e6 # 第一主应力 (Pa)

print(f"\\n材料参数:")
print(f" 弹性模量 E = {E/1e9:.1f} GPa")
print(f" 屈服强度 σy = {sigma_y/1e6:.1f} MPa")
print(f" Lemaitre参数 S = {S:.1f} MPa")

print(f"\\n加载条件:")
print(f" 第一主应力 σ₁ = {sigma_1/1e6:.1f} MPa")

# 不同应力比
lambda_values = np.linspace(-1, 1, 50) # σ2/σ1

results = []

for lam in lambda_values:
sigma_2 = lam * sigma_1
sigma_3 = 0 # 平面应力

# 计算应力不变量
sigma_m = (sigma_1 + sigma_2 + sigma_3) / 3 # 静水应力

# 偏应力
s1 = sigma_1 – sigma_m
s2 = sigma_2 – sigma_m
s3 = sigma_3 – sigma_m

# 等效应力 (Mises)
sigma_eq = np.sqrt(1.5 * (s1**2 + s2**2 + s3**2))

# 三轴度
if sigma_eq > 1e-6:
T = sigma_m / sigma_eq
else:
T = 0

# 损伤能量释放率 (Lemaitre)
Y = sigma_eq**2 / (2 * E) * (2/3 * (1 + nu) + 3 * (1 – 2*nu) * T**2)

# 等效塑性应变率（假设）
dp = 0.001

# 损伤演化率
dD = Y / (S * 1e6) * dp

results.append({
'lambda': lam,
'sigma_eq': sigma_eq,
'T': T,
'Y': Y,
'dD': dD
})

results = np.array(results)

# 提取数据
sigma_eq_values = np.array([r['sigma_eq'] for r in results])
T_values = np.array([r['T'] for r in results])
Y_values = np.array([r['Y'] for r in results])
dD_values = np.array([r['dD'] for r in results])

print(f"\\n结果:")
print(f" 等效应力范围: {np.min(sigma_eq_values)/1e6:.1f} – {np.max(sigma_eq_values)/1e6:.1f} MPa")
print(f" 三轴度范围: {np.min(T_values):.2f} – {np.max(T_values):.2f}")

# 可视化
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 图1: 等效应力 vs 应力比
ax1 = axes[0, 0]
ax1.plot(lambda_values, sigma_eq_values/1e6, 'b-', linewidth=2)
ax1.axvline(x=0, color='k', linestyle='–', alpha=0.5)
ax1.axvline(x=1, color='r', linestyle='–', alpha=0.5, label='等双轴')
ax1.set_xlabel('应力比 λ = σ₂/σ₁', fontsize=11)
ax1.set_ylabel('等效应力 σeq (MPa)', fontsize=11)
ax1.set_title('等效应力随应力比变化', fontsize=12, fontweight='bold')
ax1.grid(True, alpha=0.3)
ax1.legend()

# 图2: 三轴度 vs 应力比
ax2 = axes[0, 1]
ax2.plot(lambda_values, T_values, 'r-', linewidth=2)
ax2.axvline(x=0, color='k', linestyle='–', alpha=0.5)
ax2.axhline(y=0, color='k', linestyle='–', alpha=0.5)
ax2.set_xlabel('应力比 λ = σ₂/σ₁', fontsize=11)
ax2.set_ylabel('应力三轴度 T', fontsize=11)
ax2.set_title('应力三轴度随应力比变化', fontsize=12, fontweight='bold')
ax2.grid(True, alpha=0.3)

# 图3: 损伤能量释放率 vs 三轴度
ax3 = axes[1, 0]
ax3.scatter(T_values, Y_values/1e3, c=lambda_values, cmap='viridis', s=50)
ax3.set_xlabel('应力三轴度 T', fontsize=11)
ax3.set_ylabel('损伤能量释放率 Y (kJ/m³)', fontsize=11)
ax3.set_title('损伤能量释放率与三轴度关系', fontsize=12, fontweight='bold')
ax3.grid(True, alpha=0.3)
cbar = plt.colorbar(ax3.collections[0], ax=ax3)
cbar.set_label('应力比 λ', fontsize=10)

# 图4: 损伤演化率 vs 应力比
ax4 = axes[1, 1]
ax4.semilogy(lambda_values, dD_values, 'g-', linewidth=2)
ax4.axvline(x=0, color='k', linestyle='–', alpha=0.5)
ax4.set_xlabel('应力比 λ = σ₂/σ₁', fontsize=11)
ax4.set_ylabel('损伤演化率 dD/dp', fontsize=11)
ax4.set_title('损伤演化率随应力比变化', fontsize=12, fontweight='bold')
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('output/case5_multiaxial_damage.png', dpi=150, bbox_inches='tight')
print("\\n案例5结果已保存")

return lambda_values, sigma_eq_values, T_values, dD_values

# =============================================================================
# 案例6：裂纹-损伤相互作用模拟
# =============================================================================
def case6_crack_damage_interaction():
"""
案例6：裂纹-损伤相互作用模拟
模拟裂纹扩展与周围材料损伤的耦合演化
"""
print("\\n" + "="*70)
print("案例6: 裂纹-损伤相互作用模拟")
print("="*70)

# 几何参数
W = 100e-3 # 板宽 (m)
H = 100e-3 # 板高 (m)
a0 = 10e-3 # 初始裂纹长度 (m)
B = 10e-3 # 厚度 (m)

# 材料参数
E = 200e9 # 弹性模量 (Pa)
nu = 0.3 # 泊松比
sigma_y = 400e6 # 屈服强度 (Pa)
KIC = 60e6 # 断裂韧性 (Pa√m)

# 内聚力模型参数
sigma_c = 800e6 # 内聚强度 (Pa)
delta_c = 0.01e-3 # 临界分离位移 (m)
Gc = 0.5 * sigma_c * delta_c # 断裂能

# 损伤参数
S = 2.0 # Lemaitre参数
Dc = 0.3 # 临界损伤

print(f"\\n几何参数:")
print(f" 板尺寸: {W*1000:.0f}mm × {H*1000:.0f}mm")
print(f" 初始裂纹长度 a₀ = {a0*1000:.1f} mm")

print(f"\\n材料参数:")
print(f" 弹性模量 E = {E/1e9:.1f} GPa")
print(f" 断裂韧性 KIC = {KIC/1e6:.1f} MPa√m")
print(f" 内聚强度 σc = {sigma_c/1e6:.1f} MPa")
print(f" 断裂能 Gc = {Gc:.1f} J/m²")

# 加载参数
sigma_applied = np.linspace(0, 400e6, 100) # 远场应力

# 裂纹扩展模拟
a = a0
crack_length = [a0]
load_history = [0]
damage_zone = [0]
displacement = [0]

for sigma in sigma_applied[1:]:
# 计算应力强度因子
Y = 1.12 # 边缘裂纹修正
KI = Y * sigma * np.sqrt(np.pi * a)

# 计算J积分
J = KI**2 / E

# 裂纹尖端损伤区（简化模型）
# 损伤区尺寸与塑性区相关
rp = (KI / sigma_y)**2 / (2 * np.pi)
damage_size = 2 * rp

# 损伤区内平均损伤
D_avg = min(J / (S * 1e6) * 10, 0.5)

# 裂纹扩展判据
# 1. KI >= KIC（脆性扩展）
# 2. J >= Gc（韧性扩展）
# 3. D >= Dc（损伤诱发扩展）

if KI >= KIC or J >= Gc or D_avg >= Dc:
# 裂纹扩展
da = 0.5e-3 # 扩展步长
a = min(a + da, W/2)

crack_length.append(a)
load_history.append(sigma)
damage_zone.append(damage_size)

# 计算位移（柔度法）
C = (1 – nu**2) / E * (2 * np.pi * a**2) / B # 简化柔度
disp = C * sigma * W * B
displacement.append(disp)

crack_length = np.array(crack_length)
load_history = np.array(load_history)
damage_zone = np.array(damage_zone)
displacement = np.array(displacement)

print(f"\\n结果:")
print(f" 最终裂纹长度 = {crack_length[-1]*1000:.2f} mm")
print(f" 裂纹扩展量 = {(crack_length[-1]-a0)*1000:.2f} mm")
print(f" 最大损伤区尺寸 = {np.max(damage_zone)*1000:.3f} mm")

# 可视化
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 图1: 裂纹长度 vs 载荷
ax1 = axes[0, 0]
ax1.plot(load_history/1e6, (crack_length – a0)*1000, 'b-', linewidth=2)
ax1.set_xlabel('远场应力 σ (MPa)', fontsize=11)
ax1.set_ylabel('裂纹扩展量 Δa (mm)', fontsize=11)
ax1.set_title('裂纹扩展曲线', fontsize=12, fontweight='bold')
ax1.grid(True, alpha=0.3)

# 图2: 载荷-位移曲线
ax2 = axes[0, 1]
ax2.plot(displacement*1000, load_history/1e6, 'r-', linewidth=2)
ax2.set_xlabel('位移 δ (mm)', fontsize=11)
ax2.set_ylabel('应力 σ (MPa)', fontsize=11)
ax2.set_title('载荷-位移曲线', fontsize=12, fontweight='bold')
ax2.grid(True, alpha=0.3)

# 图3: 损伤区演化
ax3 = axes[1, 0]
ax3.plot(load_history/1e6, damage_zone*1000, 'g-', linewidth=2)
ax3.set_xlabel('远场应力 σ (MPa)', fontsize=11)
ax3.set_ylabel('损伤区尺寸 (mm)', fontsize=11)
ax3.set_title('裂纹尖端损伤区演化', fontsize=12, fontweight='bold')
ax3.grid(True, alpha=0.3)

# 图4: 裂纹尖端场可视化（概念图）
ax4 = axes[1, 1]

# 创建网格
x = np.linspace(-20e-3, 20e-3, 100)
y = np.linspace(-20e-3, 20e-3, 100)
X, Y = np.meshgrid(x, y)

# 极坐标
r = np.sqrt(X**2 + Y**2)
theta = np.arctan2(Y, X)

# 裂纹尖端应力场（简化）
KI_final = 1.12 * load_history[-1] * np.sqrt(np.pi * crack_length[-1])
sigma_field = KI_final / np.sqrt(2 * np.pi * (r + 1e-6)) * np.cos(theta/2) * \\
(1 + np.sin(theta/2)**2)

# 掩蔽裂纹面
sigma_field[X < 0] = np.nan

# 绘制
im = ax4.contourf(X*1000, Y*1000, sigma_field/1e6, levels=20, cmap='jet')
ax4.contour(X*1000, Y*1000, sigma_field/1e6, levels=5, colors='k', linewidths=0.5)

# 绘制裂纹
crack_x = np.linspace(-crack_length[-1]*1000, 0, 50)
crack_y = np.zeros_like(crack_x)
ax4.plot(crack_x, crack_y, 'k-', linewidth=3)
ax4.plot(crack_x, crack_y – 0.5, 'k-', linewidth=1) # 裂纹开口
ax4.plot(crack_x, crack_y + 0.5, 'k-', linewidth=1)

ax4.set_xlabel('x (mm)', fontsize=11)
ax4.set_ylabel('y (mm)', fontsize=11)
ax4.set_title('裂纹尖端应力场', fontsize=12, fontweight='bold')
ax4.set_aspect('equal')
plt.colorbar(im, ax=ax4, label='σ (MPa)')

plt.tight_layout()
plt.savefig('output/case6_crack_damage_interaction.png', dpi=150, bbox_inches='tight')
print("\\n案例6结果已保存")

return load_history, crack_length, displacement, damage_zone

# =============================================================================
# 主程序
# =============================================================================
if __name__ == "__main__":
# 运行所有案例
print("\\n开始运行所有案例…\\n")

# 案例1: Griffith脆性断裂分析
case1_griffith_fracture()

# 案例2: J积分计算
case2_j_integral()

# 案例3: Paris疲劳裂纹扩展
case3_paris_fatigue()

# 案例4: 蠕变损伤演化
case4_creep_damage()

# 案例5: 多轴应力状态下的损伤耦合
case5_multiaxial_damage()

# 案例6: 裂纹-损伤相互作用
case6_crack_damage_interaction()

print("\\n" + "="*70)
print("所有案例运行完成！")
print("结果图片保存在 output/ 目录")
print("="*70)