第11章：数学优化与计算设计

本章深入探讨3D打印中的数学优化方法，从拓扑优化的理论基础到实际可制造性约束的处理。我们将建立严格的数学框架，理解各种优化算法的原理，并学习如何将优化结果转化为可打印的实际设计。通过掌握这些计算设计技术，您将能够创建轻量化、高性能的3D打印结构。

11.1 拓扑优化理论与SIMP方法

拓扑优化是在给定设计域内寻找最优材料分布的数学方法，广泛应用于轻量化设计。SIMP（Solid Isotropic Material with Penalization）方法是最成熟和应用最广的拓扑优化方法之一。

11.1.1 结构优化的数学框架

结构优化问题可以分为三个层次：尺寸优化、形状优化和拓扑优化。拓扑优化是最一般的形式，允许在设计域内任意分布材料。

标准的拓扑优化问题可以表述为：

$$\min_{\rho} c(\rho) = \mathbf{U}^T\mathbf{K}\mathbf{U} = \sum_{e=1}^{N} \rho_e^p \mathbf{u}_e^T \mathbf{k}_0 \mathbf{u}_e$$ 约束条件： $$\begin{cases} \mathbf{K}\mathbf{U} = \mathbf{F} \\ V(\rho) = \sum_{e=1}^{N} v_e \rho_e \leq V_0 \\ 0 < \rho_{min} \leq \rho_e \leq 1 \end{cases}$$ 其中：

$\rho_e$ 是单元 $e$ 的相对密度（设计变量）
$c(\rho)$ 是柔度（compliance），即结构柔性的度量
$\mathbf{K}$ 是全局刚度矩阵
$\mathbf{U}$ 是位移向量
$\mathbf{F}$ 是外力向量
$V_0$ 是允许的材料体积
$p$ 是惩罚因子（通常取3）

11.1.2 SIMP方法的核心思想

SIMP方法通过引入惩罚因子，使中间密度值（灰度单元）在优化过程中变得不利，从而推动解向0-1分布（黑白设计）收敛。

材料插值模型： $$E_e(\rho_e) = \rho_e^p E_0$$ 其中 $E_0$ 是实体材料的杨氏模量。当 $p > 1$ 时，中间密度的材料刚度相对其质量贡献被削弱，优化算法倾向于产生明确的实体或空洞区域。

密度-刚度关系图示（p=3）：
      E/E₀
      1.0 |                    ●
          |                 ●
      0.5 |            ●
          |        ●
      0.0 |●●●●●
          +--------------------→ ρ
          0    0.5    1.0

11.1.3 灵敏度分析与优化算法

优化求解需要计算目标函数对设计变量的灵敏度。对于柔度最小化问题： $$\frac{\partial c}{\partial \rho_e} = -p\rho_e^{p-1} \mathbf{u}_e^T \mathbf{k}_0 \mathbf{u}_e$$ 这个表达式表明，灵敏度与单元应变能成正比。物理意义是：在应变能高的区域增加材料能有效降低结构柔度。

常用的优化算法包括：

优化准则法（OC）：基于Kuhn-Tucker条件的启发式更新 $$\rho_e^{new} = \begin{cases} \max(\rho_{min}, \rho_e - m) & \text{if } \rho_e B_e^{\eta} \leq \max(\rho_{min}, \rho_e - m) \\ \min(1, \rho_e + m) & \text{if } \rho_e B_e^{\eta} \geq \min(1, \rho_e + m) \\ \rho_e B_e^{\eta} & \text{otherwise} \end{cases}$$
移动渐近线法（MMA）：更稳定的数学规划方法，适合处理多约束问题
水平集方法：隐式表示拓扑边界，自然满足0-1约束

11.1.4 过滤技术与数值稳定性

直接求解拓扑优化问题常遇到两个数值问题：

棋盘格现象：相邻单元密度交替变化
网格依赖性：细化网格产生完全不同的拓扑

密度过滤是最常用的解决方案： $$\tilde{\rho}_e = \frac{\sum_{i \in N_e} H_{ei} v_i \rho_i}{\sum_{i \in N_e} H_{ei} v_i}$$ 其中权重函数： $$H_{ei} = \max(0, r_{min} - \text{dist}(e, i))$$ 过滤半径 $r_{min}$ 控制最小特征尺寸，通常取2-3倍单元尺寸。

过滤示意图：
原始密度场          过滤后密度场
█ □ █ □           ▓ ▓ ▓ ▒
□ █ □ █    →      ▓ ▓ ▓ ▒
█ □ █ □           ▓ ▓ ▓ ▒
□ █ □ █           ▒ ▒ ▒ ░

11.1.5 投影方法与离散化

为获得清晰的0-1设计，可在密度过滤后应用投影函数： $$\bar{\rho}_e = \frac{\tanh(\beta \eta) + \tanh(\beta(\tilde{\rho}_e - \eta))}{\tanh(\beta \eta) + \tanh(\beta(1 - \eta))}$$ 其中 $\beta$ 控制投影陡度，$\eta$ 是阈值（通常取0.5）。随优化进程逐渐增大 $\beta$，实现从连续到离散的平滑过渡。

11.1.6 约束处理与多物理场

实际问题常包含多种约束：

应力约束：防止局部应力集中 $$\sigma_{vm,e} \leq \sigma_{yield} / SF$$ 挑战：应力是高度局部化的，需要大量约束（每个单元一个）
位移约束：限制特定点的变形 $$|u_i| \leq u_{max}$$
频率约束：避免共振 $$\omega_1 \geq \omega_{min}$$
屈曲约束：防止失稳 $$\lambda_{cr} \geq \lambda_{min}$$ 多物理场问题需要考虑热-力耦合： $$\min_{\rho} \alpha_1 c_{mech}(\rho) + \alpha_2 c_{thermal}(\rho)$$ 其中热柔度： $$c_{thermal} = \mathbf{T}^T\mathbf{K}_{thermal}\mathbf{T}$$

11.2 晶格结构设计与均质化

