多场耦合优化-主题057-增材制造过程仿真

主题057：增材制造过程仿真

1. 增材制造技术概述

1.1 增材制造的定义与特点

增材制造（Additive Manufacturing, AM），又称3D打印，是一种通过逐层堆积材料来制造三维实体的技术。与传统的减材制造（切削、磨削）和等材制造（铸造、锻造）不同，增材制造基于"自下而上"的制造理念。

增材制造的核心优势：

• 设计自由度高：可制造复杂几何形状、内部空腔、晶格结构
• 材料利用率高：近净成形，减少材料浪费
• 生产周期短：无需模具，快速原型制造
• 个性化定制：小批量、多品种生产的经济性
• 功能集成：一体化制造，减少装配工序

主要技术分类：

1.2 金属增材制造的应用领域

航空航天：

• 发动机复杂构件（燃烧室、涡轮叶片）
• 结构件减重优化（拓扑优化结构）
• 快速维修与备件制造

医疗器械：

• 个性化植入物（髋臼杯、椎间融合器）
• 手术导板与模型
• 多孔骨整合结构

模具制造：

• 随形冷却水道模具
• 复杂型腔模具
• 模具修复与再制造

汽车工业：

• 轻量化结构件
• 定制零部件
• 小批量试制

1.3 增材制造的多物理场耦合特征

金属增材制造过程涉及复杂的多物理场耦合：

热场：

• 高能量密度热源（激光、电子束）
• 快速加热与冷却（10³-10⁶ K/s）
• 温度梯度引起的非均匀热膨胀

流场：

• 熔池内的Marangoni对流
• 表面张力驱动的流动
• 蒸发引起的反冲压力

应力场：

• 热应力与相变应力
• 残余应力积累
• 变形与开裂风险

相变：

• 熔化与凝固
• 固态相变
• 组织演化

质量传递：

• 元素蒸发
• 粉末输送
• 熔池与基板的物质交换

1.4 增材制造仿真的重要性

工艺开发：

• 工艺参数优化（功率、速度、扫描策略）
• 可制造性评估
• 新材料工艺窗口确定

质量控制：

• 缺陷预测（气孔、未熔合、裂纹）
• 微观结构预测
• 力学性能预测

变形控制：

• 残余应力分布预测
• 变形量预测
• 支撑结构优化

成本降低：

• 减少试错成本
• 缩短工艺开发周期
• 提高一次成功率

2. 金属增材制造工艺原理

2.1 选择性激光熔化（SLM）

工艺原理：

SLM使用高功率激光束（通常为光纤激光器，功率100-1000W）选择性熔化金属粉末层，逐层堆积形成致密金属零件。

工艺流程：

关键工艺参数：

能量密度：
$$E_v=\frac{P}{v\cdot h\cdot t}$$

其中：

• $P$：激光功率
• $v$：扫描速度
• $h$：扫描间距
• $t$：层厚

2.2 电子束熔化（EBM）

工艺特点：

EBM在真空环境中使用高能电子束（功率可达3kW以上）熔化金属粉末，主要优势包括：

• 高预热温度（可达1000°C以上），降低残余应力
• 高能量密度，可加工高熔点金属
• 真空环境，避免氧化

与SLM的比较：

2.3 激光近净成形（LENS/DED）

工艺原理：

LENS（Laser Engineered Net Shaping）或DED（Directed Energy Deposition）通过同轴或侧向送粉/送丝，激光熔化沉积材料，适用于：

• 大型零件制造
• 零件修复与再制造
• 梯度材料制造
• 多材料制造

送粉方式：

• 同轴送粉：粉末从喷嘴中心或环形送出
• 侧向送粉：粉末从侧面送入熔池

工艺优势：

• 高沉积率（可达1kg/h以上）
• 可修复现有零件
• 可制造梯度功能材料

2.4 粉末特性与铺粉质量

粉末特性要求：

铺粉质量评估：

• 粉末层厚度均匀性
• 表面平整度
• 孔隙率
• 粉末密实度

3. 熔池动力学与传热分析

3.1 热源模型

高斯热源模型：

激光/电子束的能量分布通常用高斯分布描述：
$$q(r)=\frac{3Q}{\pi r_b^2}\exp \left(-\frac{3r^2}{r_b^2}\right)$$

其中：

• $Q=\eta P$：有效热输入功率
• $\eta$：吸收率
• $r_b$：光束半径
• $r$：距光束中心的距离

体积热源模型：

考虑激光穿透效应的体积热源：
$$q(x,y,z)=\frac{3Q}{\pi r_b^{2}H}\exp \left(-\frac{3(x^2+y^2)}{r_b^2}\right)\exp \left(-\frac{3z}{H}\right)$$

其中 $H$ 为穿透深度。

双椭球热源模型：

更适合描述深熔小孔的Rosenthal双椭球模型：
$$q(x,y,z)=\frac{6\sqrt{3}Q}{\pi \sqrt{\pi }abc}\exp \left(-\frac{3x^2}{a^2}-\frac{3y^2}{b^2}-\frac{3z^2}{c^2}\right)$$

3.2 熔池传热分析

控制方程：

能量守恒方程：
$$\rho c_p\frac{\partial T}{\partial t}=\nabla \cdot (k\nabla T)+{\dot{q}}_{laser}-{\dot{q}}_{evap}-{\dot{q}}_{rad}$$

其中：

• ${\dot{q}}_{laser}$：激光热输入
• ${\dot{q}}_{evap}$：蒸发散热
• ${\dot{q}}_{rad}$：辐射散热

边界条件：

相变潜热处理：

采用等效热容法处理熔化潜热：
$$c_{p,eff}=c_p+L_f\frac{d{f}_L}{dT}$$

其中 $f_L$ 为液相分数。

3.3 熔池流动分析

熔池内的驱动力：

表面张力温度系数 $\partial \gamma /\partial T$ 通常为负值（大多数金属），导致熔池表面从中心向外流动。

Navier-Stokes方程：

熔池内的流动控制方程：
$$\rho \frac{D\mathbf{u}}{Dt}=-\nabla p+\mu {\nabla }^2\mathbf{u}+{\mathbf{F}}_b+{\mathbf{F}}_{ST}$$
$$\nabla \cdot \mathbf{u}=0$$

其中 ${\mathbf{F}}_{ST}$ 为表面张力项。

熔池几何特征：

熔池宽度：
$$W_{pool}\approx 2\sqrt{\frac{\alpha P}{\pi v\rho c_p(T_m-T_0)}}$$

熔池深度：
$$D_{pool}\approx \sqrt{\frac{\alpha P}{\pi vk}}$$

熔池长度：
$$L_{pool}\approx \frac{P}{2\pi k(T_m-T_0)}$$

3.4 蒸发与等离子体效应

蒸发速率：

根据Langmuir方程：
$${\dot{m}}_{evap}=\sqrt{\frac{M}{2\pi RT}}{P}_{sat}(T)$$

其中饱和蒸气压：
$$P_{sat}(T)=P_0\exp \left(-\frac{\Delta H_{evap}}{RT}\right)$$

反冲压力：

蒸发产生的反冲压力：
$$P_{recoil}=\frac{{\dot{m}}_{evap}^{2}RT}{M}$$

小孔效应（Keyholing）：

当能量密度足够高时，反冲压力克服表面张力形成深熔小孔：
$$P_{recoil}>\frac{2\gamma }{R_{keyhole}}$$

小孔模式的优势：

• 增加熔深
• 提高能量吸收率
• 改善侧壁熔合

小孔模式的风险：

• 小孔不稳定导致气孔
• 蒸发损失
• 飞溅

4. 凝固组织与微观结构预测

4.1 凝固理论

凝固条件：

增材制造的凝固特点：

• 极高的冷却速率：10³-10⁸ K/s
• 极高的温度梯度：10⁶-10⁸ K/m
• 定向凝固特征明显
• 非平衡凝固

G-R图（温度梯度-生长速率）：

凝固组织取决于温度梯度 $G$ 和凝固速率 $R$：

• 平面晶：高 $G/R$ 比值
• 胞状晶：中等 $G/R$ 比值
• 柱状枝晶：低 $G/R$ 比值
• 等轴晶：低 $G$ 或高过冷

枝晶间距：

一次枝晶臂间距：
$${\lambda }_1=A(G^{-0.5}{R}^{-0.25})$$

二次枝晶臂间距：
$${\lambda }_2=B(\dot{T}{)}^{-0.33}$$

其中 $\dot{T}=G\cdot R$ 为冷却速率。

4.2 晶粒生长模型

元胞自动机（CA）模型：

CA模型模拟晶粒形核与生长：

形核模型：

连续形核模型：
$$n(\Delta T)=\frac{n_{max}}{\sqrt{2\pi }\Delta T_{\sigma }}\exp \left(-\frac{(\Delta T-\Delta T_N{)}^2}{2\Delta T_{\sigma }^2}\right)$$

生长动力学：

枝晶尖端生长速率：
$$V=a_2\Delta T^2+a_3\Delta T^3$$

其中 $a_2,a_3$ 为材料常数。

4.3 相场模型

相场方程：

相场变量 $\varphi$ 描述固液界面（$\varphi =-1$ 固相，$\varphi =1$ 液相）：
$$\tau \frac{\partial \varphi }{\partial t}=W^2{\nabla }^2\varphi +\varphi -{\varphi }^3-\lambda (1-{\varphi }^2{)}^2(T-T_m)$$

热扩散方程：
$$\frac{\partial T}{\partial t}=\alpha {\nabla }^{2}T+\frac{L_f}{c_p}\frac{\partial f_S}{\partial t}$$

溶质扩散：

对于合金，考虑溶质再分配：
$$\frac{\partial C}{\partial t}=\nabla \cdot (D\nabla C)+\nabla \cdot [DC(1-k)\frac{\nabla \varphi }{1-\varphi +k(1+\varphi )}]$$

4.4 微观结构-性能关系

Hall-Petch关系：

晶粒尺寸与屈服强度的关系：
$${\sigma }_y={\sigma }_0+k_y{d}^{-0.5}$$

其中：

• ${\sigma }_0$：摩擦应力
• $k_y$：Hall-Petch系数
• $d$：平均晶粒尺寸

增材制造金属的典型组织特征：

5. 残余应力与变形分析

5.1 残余应力的形成机理

热应力：

温度梯度引起的不均匀热膨胀：
$${\sigma }_{thermal}=E\alpha \Delta T$$

其中：

• $E$：弹性模量
• $\alpha$：热膨胀系数
• $\Delta T$：温度变化

相变应力：

固态相变伴随的体积变化：
$${\sigma }_{phase}=E{\epsilon }_{phase}$$

塑性应变：

超过屈服强度时产生的永久变形：
$${\epsilon }_{plastic}={\epsilon }_{total}-{\epsilon }_{elastic}-{\epsilon }_{thermal}$$

残余应力的分布特征：

• 纵向应力（平行于扫描方向）：通常为拉应力
• 横向应力（垂直于扫描方向）：变化复杂
• 厚度方向应力：表面压应力、心部拉应力

5.2 热-弹-塑性分析

本构方程：

增量形式的热-弹-塑性本构关系：
$$d{\sigma }_{ij}=D_{ijkl}^{ep}(d{\epsilon }_{kl}-d{\epsilon }_{kl}^{th}-d{\epsilon }_{kl}^{phase})$$

其中：

• $D_{ijkl}^{ep}$：弹塑性刚度张量
• $d{\epsilon }_{kl}^{th}=\alpha dT{\delta }_{kl}$：热应变增量
• $d{\epsilon }_{kl}^{phase}$：相变应变增量

屈服准则：

Von Mises屈服准则：
$$f=\sqrt{\frac{3}{2}{s}_{ij}{s}_{ij}}-{\sigma }_y(T,\dot{\epsilon })=0$$

其中 $s_{ij}$ 为偏应力张量。

流动法则：

关联流动法则：
$$d{\epsilon }_{ij}^p=d\lambda \frac{\partial f}{\partial {\sigma }_{ij}}$$

5.3 固有应变法

基本原理：

固有应变法基于弹性分析，通过引入等效载荷来考虑塑性变形和相变的影响。

固有应变分量：
$${\epsilon }^{\ast }={\epsilon }^p+{\epsilon }^{th}+{\epsilon }^{phase}$$

等效载荷：
$${\mathbf{F}}^{\ast }={\int }_V{\mathbf{B}}^T\mathbf{D}{\epsilon }^{\ast }dV$$

固有应变法的优势：

• 计算效率高
• 适合大型构件
• 便于参数化研究

5.4 变形预测与控制

变形类型：

变形预测模型：

基于弹性变形的简化模型：
$$\delta =\frac{\alpha \Delta T{L}^2}{h}$$

其中：