晶格结构通过周期性微结构实现优异的比强度和比刚度，是增材制造独特优势的体现。均质化理论提供了从微观到宏观的严格数学桥梁。

11.2.1 周期性微结构的数学描述

考虑具有特征长度 $\epsilon$ 的周期性结构，其位移场可展开为： $$\mathbf{u}^{\epsilon}(\mathbf{x}) = \mathbf{u}^0(\mathbf{x}) + \epsilon \mathbf{u}^1(\mathbf{x}, \mathbf{y}) + \epsilon^2 \mathbf{u}^2(\mathbf{x}, \mathbf{y}) + ...$$ 其中 $\mathbf{y} = \mathbf{x}/\epsilon$ 是快变量，描述微观尺度的变化。

单胞（unit cell）定义域： $$Y = [0, 1]^d, \quad d \in \{2, 3\}$$ 周期性边界条件： $$\mathbf{u}(\mathbf{y}) \text{ 在 } \partial Y \text{ 上满足周期性}$$

11.2.2 均质化方法

通过渐近展开和尺度分离，可导出均质化的弹性张量： $$C_{ijkl}^H = \frac{1}{|Y|} \int_Y C_{ijkl}(\mathbf{y}) \left( \delta_{ik}\delta_{jl} - \frac{\partial \chi^{kl}_i}{\partial y_j} \right) d\mathbf{y}$$ 其中 $\chi^{kl}$ 是特征位移场，满足单胞问题： $$\frac{\partial}{\partial y_j} \left( C_{ijmn}(\mathbf{y}) \frac{\partial \chi^{kl}_m}{\partial y_n} \right) = \frac{\partial C_{ijkl}}{\partial y_j} \quad \text{in } Y$$ 对于各向同性材料，均质化张量可简化为两个独立参数： $$C^H = \begin{bmatrix} \lambda + 2\mu & \lambda & 0 \\ \lambda & \lambda + 2\mu & 0 \\ 0 & 0 & \mu \end{bmatrix}$$

11.2.3 典型晶格类型与性能

常见晶格结构的相对密度-性能关系：

立方晶格（Cubic）： - 相对弹性模量：$E/E_s = \rho^2$ - 屈服强度：$\sigma_y/\sigma_{ys} = \rho^{1.5}$
八面体晶格（Octet-truss）： - 相对弹性模量：$E/E_s = \rho$ - 屈服强度：$\sigma_y/\sigma_{ys} = \rho$ - 拉伸主导，接近理论极限
Gyroid TPMS： - 隐式方程：$\sin(2\pi x/L)\cos(2\pi y/L) + \sin(2\pi y/L)\cos(2\pi z/L) + \sin(2\pi z/L)\cos(2\pi x/L) = t$ - 各向同性，光滑连续 - 相对弹性模量：$E/E_s \approx 0.3\rho^{1.6}$

晶格类型对比图：
          拉伸主导              弯曲主导
    Octet  Diamond        Cubic   Kelvin
E/Es  ρ     ρ^1.2         ρ^2     ρ^2.3
σy    ρ     ρ^1.2         ρ^1.5   ρ^1.8

11.2.4 梯度晶格设计

变密度晶格通过空间变化的相对密度实现性能梯度： $$\rho(\mathbf{x}) = \rho_{min} + (\rho_{max} - \rho_{min}) \cdot f(\mathbf{x})$$ 梁径调制（对于strut-based晶格）： $$d(\mathbf{x}) = d_0 \sqrt{\rho(\mathbf{x})/\rho_0}$$ TPMS厚度调制： $$t(\mathbf{x}) = t_0 \cdot g(\rho(\mathbf{x}))$$ 其中 $g$ 是通过数值拟合得到的映射函数。

11.2.5 多尺度优化

结合拓扑优化和晶格设计的多尺度方法：

两尺度优化问题： $$\min_{\rho, \mathbf{d}} c(\rho, \mathbf{d}) = \mathbf{U}^T\mathbf{K}(\rho, C^H(\mathbf{d}))\mathbf{U}$$ 其中：

$\rho$ 控制宏观材料分布
$\mathbf{d}$ 是微结构设计变量（如支柱直径、单胞尺寸）

并发优化策略：

固定微结构，优化宏观拓扑
固定拓扑，优化微结构参数
迭代直至收敛

11.2.6 晶格的失效模式

晶格结构的失效不仅包括材料屈服，还有特殊的失稳模式：

局部屈曲临界应力（对于细长支柱）： $$\sigma_{cr} = \frac{\pi^2 E I}{(KL)^2 A}$$ 其中 $K$ 是有效长度系数（取决于边界条件）。

整体失稳：需要考虑晶格的等效连续体模型，进行屈曲分析： $$(\mathbf{K} - \lambda \mathbf{K}_G)\boldsymbol{\phi} = 0$$ 疲劳性能：应力集中系数 $$K_t = 1 + 2\sqrt{a/\rho}$$ 其中 $a$ 是节点处的缺口半径，$\rho$ 是曲率半径。

11.3 形状优化与灵敏度分析

形状优化通过调整结构边界来改善性能，保持拓扑不变。这种方法特别适合精细调整已有设计，是拓扑优化后的重要后处理步骤。

11.3.1 形状参数化方法

形状优化的关键是选择合适的参数化方法来描述和控制边界变化。

1. 基于CAD参数的方法：直接使用CAD模型中的设计参数（如圆角半径、厚度等）： $$\mathbf{x}(\mathbf{p}) = \sum_{i=1}^n p_i \mathbf{f}_i(\mathbf{u}, \mathbf{v})$$ 其中 $\mathbf{p}$ 是设计参数向量，$\mathbf{f}_i$ 是基函数。

2. 基于网格节点的方法：直接移动边界节点： $$\mathbf{x}_{new} = \mathbf{x}_{old} + \delta \mathbf{x}$$ 需要平滑约束防止网格扭曲： $$\nabla^2 \delta \mathbf{x} = 0 \quad \text{(Laplacian smoothing)}$$