• $L$：特征长度
• $h$：构件厚度

变形控制策略：

6. 扫描策略与工艺优化

6.1 扫描策略类型

单向扫描（Unidirectional）：

• 所有层沿同一方向扫描
• 优点：简单、效率高
• 缺点：各向异性明显

双向扫描（Bidirectional/Zigzag）：

• 相邻扫描线反向
• 优点：减少扫描时间
• 缺点： turnaround点质量差

条纹扫描（Stripe）：

• 将区域分成条带分别扫描
• 优点：减少残余应力
• 缺点：条带间搭接问题

棋盘扫描（Chessboard/Island）：

• 将区域分成小块交替扫描
• 优点：均匀化热分布
• 缺点：复杂、效率低

轮廓扫描（Contour）：

• 先扫描轮廓再填充内部
• 优点：表面质量好
• 缺点：轮廓-填充搭接问题

6.2 扫描参数优化

能量密度优化：

最优能量密度范围：
$$E_v^{opt}=50-150{\text{ J/mm}}^3$$

能量密度过低：未熔合、气孔
能量密度过高：球化、蒸发损失

扫描间距优化：

搭接率计算：
$${\eta }_{overlap}=\frac{2r_b-h}{2r_b}\times 100\%$$

推荐搭接率：30-50%

层厚优化：

层厚与粉末粒径的关系：
$$t_{layer}=(0.5-0.8)d_{powder}$$

6.3 支撑结构设计

支撑的功能：

• 防止翘曲变形
• 传导热量
• 支撑悬空结构
• 防止坍塌

支撑类型：

支撑设计原则：

• 支撑角度：通常>45°需要支撑
• 支撑密度：30-70%
• 支撑厚度：0.5-2mm

6.4 工艺优化方法

实验设计（DOE）：

通过正交实验或响应面法优化工艺参数：

• 因素：功率、速度、扫描间距、层厚
• 响应：致密度、表面粗糙度、硬度

数值优化：

基于仿真的优化：

• 目标函数：最小化残余应力或变形
• 约束条件：致密度>99.5%
• 优化算法：遗传算法、粒子群算法

机器学习优化：

利用历史数据训练模型：

• 输入：工艺参数、材料特性
• 输出：质量指标
• 模型：神经网络、随机森林

7. 数值方法与求解策略

7.1 有限元法（FEM）

热分析有限元方程：
$$\mathbf{C}\dot{\mathbf{T}}+\mathbf{KT}=\mathbf{Q}$$

其中：

• $\mathbf{C}$：热容矩阵
• $\mathbf{K}$：热传导矩阵
• $\mathbf{Q}$：热载荷向量

时间积分：

向后Euler法（无条件稳定）：
$$\mathbf{C}\frac{{\mathbf{T}}^{n+1}-{\mathbf{T}}^n}{\Delta t}+\mathbf{K}{\mathbf{T}}^{n+1}={\mathbf{Q}}^{n+1}$$

生死单元技术：

逐步激活单元模拟增材过程：

• 未激活单元：极低刚度
• 激活单元：正常材料属性
• 激活准则：扫描前沿到达

7.2 计算流体动力学（CFD）

熔池流动求解：

采用VOF（Volume of Fluid）方法追踪自由表面：
$$\frac{\partial F}{\partial t}+\mathbf{u}\cdot \nabla F=0$$

其中 $F$ 为流体体积分数。

表面张力处理：

CSF（Continuum Surface Force）模型：
$${\mathbf{F}}_{ST}=\gamma \kappa \nabla F$$

其中 $\kappa$ 为界面曲率。

7.3 多尺度仿真

宏观-微观耦合：

尺度桥接方法：

• 顺序耦合：先宏观后微观
• 并发耦合：同时求解
• 代表性体积元（RVE）：提取等效性能

7.4 高性能计算

并行计算策略：

• 区域分解：将计算域分配给不同处理器
• 负载均衡：平衡各处理器工作量
• 通信优化：减少处理器间数据交换

GPU加速：

利用GPU的并行计算能力：

• 矩阵运算
• 迭代求解
• 图像处理

8. 工程案例分析

8.1 案例一：316L不锈钢SLM工艺优化

问题描述：

优化316L不锈钢的SLM工艺参数，获得高致密度（>99.5%）和良好的表面质量。

材料参数：

• 密度：7980 kg/m³
• 熔点：1400°C
• 热导率：15 W/(m·K)
• 比热：500 J/(kg·K)

工艺参数范围：

• 激光功率：200-400 W
• 扫描速度：800-1500 mm/s
• 扫描间距：80-120 μm
• 层厚：30-50 μm

仿真分析：

建立三维瞬态热传导模型，分析：

优化结果：

最优参数：

• 功率：280 W
• 速度：1000 mm/s
• 扫描间距：100 μm
• 层厚：40 μm
• 能量密度：70 J/mm³

预测致密度：99.7%

8.2 案例二：Ti-6Al-4V EBM残余应力分析

问题描述：

分析Ti-6Al-4V在EBM过程中的残余应力分布，评估开裂风险。

材料特性：

• 熔点：1668°C
• 热膨胀系数：9×10⁻⁶ K⁻¹
• 弹性模量：110 GPa
• 屈服强度：900 MPa

工艺条件：

• 预热温度：650°C
• 电子束功率：3 kW
• 扫描速度：500 mm/s

仿真模型：

热-弹-塑性分析，考虑：

• 温度依赖的材料性能
• 马氏体相变
• 预热温度的影响

结果分析：

• 最大残余应力：600 MPa（表面拉应力）
• 变形量：0.3 mm（100mm构件）
• 开裂风险评估：低风险

优化建议：

• 提高预热温度至750°C
• 采用棋盘扫描策略
• 增加支撑结构密度

8.3 案例三：叶片修复的LENS过程仿真

问题描述：

仿真航空发动机叶片的LENS修复过程，优化修复参数。

修复区域：

• 磨损区域尺寸：20×10×3 mm
• 基体材料：Inconel 718
• 修复材料：Inconel 718粉末

工艺参数：

• 激光功率：2 kW
• 送粉速率：5 g/min
• 扫描速度：300 mm/min
• 层厚：0.5 mm

多物理场分析：

结果：

• 熔池深度：1.2 mm
• 稀释率：15%
• 修复区硬度：与基体匹配
• 残余应力：可接受范围

8.4 案例四：晶格结构的SLM制造仿真

问题描述：

优化晶格结构的SLM制造工艺，预测变形并设计支撑。

晶格结构：

• 类型：体心立方（BCC）
• 单胞尺寸：5 mm
• 杆径：0.8 mm
• 相对密度：15%

挑战：

• 细长杆件易变形
• 热量积聚
• 支撑去除困难

仿真内容：

优化结果：

• 采用点支撑策略
• 支撑密度：50%
• 预测最大变形：0.15 mm
• 实际变形：0.12 mm（验证）

"""
主题057：增材制造过程仿真
Additive Manufacturing Process Simulation

本程序演示增材制造过程的多物理场仿真方法，包含以下4个案例：
1. 熔池温度场仿真（移动高斯热源）
2. 扫描策略对比（不同扫描路径的热积累）
3. 残余应力分析（热-弹-塑性模型）
4. 凝固组织预测（晶粒生长模型）
"""

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation, PillowWriter
from matplotlib.patches import Rectangle
from scipy.integrate import odeint
import os

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

# 创建输出目录
OUTPUT_DIR = 'd:\\\\文档\\\\500仿真领域\\\\工程仿真\\\\多场耦合优化\\\\主题057\\\\output\\\\'
os.makedirs(OUTPUT_DIR, exist_ok=True)

print("="*60)
print("主题057：增材制造过程仿真")
print("Additive Manufacturing Process Simulation")
print("="*60)
print("\\n本程序演示增材制造过程的仿真方法，包含以下4个案例：")
print("1. 熔池温度场仿真（移动高斯热源）")
print("2. 扫描策略对比（不同扫描路径的热积累）")
print("3. 残余应力分析（热-弹-塑性模型）")
print("4. 凝固组织预测（晶粒生长模型）")
print("="*60)

def case1_melt_pool_temperature():
"""
案例1: 熔池温度场仿真
模拟移动高斯热源作用下的温度场分布
"""
print("\\n" + "="*60)
print("案例1: 熔池温度场仿真")
print("="*60)

# 316L不锈钢材料参数
material = {
'name': '316L Stainless Steel',
'rho': 7980, # 密度 (kg/m³)
'cp': 500, # 比热 (J/(kg·K))
'k': 15, # 热导率 (W/(m·K))
'Tm': 1400, # 熔点 (°C)
'T0': 25, # 初始温度 (°C)
'L_f': 2.7e5, # 熔化潜热 (J/kg)
}

# 激光工艺参数
laser = {
'P': 300, # 功率 (W)
'eta': 0.35, # 吸收率
'r_b': 0.1e-3, # 光束半径 (m)
'v': 1.0, # 扫描速度 (m/s)
}

# 计算热扩散系数
alpha = material['k'] / (material['rho'] * material['cp'])

print(f"材料: {material['name']}")
print(f"激光功率: {laser['P']} W")
print(f"扫描速度: {laser['v']*1000:.0f} mm/s")
print(f"光束半径: {laser['r_b']*1e6:.0f} μm")
print(f"热扩散系数: {alpha:.2e} m²/s")

# 创建计算域
Lx, Ly = 4e-3, 2e-3 # 4mm x 2mm
Nx, Ny = 200, 100
x = np.linspace(-Lx/2, Lx/2, Nx)
y = np.linspace(0, Ly, Ny)
X, Y = np.meshgrid(x, y)

# 高斯热源模型（准稳态）
def gaussian_heat_source(X, Y, laser_params, material_params):
"""计算移动高斯热源作用下的温度场（Rosenthal解的简化形式）"""
P_eff = laser_params['P'] * laser_params['eta']
r_b = laser_params['r_b']
v = laser_params['v']
k = material_params['k']
T0 = material_params['T0']

# 无量纲坐标（随热源移动）
xi = v * X / (2 * alpha)
eta_coord = v * Y / (2 * alpha)

# 温度场（简化解析解）
r = np.sqrt(X**2 + Y**2)
theta = P_eff / (2 * np.pi * k * r_b) * np.exp(-v * (X + r) / (2 * alpha))
T = T0 + theta * np.exp(-(X**2 + Y**2) / (2 * r_b**2))

return T

# 计算温度场
T_field = gaussian_heat_source(X, Y, laser, material)

# 计算熔池尺寸
melt_pool_mask = T_field > material['Tm']
if np.any(melt_pool_mask):
melt_indices = np.where(melt_pool_mask)
melt_width = (np.max(X[melt_indices]) – np.min(X[melt_indices])) * 1e6 # μm
melt_depth = (np.max(Y[melt_indices]) – np.min(Y[melt_indices])) * 1e6 # μm
else:
melt_width, melt_depth = 0, 0

# 计算冷却速率
cooling_rate = laser['v'] * (material['Tm'] – material['T0']) / (2 * laser['r_b'])

print(f"\\n熔池尺寸:")
print(f" 熔池宽度: {melt_width:.0f} μm")
print(f" 熔池深度: {melt_depth:.0f} μm")
print(f" 冷却速率: {cooling_rate/1e3:.1f} K/ms")

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

# 1. 温度场分布
ax1 = axes[0, 0]
levels = np.linspace(25, 1500, 50)
cf1 = ax1.contourf(X*1e3, Y*1e3, T_field, levels=levels, cmap='hot')
ax1.contour(X*1e3, Y*1e3, T_field, levels=[material['Tm']], colors='blue', linewidths=2)
ax1.set_xlabel('X (mm)', fontsize=11)
ax1.set_ylabel('Y (mm)', fontsize=11)
ax1.set_title('Temperature Field (Moving Gaussian Heat Source)', fontsize=12)
ax1.set_aspect('equal')
cbar1 = plt.colorbar(cf1, ax=ax1)
cbar1.set_label('Temperature (°C)', fontsize=10)

# 2. 中心线温度分布
ax2 = axes[0, 1]
center_idx = Ny // 2
ax2.plot(x*1e3, T_field[center_idx, :], 'r-', linewidth=2)
ax2.axhline(y=material['Tm'], color='b', linestyle='–', label=f"Melting Point: {material['Tm']}°C")
ax2.set_xlabel('X (mm)', fontsize=11)
ax2.set_ylabel('Temperature (°C)', fontsize=11)
ax2.set_title('Temperature Profile along Centerline', fontsize=12)
ax2.legend()
ax2.grid(True, alpha=0.3)