3. 基于基函数的方法：使用B样条或NURBS描述边界： $$\mathbf{x}(t) = \sum_{i=0}^n N_i^p(t) \mathbf{P}_i$$ 其中 $N_i^p$ 是p次B样条基函数，$\mathbf{P}_i$ 是控制点。

11.3.2 形状灵敏度分析

考虑一般目标函数 $J(\Omega, u(\Omega))$，其中 $\Omega$ 是设计域，$u$ 是状态变量（如位移场）。

物质导数法：形状导数可表示为边界积分： $$\frac{dJ}{d\Omega} = \int_{\partial \Omega} g(\mathbf{x}) V_n(\mathbf{x}) ds$$ 其中 $V_n$ 是边界法向速度，$g$ 是形状梯度密度。

对于柔度最小化： $$g = -\frac{1}{2} \boldsymbol{\sigma}:\boldsymbol{\varepsilon}$$ 伴随法：为避免计算每个设计变量的状态灵敏度，引入伴随方程： $$\mathbf{K}^T \boldsymbol{\lambda} = -\frac{\partial J}{\partial \mathbf{u}}$$ 则灵敏度： $$\frac{dJ}{dp_i} = \frac{\partial J}{\partial p_i} + \boldsymbol{\lambda}^T \frac{\partial \mathbf{K}}{\partial p_i} \mathbf{u}$$

11.3.3 速度场计算与正则化

从形状梯度到速度场需要求解扩展方程： $$-\alpha \nabla^2 \mathbf{V} + \mathbf{V} = -g \mathbf{n} \quad \text{on } \partial \Omega$$ $$\mathbf{V} = 0 \quad \text{in } \Omega$$ 其中 $\alpha$ 是正则化参数，控制速度场的光滑度。

Hilbertian方法：将 $L^2$ 内积替换为 $H^1$ 内积： $$\langle \mathbf{V}, \mathbf{W} \rangle_{H^1} = \int_{\partial \Omega} (\mathbf{V} \cdot \mathbf{W} + \epsilon^2 \nabla \mathbf{V} : \nabla \mathbf{W}) ds$$

11.3.4 网格变形技术

边界移动后需要更新内部网格，常用方法：

1. 弹性类比法：将网格视为弹性体： $$\nabla \cdot (E(\mathbf{x}) \nabla \mathbf{u}) = 0$$ 其中弹性模量与单元质量成反比： $$E = \frac{1}{V_e^p}$$

2. RBF插值： $$\mathbf{u}(\mathbf{x}) = \sum_{i=1}^n \alpha_i \phi(||\mathbf{x} - \mathbf{x}_i||) + \mathbf{p}(\mathbf{x})$$ 常用径向基函数：

Wendland函数：$\phi(r) = (1-r)^4_+(4r+1)$
薄板样条：$\phi(r) = r^2 \log r$

11.3.5 约束处理

形状优化常见约束：

1. 体积约束： $$V(\Omega) = \int_{\Omega} d\Omega \leq V_{max}$$ 灵敏度： $$\frac{dV}{d\Omega} = \int_{\partial \Omega} V_n ds$$

2. 周长约束（2D）或表面积约束（3D）： $$P(\partial \Omega) = \int_{\partial \Omega} ds \leq P_{max}$$ 灵敏度（2D）： $$\frac{dP}{d\Omega} = \int_{\partial \Omega} \kappa V_n ds$$ 其中 $\kappa$ 是曲率。

3. 最小/最大厚度约束：使用骨架法（medial axis）计算局部厚度： $$t_{min} \leq t(\mathbf{x}) \leq t_{max}$$

11.3.6 形状优化算法

1. 最速下降法： $$\Omega_{k+1} = \Omega_k - \alpha_k \nabla J(\Omega_k)$$ 步长 $\alpha_k$ 通过线搜索确定。

2. 拟牛顿法（L-BFGS）：近似Hessian矩阵，收敛更快： $$\mathbf{H}_{k+1} = \mathbf{H}_k + \text{rank-2 update}$$

3. 水平集方法：通过Hamilton-Jacobi方程演化水平集函数： $$\frac{\partial \phi}{\partial t} + V_n |\nabla \phi| = 0$$ 其中 $\phi = 0$ 定义边界。

11.4 多目标优化与Pareto前沿

在实际3D打印设计中，我们往往需要同时优化多个相互冲突的目标：重量、刚度、成本、打印时间等。多目标优化提供了系统化的方法来处理这些权衡关系。

11.4.1 多目标优化问题的数学表述

标准多目标优化问题（MOP）可表述为： $$\min_{\mathbf{x} \in \mathcal{X}} \mathbf{f}(\mathbf{x}) = [f_1(\mathbf{x}), f_2(\mathbf{x}), ..., f_m(\mathbf{x})]^T$$ 约束条件： $$\begin{cases} g_i(\mathbf{x}) \leq 0, \quad i = 1, ..., p \\ h_j(\mathbf{x}) = 0, \quad j = 1, ..., q \\ \mathbf{x}_L \leq \mathbf{x} \leq \mathbf{x}_U \end{cases}$$ 在3D打印优化中，典型的目标函数包括：

$f_1$：结构质量 $M = \int_{\Omega} \rho dV$
$f_2$：柔度 $C = \mathbf{U}^T\mathbf{K}\mathbf{U}$
$f_3$：最大应力 $\sigma_{max} = \max_{\Omega} \sigma_{vm}$
$f_4$：打印时间 $T = V_{material}/\dot{V} + N_{layers} \cdot t_{layer}$
$f_5$：支撑材料量 $V_{support}$

11.4.2 Pareto最优性理论

Pareto支配定义：解 $\mathbf{x}^1$ 支配 $\mathbf{x}^2$（记作 $\mathbf{x}^1 \prec \mathbf{x}^2$）当且仅当： $$\begin{cases} f_i(\mathbf{x}^1) \leq f_i(\mathbf{x}^2), \quad \forall i \in \{1, ..., m\} \\ \exists j \in \{1, ..., m\}: f_j(\mathbf{x}^1) < f_j(\mathbf{x}^2) \end{cases}$$ Pareto最优解：如果不存在 $\mathbf{x}' \in \mathcal{X}$ 使得 $\mathbf{x}' \prec \mathbf{x}^*$，则 $\mathbf{x}^*$ 是Pareto最优解。