# 3. 深度方向温度分布
ax3 = axes[1, 0]
center_x_idx = Nx // 2
ax3.plot(T_field[:, center_x_idx], y*1e3, 'b-', linewidth=2)
ax3.axvline(x=material['Tm'], color='r', linestyle='–', label=f"Melting Point: {material['Tm']}°C")
ax3.set_xlabel('Temperature (°C)', fontsize=11)
ax3.set_ylabel('Depth (mm)', fontsize=11)
ax3.set_title('Temperature Profile along Depth', fontsize=12)
ax3.legend()
ax3.grid(True, alpha=0.3)
ax3.invert_yaxis()

# 4. 工艺参数影响
ax4 = axes[1, 1]
powers = np.linspace(100, 500, 5)
melt_widths = []
melt_depths = []

for P in powers:
laser_temp = laser.copy()
laser_temp['P'] = P
T_temp = gaussian_heat_source(X, Y, laser_temp, material)
mask = T_temp > material['Tm']
if np.any(mask):
indices = np.where(mask)
mw = (np.max(X[indices]) – np.min(X[indices])) * 1e6
md = (np.max(Y[indices]) – np.min(Y[indices])) * 1e6
else:
mw, md = 0, 0
melt_widths.append(mw)
melt_depths.append(md)

ax4.plot(powers, melt_widths, 'ro-', linewidth=2, markersize=8, label='Melt Pool Width')
ax4.plot(powers, melt_depths, 'bs-', linewidth=2, markersize=8, label='Melt Pool Depth')
ax4.set_xlabel('Laser Power (W)', fontsize=11)
ax4.set_ylabel('Melt Pool Size (μm)', fontsize=11)
ax4.set_title('Effect of Laser Power on Melt Pool Size', fontsize=12)
ax4.legend()
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(OUTPUT_DIR + 'case1_melt_pool_temperature.png', dpi=150, bbox_inches='tight')
plt.close()

print(f"案例1完成！图片已保存至 {OUTPUT_DIR}case1_melt_pool_temperature.png")

return T_field, melt_width, melt_depth

def case2_scanning_strategies():
"""
案例2: 扫描策略对比
比较不同扫描路径的热积累效应
"""
print("\\n" + "="*60)
print("案例2: 扫描策略对比")
print("="*60)

# 基板尺寸
Lx, Ly = 10e-3, 10e-3 # 10mm x 10mm
Nx, Ny = 100, 100
x = np.linspace(0, Lx, Nx)
y = np.linspace(0, Ly, Ny)
dx, dy = x[1] – x[0], y[1] – y[0]

# 材料参数（316L）
material = {
'rho': 7980,
'cp': 500,
'k': 15,
'alpha': 15 / (7980 * 500), # 热扩散系数
}

# 激光参数
laser = {
'P': 200, # W
'eta': 0.35,
'r_b': 0.1e-3, # m
'v': 0.5, # m/s
'dt': 1e-4, # s
}

# 初始化温度场
T0 = 25 # 初始温度

def simulate_scanning(scan_path, T_initial, n_layers=1):
"""模拟扫描过程的温度场演化"""
T = T_initial.copy()
T_history = [T.copy()]

for layer in range(n_layers):
for point in scan_path:
# 激光位置
x_laser, y_laser = point

# 计算距离
dist_sq = (X – x_laser)**2 + (Y – y_laser)**2

# 高斯热源
q = laser['P'] * laser['eta'] / (np.pi * laser['r_b']**2) * \\
np.exp(-dist_sq / laser['r_b']**2)

# 热传导（简化显式格式）
dT = material['alpha'] * laser['dt'] * (
(np.roll(T, 1, axis=0) – 2*T + np.roll(T, -1, axis=0)) / dy**2 +
(np.roll(T, 1, axis=1) – 2*T + np.roll(T, -1, axis=1)) / dx**2
) + q / (material['rho'] * material['cp']) * laser['dt']

T = T + dT

# 记录温度历史
if len(T_history) % 10 == 0:
T_history.append(T.copy())

return T, T_history

# 定义扫描路径
# 1. 单向扫描
path_unidirectional = []
for i in range(10):
y_pos = (i + 0.5) * Ly / 10
if i % 2 == 0:
for j in range(20):
x_pos = (j + 0.5) * Lx / 20
path_unidirectional.append((x_pos, y_pos))
else:
for j in range(19, -1, -1):
x_pos = (j + 0.5) * Lx / 20
path_unidirectional.append((x_pos, y_pos))

# 2. 棋盘扫描
path_chessboard = []
block_size = 2
for i in range(0, 10, block_size):
for j in range(0, 10, block_size):
if (i // block_size + j // block_size) % 2 == 0:
y_start = (i + 0.5) * Ly / 10
for ii in range(block_size):
y_pos = y_start + ii * Ly / 10
for jj in range(block_size):
x_pos = (j + jj + 0.5) * Lx / 10
path_chessboard.append((x_pos, y_pos))

# 创建网格
X, Y = np.meshgrid(x, y)

# 简化仿真：只计算最终温度分布
print("模拟单向扫描…")
T_uni = T0 * np.ones((Ny, Nx))
for point in path_unidirectional[:50]: # 简化：只扫描部分路径
x_laser, y_laser = point
dist_sq = (X – x_laser)**2 + (Y – y_laser)**2
q = laser['P'] * laser['eta'] / (np.pi * laser['r_b']**2) * \\
np.exp(-dist_sq / laser['r_b']**2)
dT = q / (material['rho'] * material['cp']) * 0.01
T_uni += dT

print("模拟棋盘扫描…")
T_chess = T0 * np.ones((Ny, Nx))
for point in path_chessboard[:50]:
x_laser, y_laser = point
dist_sq = (X – x_laser)**2 + (Y – y_laser)**2
q = laser['P'] * laser['eta'] / (np.pi * laser['r_b']**2) * \\
np.exp(-dist_sq / laser['r_b']**2)
dT = q / (material['rho'] * material['cp']) * 0.01
T_chess += dT

# 计算热积累指标
max_temp_uni = np.max(T_uni)
max_temp_chess = np.max(T_chess)
avg_temp_uni = np.mean(T_uni)
avg_temp_chess = np.mean(T_chess)

print(f"\\n热积累对比:")
print(f" 单向扫描 – 最高温度: {max_temp_uni:.1f}°C, 平均温度: {avg_temp_uni:.1f}°C")
print(f" 棋盘扫描 – 最高温度: {max_temp_chess:.1f}°C, 平均温度: {avg_temp_chess:.1f}°C")

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

# 1. 单向扫描温度场
ax1 = axes[0, 0]
cf1 = ax1.contourf(X*1e3, Y*1e3, T_uni, levels=20, cmap='hot')
ax1.set_xlabel('X (mm)', fontsize=11)
ax1.set_ylabel('Y (mm)', fontsize=11)
ax1.set_title(f'Unidirectional Scanning\\nMax Temp: {max_temp_uni:.1f}°C', fontsize=12)
ax1.set_aspect('equal')
plt.colorbar(cf1, ax=ax1)

# 2. 棋盘扫描温度场
ax2 = axes[0, 1]
cf2 = ax2.contourf(X*1e3, Y*1e3, T_chess, levels=20, cmap='hot')
ax2.set_xlabel('X (mm)', fontsize=11)
ax2.set_ylabel('Y (mm)', fontsize=11)
ax2.set_title(f'Chessboard Scanning\\nMax Temp: {max_temp_chess:.1f}°C', fontsize=12)
ax2.set_aspect('equal')
plt.colorbar(cf2, ax=ax2)

# 3. 扫描路径可视化
ax3 = axes[1, 0]
path_x = [p[0]*1e3 for p in path_unidirectional[:100]]
path_y = [p[1]*1e3 for p in path_unidirectional[:100]]
ax3.plot(path_x, path_y, 'b-', linewidth=1, alpha=0.7)
ax3.scatter(path_x[::10], path_y[::10], c='red', s=20)
ax3.set_xlabel('X (mm)', fontsize=11)
ax3.set_ylabel('Y (mm)', fontsize=11)
ax3.set_title('Unidirectional Scan Path', fontsize=12)
ax3.set_aspect('equal')
ax3.grid(True, alpha=0.3)

# 4. 温度对比
ax4 = axes[1, 1]
strategies = ['Unidirectional', 'Chessboard']
max_temps = [max_temp_uni, max_temp_chess]
avg_temps = [avg_temp_uni, avg_temp_chess]

x_pos = np.arange(len(strategies))
width = 0.35

bars1 = ax4.bar(x_pos – width/2, max_temps, width, label='Max Temperature', color='red', alpha=0.7)
bars2 = ax4.bar(x_pos + width/2, avg_temps, width, label='Avg Temperature', color='blue', alpha=0.7)

ax4.set_ylabel('Temperature (°C)', fontsize=11)
ax4.set_title('Temperature Comparison between Scanning Strategies', fontsize=12)
ax4.set_xticks(x_pos)
ax4.set_xticklabels(strategies)
ax4.legend()
ax4.grid(True, alpha=0.3, axis='y')

# 添加数值标签
for bars in [bars1, bars2]:
for bar in bars:
height = bar.get_height()
ax4.text(bar.get_x() + bar.get_width()/2., height,
f'{height:.1f}', ha='center', va='bottom', fontsize=9)

plt.tight_layout()
plt.savefig(OUTPUT_DIR + 'case2_scanning_strategies.png', dpi=150, bbox_inches='tight')
plt.close()

print(f"案例2完成！图片已保存至 {OUTPUT_DIR}case2_scanning_strategies.png")

return T_uni, T_chess

def case3_residual_stress():
"""
案例3: 残余应力分析
基于热-弹-塑性模型预测残余应力分布
"""
print("\\n" + "="*60)
print("案例3: 残余应力分析")
print("="*60)

# Ti-6Al-4V材料参数
material = {
'name': 'Ti-6Al-4V',
'E': 110e9, # 弹性模量 (Pa)
'nu': 0.31, # 泊松比
'alpha': 9e-6, # 热膨胀系数 (1/K)
'sigma_y': 900e6, # 屈服强度 (Pa)
'Tm': 1668, # 熔点 (°C)
'T0': 25, # 室温 (°C)
}

# 构件尺寸
L = 50e-3 # 50mm
H = 5e-3 # 5mm

# 离散
Nx, Ny = 100, 20
x = np.linspace(0, L, Nx)
y = np.linspace(0, H, Ny)
X, Y = np.meshgrid(x, y)

print(f"材料: {material['name']}")
print(f"构件尺寸: {L*1e3:.0f} mm x {H*1e3:.0f} mm")
print(f"弹性模量: {material['E']/1e9:.0f} GPa")
print(f"屈服强度: {material['sigma_y']/1e6:.0f} MPa")

# 简化温度场（模拟多层沉积后的温度分布）
# 假设表面温度高，内部温度低
T_max = 800 # 表面最高温度
T_min = 200 # 内部最低温度

T_field = T_min + (T_max – T_min) * np.exp(-Y / (H/3))

# 热应变
epsilon_th = material['alpha'] * (T_field – material['T0'])

# 简化应力计算（假设完全约束）
# 弹性应力 = E * epsilon_th
sigma_elastic = material['E'] * epsilon_th

# 考虑塑性（简化处理）
sigma_residual = np.minimum(sigma_elastic, material['sigma_y'])
sigma_residual = np.maximum(sigma_residual, -material['sigma_y'])

# 计算变形（简化模型）
# 假设构件两端自由，中心受约束
delta_T_avg = np.mean(T_field) – material['T0']
deformation = material['alpha'] * delta_T_avg * L

print(f"\\n残余应力分析结果:")
print(f" 最大残余应力: {np.max(sigma_residual)/1e6:.1f} MPa")
print(f" 最小残余应力: {np.min(sigma_residual)/1e6:.1f} MPa")
print(f" 平均温度: {np.mean(T_field):.1f}°C")
print(f" 预测变形: {deformation*1e6:.1f} μm")

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

# 1. 温度场
ax1 = axes[0, 0]
cf1 = ax1.contourf(X*1e3, Y*1e3, T_field, levels=20, cmap='hot')
ax1.set_xlabel('X (mm)', fontsize=11)
ax1.set_ylabel('Y (mm)', fontsize=11)
ax1.set_title('Temperature Distribution after Deposition', fontsize=12)
plt.colorbar(cf1, ax=ax1, label='Temperature (°C)')

# 2. 残余应力分布
ax2 = axes[0, 1]
cf2 = ax2.contourf(X*1e3, Y*1e3, sigma_residual/1e6, levels=20, cmap='RdBu_r')
ax2.set_xlabel('X (mm)', fontsize=11)
ax2.set_ylabel('Y (mm)', fontsize=11)
ax2.set_title('Residual Stress Distribution', fontsize=12)
plt.colorbar(cf2, ax=ax2, label='Stress (MPa)')

# 3. 沿厚度方向的应力分布
ax3 = axes[1, 0]
center_x = Nx // 2
ax3.plot(sigma_residual[:, center_x]/1e6, y*1e3, 'b-', linewidth=2)
ax3.axvline(x=0, color='k', linestyle='–', alpha=0.5)
ax3.axvline(x=material['sigma_y']/1e6, color='r', linestyle='–', alpha=0.5, label='Yield Strength')
ax3.axvline(x=-material['sigma_y']/1e6, color='r', linestyle='–', alpha=0.5)
ax3.set_xlabel('Residual Stress (MPa)', fontsize=11)
ax3.set_ylabel('Depth (mm)', fontsize=11)
ax3.set_title('Residual Stress along Thickness', fontsize=12)
ax3.legend()
ax3.grid(True, alpha=0.3)
ax3.invert_yaxis()

# 4. 不同预热温度的影响
ax4 = axes[1, 1]
preheat_temps = [25, 200, 400, 600, 800]
max_stresses = []
deformations = []

for T_pre in preheat_temps:
# 简化计算：预热降低温度梯度
T_eff = T_field – T_pre * 0.5
delta_T = np.mean(T_eff) – material['T0']
sigma = material['E'] * material['alpha'] * delta_T
sigma = min(sigma, material['sigma_y'])
deform = material['alpha'] * delta_T * L
max_stresses.append(sigma / 1e6)
deformations.append(deform * 1e6)

ax4_twin = ax4.twinx()
line1 = ax4.plot(preheat_temps, max_stresses, 'ro-', linewidth=2, markersize=8, label='Max Stress')
line2 = ax4_twin.plot(preheat_temps, deformations, 'bs-', linewidth=2, markersize=8, label='Deformation')

ax4.set_xlabel('Preheat Temperature (°C)', fontsize=11)
ax4.set_ylabel('Max Residual Stress (MPa)', fontsize=11, color='red')
ax4_twin.set_ylabel('Deformation (μm)', fontsize=11, color='blue')
ax4.set_title('Effect of Preheat Temperature', fontsize=12)
ax4.grid(True, alpha=0.3)

# 合并图例
lines = line1 + line2
labels = [l.get_label() for l in lines]
ax4.legend(lines, labels, loc='best')

plt.tight_layout()
plt.savefig(OUTPUT_DIR + 'case3_residual_stress.png', dpi=150, bbox_inches='tight')
plt.close()

print(f"案例3完成！图片已保存至 {OUTPUT_DIR}case3_residual_stress.png")

return sigma_residual, deformation

def case4_solidification_microstructure():
"""
案例4: 凝固组织预测
基于元胞自动机模型预测晶粒生长
"""
print("\\n" + "="*60)
print("案例4: 凝固组织预测")
print("="*60)

# 模拟区域
N = 200
L = 1e-3 # 1mm
dx = L / N

# 材料参数（316L）
material = {
'Tm': 1400, # 熔点 (°C)
'T0': 25, # 初始温度 (°C)
'G': 5e6, # 温度梯度 (K/m)
'R': 1e-3, # 凝固速率 (m/s)
}

print(f"模拟区域: {L*1e3:.1f} mm x {L*1e3:.1f} mm")
print(f"温度梯度: {material['G']/1e6:.1f} K/mm")
print(f"凝固速率: {material['R']*1e3:.1f} mm/s")

# 温度场（线性梯度）
x = np.linspace(0, L, N)
y = np.linspace(0, L, N)
X, Y = np.meshgrid(x, y)

# 固液界面位置（随时间移动）
interface_pos = L / 2

# 温度场
T_field = material['Tm'] + material['G'] * (Y – interface_pos)

# 元胞自动机参数
n_grains = 50 # 晶粒数量
np.random.seed(42)

# 初始化晶粒（随机形核位置）
nuclei_x = np.random.randint(0, N, n_grains)
nuclei_y = np.random.randint(N//2, N, n_grains) # 只在液相区形核
nuclei_ori = np.random.rand(n_grains) * 2 * np.pi # 随机取向

# 晶粒ID矩阵
grain_ids = np.zeros((N, N), dtype=int)
grain_angles = np.zeros(n_grains)

# 分配晶粒ID
for i in range(n_grains):
grain_ids[nuclei_y[i], nuclei_x[i]] = i + 1
grain_angles[i] = nuclei_ori[i]

# 简化晶粒生长模拟
# 使用距离变换模拟晶粒生长
from scipy.ndimage import distance_transform_edt

# 创建形核点掩码
nuclei_mask = grain_ids > 0

# 距离变换
distances, indices = distance_transform_edt(~nuclei_mask, return_indices=True)

# 根据距离分配晶粒（简化）
grain_structure = np.zeros((N, N))
for i in range(N):
for j in range(N):
if nuclei_mask[i, j]:
grain_structure[i, j] = grain_ids[i, j]
else:
# 找到最近的形核点
ni, nj = indices[0, i, j], indices[1, j, j]
if 0 <= ni < N and 0 <= nj < N:
grain_structure[i, j] = grain_ids[ni, nj]

# 计算晶粒尺寸统计
unique_grains = np.unique(grain_structure[grain_structure > 0])
grain_sizes = []
for gid in unique_grains:
size = np.sum(grain_structure == gid) * dx * 1e6 # μm²
grain_sizes.append(np.sqrt(size)) # 等效直径

avg_grain_size = np.mean(grain_sizes)
std_grain_size = np.std(grain_sizes)

# 计算冷却速率
cooling_rate = material['G'] * material['R'] # K/s

# 预测二次枝晶臂间距
SDAS = 50 * (cooling_rate)**(-0.33) # μm

print(f"\\n凝固组织预测结果:")
print(f" 晶粒数量: {len(unique_grains)}")
print(f" 平均晶粒尺寸: {avg_grain_size:.1f} μm")
print(f" 晶粒尺寸标准差: {std_grain_size:.1f} μm")
print(f" 冷却速率: {cooling_rate:.2e} K/s")
print(f" 预测SDAS: {SDAS:.2f} μm")

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

# 1. 温度场
ax1 = axes[0, 0]
cf1 = ax1.contourf(X*1e6, Y*1e6, T_field, levels=20, cmap='hot')
ax1.contour(X*1e6, Y*1e6, T_field, levels=[material['Tm']], colors='blue', linewidths=2)
ax1.set_xlabel('X (μm)', fontsize=11)
ax1.set_ylabel('Y (μm)', fontsize=11)
ax1.set_title('Temperature Field with Solidification Front', fontsize=12)
ax1.set_aspect('equal')
plt.colorbar(cf1, ax=ax1, label='Temperature (°C)')

# 2. 晶粒结构
ax2 = axes[0, 1]
# 为不同晶粒分配不同颜色
grain_structure_plot = grain_structure % 20 # 循环使用20种颜色
cf2 = ax2.imshow(grain_structure_plot, extent=[0, L*1e6, 0, L*1e6],
cmap='tab20', origin='lower')
ax2.set_xlabel('X (μm)', fontsize=11)
ax2.set_ylabel('Y (μm)', fontsize=11)
ax2.set_title('Predicted Grain Structure', fontsize=12)
ax2.set_aspect('equal')

# 3. 晶粒尺寸分布
ax3 = axes[1, 0]
ax3.hist(grain_sizes, bins=20, color='steelblue', alpha=0.7, edgecolor='black')
ax3.axvline(x=avg_grain_size, color='red', linestyle='–', linewidth=2,
label=f'Mean: {avg_grain_size:.1f} μm')
ax3.set_xlabel('Grain Size (μm)', fontsize=11)
ax3.set_ylabel('Frequency', fontsize=11)
ax3.set_title('Grain Size Distribution', fontsize=12)
ax3.legend()
ax3.grid(True, alpha=0.3, axis='y')

# 4. 工艺参数对组织的影响
ax4 = axes[1, 1]
G_values = np.linspace(1e6, 10e6, 10) # K/m
R_values = np.linspace(0.5e-3, 5e-3, 10) # m/s

SDAS_values = []
for G in G_values:
for R in R_values:
CR = G * R
sdas = 50 * (CR)**(-0.33)
SDAS_values.append(sdas)

# 绘制SDAS随冷却速率的变化
cooling_rates = G_values * 1e-3 # 简化
sdas_plot = [50 * (cr)**(-0.33) for cr in cooling_rates]

ax4.loglog(cooling_rates, sdas_plot, 'bo-', linewidth=2, markersize=6)
ax4.set_xlabel('Cooling Rate (K/s)', fontsize=11)
ax4.set_ylabel('SDAS (μm)', fontsize=11)
ax4.set_title('Effect of Cooling Rate on SDAS', fontsize=12)
ax4.grid(True, alpha=0.3, which='both')

plt.tight_layout()
plt.savefig(OUTPUT_DIR + 'case4_solidification_microstructure.png', dpi=150, bbox_inches='tight')
plt.close()

print(f"案例4完成！图片已保存至 {OUTPUT_DIR}case4_solidification_microstructure.png")

return grain_structure, avg_grain_size, SDAS

def create_am_animation():
"""创建增材制造过程动画"""
print("\\n" + "="*60)
print("创建增材制造过程动画…")
print("="*60)

# 逐层沉积动画
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))