Pareto前沿：所有Pareto最优解在目标空间的映射： $$PF = \{\mathbf{f}(\mathbf{x}^*) | \mathbf{x}^* \text{ is Pareto optimal}\}$$

双目标Pareto前沿示意图：
    f₂ ↑
       |  × 被支配解
    4  |× × 
       |  × ●---●
    3  | ×  ●     ●  ← Pareto前沿
       |   ●        ●
    2  |  ●          ●
       | ●            ●
    1  |●              ●
       +----------------→ f₁
       1  2  3  4  5  6

11.4.3 多目标优化方法

1. 权重法（Weighted Sum Method）：将多目标转化为单目标： $$f_{weighted} = \sum_{i=1}^m w_i f_i(\mathbf{x}), \quad \sum_{i=1}^m w_i = 1, w_i \geq 0$$ 局限性：无法获得非凸Pareto前沿的凹部分。

2. ε-约束法：优化一个主要目标，其他目标作为约束： $$\begin{cases} \min f_k(\mathbf{x}) \\ f_i(\mathbf{x}) \leq \epsilon_i, \quad i \neq k \end{cases}$$ 通过系统变化 $\epsilon_i$ 可获得完整Pareto前沿。

3. 目标规划法：最小化与理想点的距离： $$\min_{\mathbf{x}} \left( \sum_{i=1}^m \left( \frac{f_i(\mathbf{x}) - f_i^{ideal}}{f_i^{nadir} - f_i^{ideal}} \right)^p \right)^{1/p}$$ 其中：

$f_i^{ideal}$：第i个目标的理想值（单独优化时的最优值）
$f_i^{nadir}$：第i个目标的最差值（在Pareto前沿上）
$p$：范数参数（$p=2$ 对应欧氏距离）

4. 法向边界交叉法（NBI）：通过构造均匀分布的法向量来获得均匀的Pareto前沿： $$\begin{cases} \max_{\mathbf{x}, t} t \\ \boldsymbol{\Phi}\mathbf{w} + t\mathbf{n} = \mathbf{f}(\mathbf{x}) \end{cases}$$ 其中 $\boldsymbol{\Phi}$ 是理想点矩阵，$\mathbf{n}$ 是法向量。

11.4.4 进化算法在多目标优化中的应用

进化算法特别适合求解多目标优化问题，因为它们维护一个解的种群，可同时逼近整个Pareto前沿。

NSGA-II算法核心机制：

非支配排序：将种群分层，第k层的解不被前k-1层的任何解支配。
拥挤度距离：衡量解在目标空间的分布密度： $$CD_i = \sum_{m=1}^M \frac{f_m^{i+1} - f_m^{i-1}}{f_m^{max} - f_m^{min}}$$
选择算子：优先选择低支配等级的解，同等级内选择拥挤度大的解。

MOEA/D（基于分解的多目标进化算法）：将MOP分解为N个标量优化子问题： $$g^{te}(\mathbf{x}|\boldsymbol{\lambda}, \mathbf{z}^*) = \max_{1 \leq i \leq m} \{\lambda_i |f_i(\mathbf{x}) - z_i^*|\}$$ 其中 $\boldsymbol{\lambda}$ 是权重向量，$\mathbf{z}^*$ 是参考点。

11.4.5 3D打印特定的多目标权衡

1. 强度-重量权衡：对于晶格填充结构： $$\begin{cases} f_1 = \rho_{rel} \cdot V_{total} \\ f_2 = \sigma_{max} / \sigma_{yield} \end{cases}$$ Pareto前沿通常遵循幂律关系： $$\sigma/\sigma_s \propto (\rho/\rho_s)^n$$ 其中 $n \in [1, 2]$ 取决于晶格类型。

2. 精度-速度权衡： $$\begin{cases} f_1 = \sqrt{\sum_{i} (x_i^{actual} - x_i^{nominal})^2} \\ f_2 = 1/t_{print} \end{cases}$$ 提高打印速度通常导致振动增加，降低精度： $$\sigma_{position} \propto v^2/a_{max}$$

3. 支撑-表面质量权衡： $$\begin{cases} f_1 = V_{support} \\ f_2 = \sum_{facets} A_i \cdot (1 - \mathbf{n}_i \cdot \mathbf{z})^+ \end{cases}$$ 其中第二项衡量需要支撑的悬垂面积。

11.4.6 交互式决策支持

获得Pareto前沿后，需要决策者选择最终解。常用方法：

1. 可视化技术： - 平行坐标图：适合高维目标空间 - 散点矩阵：显示目标间的两两关系 - 雷达图：直观显示单个解的多目标性能

2. 偏好引导：引入决策者偏好信息： $$u(\mathbf{f}) = \prod_{i=1}^m \left(1 - e^{-\alpha_i(f_i^{max} - f_i)}\right)$$ 其中 $\alpha_i$ 反映对第i个目标的重视程度。

3. 模糊决策：使用隶属函数量化满意度： $$\mu_i(f_i) = \begin{cases} 1 & f_i \leq f_i^{best} \\ \frac{f_i^{worst} - f_i}{f_i^{worst} - f_i^{best}} & f_i^{best} < f_i < f_i^{worst} \\ 0 & f_i \geq f_i^{worst} \end{cases}$$ 综合满意度： $$\mu_{total} = \min_{i} \mu_i(f_i) \quad \text{(保守策略)}$$ 或 $$\mu_{total} = \prod_{i} \mu_i(f_i)^{w_i} \quad \text{(折衷策略)}$$

11.5 有限元分析集成

有限元分析（FEA）是验证和优化3D打印设计的关键工具。将FEA集成到设计流程中，可以在打印前预测结构性能，避免昂贵的试错过程。

11.5.1 3D打印结构的有限元建模

3D打印结构具有独特的特征，需要特殊的建模考虑：

1. 各向异性材料模型： FDM打印件在不同方向具有不同的力学性能。正交各向异性本构关系： $$\begin{bmatrix} \sigma_{11} \\ \sigma_{22} \\ \sigma_{33} \\ \tau_{23} \\ \tau_{13} \\ \tau_{12} \end{bmatrix} = \begin{bmatrix} C_{11} & C_{12} & C_{13} & 0 & 0 & 0 \\ C_{12} & C_{22} & C_{23} & 0 & 0 & 0 \\ C_{13} & C_{23} & C_{33} & 0 & 0 & 0 \\ 0 & 0 & 0 & C_{44} & 0 & 0 \\ 0 & 0 & 0 & 0 & C_{55} & 0 \\ 0 & 0 & 0 & 0 & 0 & C_{66} \end{bmatrix} \begin{bmatrix} \varepsilon_{11} \\ \varepsilon_{22} \\ \varepsilon_{33} \\ \gamma_{23} \\ \gamma_{13} \\ \gamma_{12} \end{bmatrix}$$ 其中：

方向1：沿打印路径方向
方向2：垂直于路径在层内
方向3：垂直于层（Z方向）

典型的强度比例（以PLA为例）： $$E_1 : E_2 : E_3 \approx 1 : 0.95 : 0.7$$ $$\sigma_{y1} : \sigma_{y2} : \sigma_{y3} \approx 1 : 0.9 : 0.5$$

2. 层间结合建模：使用内聚力单元（cohesive elements）模拟层间界面： $$\mathbf{T} = \begin{bmatrix} T_n \\ T_s \\ T_t \end{bmatrix} = \mathbf{K} \begin{bmatrix} \delta_n \\ \delta_s \\ \delta_t \end{bmatrix}$$ 双线性牵引-分离法则： $$T = \begin{cases} K \delta & \delta < \delta_0 \\ K \delta_0 \frac{\delta_f - \delta}{\delta_f - \delta_0} & \delta_0 \leq \delta < \delta_f \\ 0 & \delta \geq \delta_f \end{cases}$$ 其中 $\delta_0$ 是损伤起始位移，$\delta_f$ 是完全失效位移。

3. 填充结构的等效模型：对于规则填充图案，使用均质化得到的等效属性： $$E_{eff} = E_{solid} \cdot \rho_{infill}^n$$ 其中指数 $n$ 取决于填充类型：

蜂窝填充：$n \approx 2.5$
直线填充：$n \approx 1.8$
三角填充：$n \approx 2.2$

11.5.2 网格生成策略

1. 体素化网格：直接从切片数据生成六面体网格：

切片层 → 2D像素化 → 垂直堆叠 → 3D体素

优点：完美匹配打印分辨率缺点：单元数量巨大

2. 自适应网格细化：在高应力梯度区域细化网格： $$h_{new} = h_{old} \cdot \left( \frac{\eta_{target}}{\eta_{element}} \right)^{1/(p+1)}$$ 其中 $\eta$ 是误差估计，$p$ 是单元阶数。

3. 多尺度网格： - 宏观尺度：整体几何 - 细观尺度：填充结构 - 微观尺度：层间界面

使用子模型技术或多尺度有限元方法连接不同尺度。

11.5.3 边界条件与载荷

残余应力初始化： 3D打印过程产生的残余应力： $$\sigma_{residual} = E \alpha \Delta T$$ 其中 $\alpha$ 是热膨胀系数，$\Delta T$ 是冷却温差。

对于翘曲变形，使用简化的层合板理论： $$\kappa = \frac{6(\alpha_1 - \alpha_2)\Delta T}{h(1 + 4m + 6m^2 + 4m^3 + m^4)}$$ 其中 $m = E_1/E_2$ 是模量比。

打印方向相关的强度准则： Tsai-Wu失效准则适合各向异性材料： $$F_i \sigma_i + F_{ij} \sigma_i \sigma_j = 1$$ 其中： $$F_1 = \frac{1}{X_t} - \frac{1}{X_c}, \quad F_{11} = \frac{1}{X_t X_c}$$ $X_t$、$X_c$ 分别是拉伸和压缩强度。

11.5.4 求解策略与收敛性

1. 非线性求解：对于大变形或材料非线性：

Newton-Raphson迭代： $$\mathbf{K}_T^{(i)} \Delta \mathbf{u}^{(i)} = \mathbf{F}_{ext} - \mathbf{F}_{int}^{(i)}$$ 其中 $\mathbf{K}_T$ 是切线刚度矩阵。

收敛准则： $$\frac{||\mathbf{R}^{(i)}||}{||\mathbf{F}_{ext}||} < \epsilon_F \quad \text{and} \quad \frac{||\Delta \mathbf{u}^{(i)}||}{||\mathbf{u}^{(i)}||} < \epsilon_u$$

2. 接触问题：装配体中的接触使用罚函数法： $$\mathbf{F}_{contact} = k_n g_n \mathbf{n}$$ 其中 $k_n$ 是接触刚度，$g_n$ 是间隙函数。

3. 动态分析：对于振动分析，考虑材料阻尼： $$\mathbf{M}\ddot{\mathbf{u}} + \mathbf{C}\dot{\mathbf{u}} + \mathbf{K}\mathbf{u} = \mathbf{F}(t)$$ Rayleigh阻尼： $$\mathbf{C} = \alpha\mathbf{M} + \beta\mathbf{K}$$

11.5.5 后处理与结果解释

1. 应力集中因子：识别潜在失效位置： $$K_t = \frac{\sigma_{max}}{\sigma_{nominal}}$$ 对于3D打印件，层间和填充边界常出现应力集中。

2. 安全系数分布：考虑各向异性的安全系数： $$SF = \min \left( \frac{X_t}{\sigma_1}, \frac{X_c}{|\sigma_1|}, \frac{Y_t}{\sigma_2}, \frac{S}{\tau_{12}} \right)$$