# 基板
Lx, Ly = 10e-3, 5e-3

# 初始化
ax1.set_xlim(0, Lx*1e3)
ax1.set_ylim(0, Ly*1e3)
ax1.set_xlabel('X (mm)', fontsize=11)
ax1.set_ylabel('Y (mm)', fontsize=11)
ax1.set_title('Layer-by-Layer Deposition', fontsize=12)
ax1.set_aspect('equal')

ax2.set_xlim(0, Lx*1e3)
ax2.set_ylim(0, Ly*1e3)
ax2.set_xlabel('X (mm)', fontsize=11)
ax2.set_ylabel('Y (mm)', fontsize=11)
ax2.set_title('Temperature Field', fontsize=12)
ax2.set_aspect('equal')

layers = []
temps = []

def init():
return []

def animate(frame):
ax1.clear()
ax2.clear()

ax1.set_xlim(0, Lx*1e3)
ax1.set_ylim(0, Ly*1e3)
ax1.set_xlabel('X (mm)', fontsize=11)
ax1.set_ylabel('Y (mm)', fontsize=11)
ax1.set_title('Layer-by-Layer Deposition', fontsize=12)
ax1.set_aspect('equal')

ax2.set_xlim(0, Lx*1e3)
ax2.set_ylim(0, Ly*1e3)
ax2.set_xlabel('X (mm)', fontsize=11)
ax2.set_ylabel('Y (mm)', fontsize=11)
ax2.set_title('Temperature Field', fontsize=12)
ax2.set_aspect('equal')

# 绘制已沉积层
n_layers = frame // 10 + 1
layer_thickness = 0.2 # mm

for i in range(n_layers):
y_bottom = i * layer_thickness
rect = Rectangle((0, y_bottom), Lx*1e3, layer_thickness,
linewidth=1, edgecolor='black', facecolor='steelblue', alpha=0.7)
ax1.add_patch(rect)

# 绘制当前扫描线
scan_progress = (frame % 10) / 10
ax1.plot([0, scan_progress * Lx*1e3],
[n_layers * layer_thickness – layer_thickness/2] * 2,
'r-', linewidth=3)

# 温度场（简化）
x_temp = np.linspace(0, Lx*1e3, 50)
y_temp = np.linspace(0, Ly*1e3, 50)
X_temp, Y_temp = np.meshgrid(x_temp, y_temp)

# 高斯温度分布
laser_x = scan_progress * Lx*1e3
laser_y = n_layers * layer_thickness – layer_thickness/2
T_field = 25 + 800 * np.exp(-((X_temp – laser_x)**2 + (Y_temp – laser_y)**2) / 10)

cf = ax2.contourf(X_temp, Y_temp, T_field, levels=10, cmap='hot', vmin=25, vmax=1000)

return []

anim = FuncAnimation(fig, animate, init_func=init, frames=50,
interval=100, blit=False)

writer = PillowWriter(fps=10)
anim.save(OUTPUT_DIR + 'am_deposition_animation.gif', writer=writer)
plt.close()

print(f"动画创建完成！已保存至 {OUTPUT_DIR}am_deposition_animation.gif")

# 主程序
if __name__ == "__main__":
# 运行所有案例
results1 = case1_melt_pool_temperature()
results2 = case2_scanning_strategies()
results3 = case3_residual_stress()
results4 = case4_solidification_microstructure()

# 创建动画
create_am_animation()

print("\\n" + "="*60)
print("所有案例完成！")
print("="*60)
print("\\n输出文件:")
print(" – case1_melt_pool_temperature.png")
print(" – case2_scanning_strategies.png")
print(" – case3_residual_stress.png")
print(" – case4_solidification_microstructure.png")
print(" – am_deposition_animation.gif")