3. 疲劳寿命预测：使用修正的S-N曲线： $$N = A \cdot (\Delta \sigma)^{-b} \cdot C_{surface} \cdot C_{size} \cdot C_{load}$$ 其中修正系数考虑3D打印特有的表面粗糙度和尺寸效应。

11.5.6 优化循环中的FEA

1. 灵敏度分析集成：直接微分法计算目标函数对设计变量的导数： $$\frac{d f}{d x_i} = \frac{\partial f}{\partial x_i} + \frac{\partial f}{\partial \mathbf{u}} \frac{d\mathbf{u}}{d x_i}$$ 通过伴随法避免计算 $d\mathbf{u}/dx_i$： $$\boldsymbol{\lambda}^T \mathbf{K} = \frac{\partial f}{\partial \mathbf{u}}$$ 则： $$\frac{d f}{d x_i} = \frac{\partial f}{\partial x_i} - \boldsymbol{\lambda}^T \frac{\partial \mathbf{K}}{\partial x_i} \mathbf{u}$$

2. 代理模型加速：使用Kriging或神经网络替代昂贵的FEA：

Kriging预测： $$\hat{y}(\mathbf{x}) = \mu + \mathbf{r}^T(\mathbf{x}) \mathbf{R}^{-1}(\mathbf{y} - \mathbf{1}\mu)$$ 其中 $\mathbf{R}$ 是相关矩阵，$\mathbf{r}$ 是相关向量。

预测不确定性： $$s^2(\mathbf{x}) = \sigma^2 \left(1 - \mathbf{r}^T \mathbf{R}^{-1} \mathbf{r} + \frac{(1 - \mathbf{1}^T \mathbf{R}^{-1} \mathbf{r})^2}{\mathbf{1}^T \mathbf{R}^{-1} \mathbf{1}} \right)$$

3. 多保真度方法：结合低保真（粗网格）和高保真（细网格）模型： $$\hat{f}_{HF}(\mathbf{x}) = \rho(\mathbf{x}) f_{LF}(\mathbf{x}) + \delta(\mathbf{x})$$ 其中 $\rho$ 是缩放因子，$\delta$ 是加性修正。

11.6 优化结果的可制造性处理

拓扑优化和形状优化的结果往往需要后处理才能满足3D打印的制造约束。本节讨论如何将理论优化结果转化为实际可打印的设计。

11.6.1 最小特征尺寸控制

优化结果可能产生过细的结构特征，无法可靠打印。控制方法包括：

1. 密度投影与形态学操作：使用腐蚀-膨胀序列消除细小特征：

腐蚀操作： $$\tilde{\rho}_{erode} = \frac{\tanh(\beta \eta) + \tanh(\beta(\hat{\rho} - \eta))}{\tanh(\beta \eta) + \tanh(\beta(1 - \eta))}$$ 膨胀操作： $$\tilde{\rho}_{dilate} = \frac{\tanh(\beta \eta) + \tanh(\beta(\tilde{\rho}_{erode} + \Delta - \eta))}{\tanh(\beta \eta) + \tanh(\beta(1 - \eta))}$$ 其中 $\Delta$ 控制最小特征尺寸： $$l_{min} = 2\Delta \cdot l_{element}$$

2. 鲁棒公式化：同时考虑腐蚀、中间和膨胀设计： $$\min_{\rho} \max(f_{erode}, f_{intermediate}, f_{dilate})$$ 这确保结构在制造误差下仍保持性能。

3. 显式几何约束：使用骨架提取算法计算局部厚度： $$d_{medial}(\mathbf{x}) = \min_{\mathbf{y} \in \partial \Omega} ||\mathbf{x} - \mathbf{y}||$$ 约束： $$d_{medial}(\mathbf{x}) \geq r_{min}, \quad \forall \mathbf{x} \in \Omega$$

11.6.2 悬垂角约束

FDM打印需要控制悬垂角以减少支撑需求：

1. 连续悬垂约束：定义可打印性场： $$p(\mathbf{x}) = \int_0^{z(\mathbf{x})} \rho(\mathbf{x}', y, z') \cdot \exp(-\lambda ||\mathbf{x}' - \mathbf{x}||) dz'$$ 其中 $\lambda$ 控制支撑锥角： $$\theta_{overhang} = \arctan(1/\lambda)$$ 约束： $$\rho(\mathbf{x}) \leq p(\mathbf{x}) + \epsilon$$

2. 离散层方法：逐层施加支撑约束： $$\rho_{i,j,k} \leq \max(\rho_{i-1,j,k-1}, \rho_{i,j,k-1}, \rho_{i+1,j,k-1})$$ 这确保每个体素都有下方支撑。

3. 打印方向优化：同时优化结构和打印方向： $$\min_{\rho, \mathbf{d}} f(\rho) + \alpha V_{support}(\rho, \mathbf{d})$$ 其中 $\mathbf{d}$ 是打印方向单位向量。

11.6.3 网格到实体转换

将优化的密度场或网格转换为可打印的CAD模型：

1. 等值面提取：使用Marching Cubes算法提取 $\rho = 0.5$ 等值面：

插值计算顶点位置： $$\mathbf{v} = \mathbf{v}_1 + \frac{\rho_{threshold} - \rho_1}{\rho_2 - \rho_1}(\mathbf{v}_2 - \mathbf{v}_1)$$

2. 网格平滑： Laplacian平滑减少阶梯效应： $$\mathbf{v}_{new} = \mathbf{v}_{old} + \lambda \sum_{i \in N(v)} w_i (\mathbf{v}_i - \mathbf{v}_{old})$$ 保持体积的Taubin平滑： $$\begin{cases} \mathbf{v}^{(k+1/2)} = \mathbf{v}^{(k)} + \lambda \mathbf{L}\mathbf{v}^{(k)} \\ \mathbf{v}^{(k+1)} = \mathbf{v}^{(k+1/2)} + \mu \mathbf{L}\mathbf{v}^{(k+1/2)} \end{cases}$$ 其中 $\lambda > 0$，$\mu < -\lambda$。

3. 网格简化：使用二次误差度量（QEM）减少三角形数量：

边收缩代价： $$\mathbf{Q} = \sum_{faces} \mathbf{n}_f \mathbf{n}_f^T$$ $$cost = \mathbf{v}^T \mathbf{Q} \mathbf{v}$$

11.6.4 分段与装配设计

大型优化结构可能超出打印尺寸，需要分段：

1. 自动分段算法：基于应力的分段：在低应力区域切割 $$\min_{\Pi} \sum_{cut} \int_{\Gamma_{cut}} \sigma_{vm} dA$$ 其中 $\Pi$ 是分段方案，$\Gamma_{cut}$ 是切割面。

2. 连接设计：优化连接刚度和强度：

螺栓连接的等效刚度： $$k_{joint} = \frac{1}{1/k_{bolt} + 1/k_{member1} + 1/k_{member2}}$$ 卡扣连接的保持力： $$F_{retention} = \sigma_{yield} \cdot A_{undercut} \cdot \tan(\alpha)$$ 其中 $\alpha$ 是卡扣角度。

3. 公差分配：统计公差叠加： $$\sigma_{assembly}^2 = \sum_{i} \left(\frac{\partial f}{\partial x_i}\right)^2 \sigma_i^2$$ 对于3D打印典型公差：

FDM：±0.2-0.5mm
SLA：±0.05-0.15mm
SLS：±0.1-0.3mm

11.6.5 多材料映射

将连续材料分布映射到离散材料选择：

1. 聚类方法： K-means聚类将连续属性离散化： $$\min_{\mathcal{C}} \sum_{i=1}^k \sum_{\mathbf{x} \in C_i} ||\mathbf{E}(\mathbf{x}) - \boldsymbol{\mu}_i||^2$$ 其中 $\mathbf{E}$ 是弹性张量，$\boldsymbol{\mu}_i$ 是第i类材料的平均属性。

2. 界面过渡设计：梯度过渡区减少应力集中： $$E(x) = E_1 + (E_2 - E_1) \cdot \frac{1}{2}\left(1 + \tanh\left(\frac{x - x_0}{w}\right)\right)$$ 其中 $w$ 控制过渡区宽度。

3. 打印路径规划：多材料切换优化： $$\min \sum_{switches} t_{purge} + l_{travel}/v_{travel}$$ 约束材料切换频率以减少浪费。

11.6.6 验证与迭代

1. 可打印性检查清单： - 最小壁厚 ≥ 2×喷嘴直径 - 最小孔径 ≥ 3×层厚 - 悬垂角 ≤ 45°（无支撑） - 桥接跨度 ≤ 50×喷嘴直径 - 长宽比 ≤ 10:1（防止翘曲）

2. 打印仿真：基于G-code的打印过程仿真：

材料沉积路径
层间等待时间
热历史预测
支撑生成验证

3. 试验验证：缩放模型测试： $$\sigma_{model} = \sigma_{full} \cdot \left(\frac{L_{model}}{L_{full}}\right)^{-n}$$

其中 $n$ 是尺寸效应指数（通常0.05-0.1）。

本章小结

本章深入探讨了3D打印中的数学优化方法，从拓扑优化的SIMP方法到多目标优化的Pareto前沿理论，从晶格结构的均质化分析到有限元分析的集成应用。关键要点包括：

拓扑优化提供了在给定约束下寻找最优材料分布的系统方法，SIMP方法通过惩罚中间密度推动解向0-1分布收敛
晶格结构通过周期性微结构实现优异的比强度，均质化理论提供了从微观到宏观的数学桥梁
形状优化通过调整边界来精细调整性能，需要careful的参数化和灵敏度分析
多目标优化处理实际设计中的多个冲突目标，Pareto前沿提供了权衡决策的理论基础
有限元分析验证优化设计，需要考虑3D打印特有的各向异性和层间结合
可制造性处理将理论优化结果转化为实际可打印设计，包括特征尺寸控制、悬垂约束和分段装配

掌握这些优化技术，您将能够设计出既满足性能要求又适合3D打印制造的创新结构。

练习题

基础题

11.1 给定一个2D悬臂梁问题，长100mm，高40mm，左端固定，右端施加向下100N集中力。使用SIMP方法，体积分数限制为0.3，惩罚因子p=3。计算第一次迭代后，距固定端20mm、距底边20mm处单元的灵敏度值。假设初始密度均匀分布为0.3，该点应变能密度为1.5 J/mm³。

提示

灵敏度公式：$\frac{\partial c}{\partial \rho_e} = -p\rho_e^{p-1} u_e^T k_0 u_e$

答案

灵敏度 = -3 × 0.3² × 1.5 = -0.405 J/mm³ 负值表示增加该处密度会降低柔度（提高刚度）

11.2 一个Octet-truss晶格结构，相对密度为0.2，基体材料杨氏模量为70 GPa（铝合金）。计算该晶格的等效杨氏模量和比刚度（E/ρ）。如果要达到等效模量14 GPa，需要多大的相对密度？

提示

Octet-truss：$E/E_s = \rho$，拉伸主导结构

答案

等效模量：E = 70 × 0.2 = 14 GPa 比刚度：E/ρ = 14/(0.2×2.7) = 25.9 GPa/(g/cm³) 要达到14 GPa：ρ = 14/70 = 0.2（恰好相同）

11.3 设计一个双目标优化问题，目标1：最小化质量 $f_1 = \rho V$，目标2：最小化柔度 $f_2 = c/c_0$。已知三个设计点：A(0.3, 2.5)，B(0.5, 1.8)，C(0.7, 1.2)。判断哪些点在Pareto前沿上，并计算点B相对于理想点(0.3, 1.2)的加权切比雪夫距离（权重相等）。

提示

检查Pareto支配关系；切比雪夫距离 = max{w₁|f₁-f₁|, w₂|f₂-f₂|}

答案

Pareto前沿：所有三点都在（无相互支配）点B的切比雪夫距离 = max{0.5×|0.5-0.3|, 0.5×|1.8-1.2|} = max{0.1, 0.3} = 0.3

挑战题

11.4 考虑一个变密度TPMS结构，使用Gyroid函数：$f(x,y,z) = \sin(2\pi x/L)\cos(2\pi y/L) + \sin(2\pi y/L)\cos(2\pi z/L) + \sin(2\pi z/L)\cos(2\pi x/L) - t(x,y,z)$。如果要实现密度从中心ρ=0.4线性过渡到边缘ρ=0.1的球形梯度（半径R=50mm），推导厚度参数t(r)与径向距离r的关系。已知对于均匀Gyroid，ρ ≈ 0.5(1 - 0.3t)当|t|<1。

提示

建立ρ(r)函数，然后反解t(ρ)，最后得到t(r)的复合函数

答案

密度分布：ρ(r) = 0.4 - 0.3r/50 （r从0到50mm）从ρ ≈ 0.5(1 - 0.3t)反解：t = (1 - 2ρ)/0.3 代入：t(r) = (1 - 2(0.4 - 0.3r/50))/0.3 = (0.2 + 0.6r/50)/0.3 = 0.67 + 0.04r 验证：r=0时t=0.67→ρ=0.4✓；r=50时t=2.67（超出有效范围，需要修正模型）

11.5 优化一个3D打印支架，同时考虑三个目标：f₁=质量(kg)，f₂=最大应力(MPa)，f₃=打印时间(小时)。使用ε-约束法，主目标为最小化应力，质量约束≤0.5kg，打印时间约束≤3小时。已知Pareto最优解集中有一解的目标值为(0.45, 25, 2.8)，灵敏度为∇f₂=[-15, 1, -8]ᵀ。如果放松质量约束到0.48kg，估计应力的改善量。

提示

使用灵敏度分析和拉格朗日乘子理论

答案

质量约束放松：Δε₁ = 0.48 - 0.45 = 0.03 kg 根据KKT条件，在最优点：∇f₂ + λ₁∇g₁ + λ₂∇g₂ = 0 其中g₁ = f₁ - ε₁（质量约束），g₂ = f₃ - ε₃（时间约束）由于约束放松，一阶近似：Δf₂ ≈ (∂f₂/∂ε₁)Δε₁ = -λ₁Δε₁ 从灵敏度：λ₁ ≈ 15 MPa/kg 应力改善：Δf₂ ≈ -15 × 0.03 = -0.45 MPa 新应力约为24.55 MPa

11.6 设计一个自支撑的拓扑优化结构，打印方向为+Z。使用连续悬垂约束，支撑锥角45°。给定密度场在某点(x,y,z)=(10,10,20)mm处ρ=0.8，计算该点所需的最小支撑密度场p值。假设支撑区域为以该点为顶点向下的圆锥，使用指数衰减核$\exp(-\lambda r)$，其中λ=1对应45°锥角。

提示

计算45°锥内的密度积分，考虑圆锥几何

答案

45°锥角→水平投影半径=垂直高度在高度z'处，积分区域半径r_max = z - z' 支撑场：p = ∫₀²⁰ ∫₀^(2π) ∫₀^(z-z') ρ(r,θ,z')exp(-r)r dr dθ dz' 假设下方均匀密度ρ₀： p = 2πρ₀ ∫₀²⁰ ∫₀^(20-z') r exp(-r) dr dz' = 2πρ₀ ∫₀²⁰ [1-(1+r)exp(-r)]|₀^(20-z') dz' ≈ 2πρ₀ × 20 × 0.865 = 109ρ₀ 要满足ρ=0.8≤p，需要ρ₀ ≥ 0.0073 实际中通常需要更高的支撑密度（~0.1-0.3）

常见陷阱与错误

优化不收敛：检查步长、惩罚参数是否合适；考虑使用continuation方法逐渐增加惩罚
棋盘格现象：使用密度过滤或灵敏度过滤；确保过滤半径大于单元尺寸
局部极小值：使用多起始点或全局优化算法；考虑问题的凸性
网格依赖性：使用正则化技术；控制最小特征尺寸
制造约束违反：在优化过程中显式考虑约束，而非后处理
多目标权重选择：进行敏感性分析；使用自适应权重方法
FEA计算成本：使用代理模型或降阶模型；并行计算
优化结果不可打印：增加鲁棒性考虑；验证关键尺寸

最佳实践检查清单

优化前准备

[ ] 明确定义设计域和非设计域
[ ] 确定所有载荷工况和边界条件
[ ] 选择appropriate目标函数和约束
[ ] 考虑制造约束（最小尺寸、悬垂角等）
[ ] 准备验证模型或实验数据

优化设置

[ ] 选择合适的优化方法（SIMP、Level-set、BESO等）
[ ] 设置reasonable收敛准则
[ ] 使用continuation策略（渐进加载参数）
[ ] 实施过滤技术防止数值不稳定
[ ] 考虑多尺度（如果涉及晶格结构）

求解过程

[ ] 监控目标函数和约束的收敛历史
[ ] 检查灵敏度计算的准确性
[ ] 定期保存中间结果
[ ] 验证物理合理性（如对称性）
[ ] 注意异常的密度分布

后处理验证

[ ] 执行详细FEA验证优化结果
[ ] 检查所有制造约束满足情况
[ ] 生成光滑的可打印几何
[ ] 评估对参数变化的鲁棒性
[ ] 考虑分段和装配（如需要）

制造准备

[ ] 导出合适的文件格式（STL、STEP等）
[ ] 添加必要的加工余量
[ ] 设计连接和装配特征
[ ] 生成支撑结构（如需要）
[ ] 准备打印参数建议

通过系统应用这些优化方法，您可以充分发挥3D打印的设计自由度，创造出传统制造难以实现的高性能结构。记住，优化只是工具，工程判断和实际验证同样重要。