主题007：热应力优化设计与可靠性分析

目录

1. 引言
2. 核心理论
3. 实例一：热应力约束优化
4. 实例二：可靠性分析
5. 实例三：拓扑优化
6. 实例四：多目标优化
7. 实例五：可靠性优化设计RBDO
8. 总结与习题

---

引言

1.1 工程背景

在现代工程设计中，结构不仅要承受机械载荷，还必须应对复杂的热环境。航空航天器再入大气层时表面温度可达上千摄氏度，核反应堆压力容器在运行中经历剧烈的温度循环，汽车发动机部件在启动-停止过程中承受反复的热冲击。这些热载荷往往比机械载荷更加危险，因为热应力具有全局性（影响整个结构）、耦合性（与材料性能相互影响）和时变性（随温度场演化）的特点。

传统设计方法采用安全系数来应对不确定性，但这种做法存在明显缺陷：

• 过度保守：统一的安全系数无法反映不同失效模式的实际风险
• 经济性差：为追求安全而过度设计，造成材料浪费
• 可靠性未知：无法量化评估设计的实际可靠度

优化设计与可靠性分析的引入为解决这些问题提供了科学途径。通过数学规划方法寻找最优设计方案，同时利用概率统计理论量化不确定性，工程师可以在性能、成本和可靠性之间取得最佳平衡。

1.2 学习目标

完成本主题学习后，读者将能够：

1. 建立热应力优化问题的数学模型
2. 使用Python实现尺寸优化、形状优化和拓扑优化
3. 运用蒙特卡洛模拟和FORM方法进行可靠性分析
4. 求解多目标优化问题并理解Pareto最优概念
5. 掌握基于可靠性的优化设计（RBDO）方法
6. 将优化理论应用于实际工程问题

1.3 本章结构

本章包含五个循序渐进的实例：

• 实例一：从简单的尺寸优化入手，建立热应力约束优化的基本框架
• 实例二：引入不确定性，学习可靠性分析的核心方法
• 实例三：拓展到拓扑优化，探索结构构型的最优设计
• 实例四：处理多目标冲突，掌握Pareto优化方法
• 实例五：综合优化与可靠性，实现RBDO设计

---

核心理论

2.1 优化设计基础

2.1.1 优化问题的数学表述

一个典型的优化问题可以表述为：

$$\begin{array}{rlrl}& \text{minimize}& & f(\mathbf{x})\\ & \text{subject to}& & g_i(\mathbf{x})\le 0,i=1,2,\dots ,m\\ & & & h_j(\mathbf{x})=0,j=1,2,\dots ,p\\ & & & x_k^{\text{min}}\le x_k\le x_k^{\text{max}},k=1,2,\dots ,n\end{array}$$

其中：

• $\mathbf{x}=[x_1,x_2,\dots ,x_n{]}^T$ 是设计变量向量
• $f(\mathbf{x})$ 是目标函数（如重量、成本、应力等）
• $g_i(\mathbf{x})$ 是不等式约束（如应力约束、变形约束）
• $h_j(\mathbf{x})$ 是等式约束
• $x_k^{\text{min}}$ 和 $x_k^{\text{max}}$ 是设计变量的边界

2.1.2 优化算法分类

基于梯度的方法：

• 序列线性规划（SLP）：将非线性问题线性化后求解
• 序列二次规划（SQP）：构建二次近似模型，收敛速度快
• 可行方向法：沿可行方向搜索，保证迭代点可行

无梯度方法：

• 遗传算法（GA）：模拟自然选择，适用于离散变量
• 粒子群优化（PSO）：模拟鸟群觅食，并行搜索能力强
• 模拟退火（SA）：模拟物理退火过程，可跳出局部最优

2.2 热应力优化理论

2.2.1 热弹性本构关系

考虑热效应的线弹性本构方程为：

$${\sigma }_{ij}=C_{ijkl}({\epsilon }_{kl}-{\alpha }_{kl}\Delta T)$$

其中：

• ${\sigma }_{ij}$ 是应力张量
• $C_{ijkl}$ 是弹性刚度张量
• ${\epsilon }_{kl}$ 是总应变张量
• ${\alpha }_{kl}$ 是热膨胀系数张量
• $\Delta T$ 是温度变化

对于各向同性材料，简化为：

$${\sigma }_{ij}=\lambda {\epsilon }_{kk}{\delta }_{ij}+2\mu {\epsilon }_{ij}-(3\lambda +2\mu )\alpha \Delta T{\delta }_{ij}$$

2.2.2 热应力约束的数学表达

热应力优化中常见的约束包括：

应力约束：
$$g(\mathbf{x})=\frac{{\sigma }_{\text{max}}(\mathbf{x})}{{\sigma }_{\text{allow}}}-1\le 0$$

变形约束：
$$g(\mathbf{x})=\frac{w_{\text{max}}(\mathbf{x})}{w_{\text{allow}}}-1\le 0$$

热应变能约束：
$$g(\mathbf{x})=\frac{U_{\text{thermal}}(\mathbf{x})}{U_{\text{allow}}}-1\le 0$$

2.3 拓扑优化理论

2.3.1 SIMP方法

固体各向同性材料惩罚（SIMP）方法是目前最广泛使用的拓扑优化方法。其核心思想是引入伪密度 ${\rho }_e\in [0,1]$ 表示单元e的材料分布：

$$E_e={\rho }_e^p{E}_0$$

其中：

• $E_e$ 是单元的弹性模量
• $E_0$ 是实体材料的弹性模量
• $p\ge 3$ 是惩罚因子（通常取3）

惩罚因子的作用是使中间密度（$0<\rho <1$）的材料变得不经济，从而推动优化向0/1的离散解收敛。

2.3.2 敏度分析

拓扑优化的核心是计算目标函数对设计变量的敏度（梯度）。对于热应变能最小化问题：

$$C=\frac{1}{2}{\mathbf{u}}^T\mathbf{Ku}$$

敏度为：

$$\frac{\partial C}{\partial {\rho }_e}=-\frac{p}{2}{\rho }_e^{p-1}{\mathbf{u}}_e^T{\mathbf{k}}_0{\mathbf{u}}_e$$

其中${\mathbf{k}}_0$是实体单元的刚度矩阵。

2.3.3 密度过滤

为避免棋盘格现象和网格依赖性，需要对敏度进行过滤：

$$\frac{\widehat{\partial C}}{\partial {\rho }_e}=\frac{{\sum }_{i\in N_e}w(r_{ei})\frac{\partial C}{\partial {\rho }_i}}{{\sum }_{i\in N_e}w(r_{ei})}$$

其中$w(r)$是权重函数，$N_e$是以单元e为中心的邻域。

2.4 可靠性分析理论

2.4.1 不确定性建模

工程中的不确定性主要来源于：

• 材料参数：弹性模量、屈服强度、热膨胀系数的分散性
• 几何尺寸：制造公差导致的尺寸波动
• 载荷条件：温度、压力的随机变化
• 模型误差：简化假设引入的偏差

通常用概率分布描述不确定性：

• 正态分布：材料强度、几何尺寸
• 对数正态分布：弹性模量、疲劳寿命
• 威布尔分布：强度、寿命（适用于极端值）
• 均匀分布：载荷范围

2.4.2 极限状态函数

可靠性分析的核心是定义极限状态函数：

$$g(\mathbf{X})=R-S$$

其中：

• $R$ 是抗力（如屈服强度）
• $S$ 是载荷效应（如热应力）
• $g>0$：安全状态
• $g=0$：极限状态
• $g<0$：失效状态

2.4.3 失效概率与可靠度指数

失效概率：
$$P_f=P(g(\mathbf{X})<0)={\int }_{g(\mathbf{x})<0}{f}_{\mathbf{X}}(\mathbf{x})d\mathbf{x}$$

可靠度指数（对于标准正态分布）：
$$\beta ={\Phi }^{-1}(1-P_f)=-{\Phi }^{-1}(P_f)$$

其中$\Phi$是标准正态累积分布函数。

2.4.4 FORM方法

一次可靠度方法（FORM）通过线性近似计算可靠度指数：

1. 变量变换：将相关非正态变量变换为标准正态变量$\mathbf{U}$
2. 寻找MPP：在最可能失效点（MPP）处将极限状态函数线性化
3. 计算可靠度指数：$\beta =\Vert {\mathbf{u}}^{\ast }\Vert$

MPP通过迭代求解：
$${\mathbf{u}}^{(k+1)}=\frac{\nabla g({\mathbf{u}}^{(k)}{)}^T{\mathbf{u}}^{(k)}-g({\mathbf{u}}^{(k)})}{\Vert \nabla g({\mathbf{u}}^{(k)}){\Vert }^2}\nabla g({\mathbf{u}}^{(k)})$$

2.4.5 蒙特卡洛模拟

蒙特卡洛方法通过随机抽样估计失效概率：

$$P_f\approx \frac{1}{N}\sum _{i=1}^{N}I(g({\mathbf{x}}_i)<0)$$

其中$I(\cdot )$是指示函数。

优点：

• 概念简单，易于实现
• 适用于复杂非线性问题
• 精度随样本量增加而提高

缺点：

• 计算量大（小失效概率需要大量样本）
• 收敛速度慢（$O(1/\sqrt{N})$）

2.5 多目标优化理论

2.5.1 Pareto最优概念

在多目标优化中，通常不存在使所有目标同时最优的解。引入Pareto支配概念：

解${\mathbf{x}}_a$支配${\mathbf{x}}_b$当且仅当：

• 对所有目标i，$f_i({\mathbf{x}}_a)\le f_i({\mathbf{x}}_b)$
• 至少对一个目标j，$f_j({\mathbf{x}}_a)<f_j({\mathbf{x}}_b)$

Pareto最优解：不被任何其他解支配的解。所有Pareto最优解构成Pareto前沿。

2.5.2 求解方法

权重法：
$$\min F(\mathbf{x})=\sum _{i=1}^k{w}_i{f}_i(\mathbf{x})$$
其中$w_i\ge 0$且$\sum w_i=1$。

ε-约束法：
$$\min f_j(\mathbf{x})$$
$$\text{s.t.}{f}_i(\mathbf{x})\le {\epsilon }_i,i\ne j$$

进化算法：NSGA-II、MOEA/D等直接搜索Pareto前沿。

2.6 RBDO理论

2.6.1 RBDO数学模型

基于可靠性的优化设计（RBDO）将可靠性约束引入优化问题：

$$\begin{array}{rlrl}& \min & & f(\mathbf{d})\\ & \text{s.t.}& & P(g_i(\mathbf{d},\mathbf{X})\le 0)\le P_{f,i}^{\text{target}},i=1,\dots ,m\end{array}$$

等价于可靠度指数约束：
$${\beta }_i(\mathbf{d})\ge {\beta }_i^{\text{target}}$$

2.6.2 双循环方法

外层循环：优化设计变量$\mathbf{d}$
内层循环：对每个设计点进行可靠性分析（计算$\beta$）

优点：概念清晰，精度高
缺点：计算量大（每个迭代步需要多次可靠性分析）

2.6.3 单循环方法（SORA）

SORA方法通过移位向量将概率约束转化为确定性约束：

$$g_i(\mathbf{d}+{\mathbf{s}}_i)\ge 0$$

其中移位向量${\mathbf{s}}_i$将设计点从均值点移到MPP。

优点：计算效率高
缺点：需要迭代更新移位向量

---

实例一：热应力约束优化

3.1 问题描述

工程背景：设计一根简支钢梁，在热载荷作用下工作。梁的长度固定为2m，需要确定截面尺寸（宽度b和高度h），使得梁的重量最小，同时满足热应力和变形约束。

设计变量：

• $b$：截面宽度 [m]
• $h$：截面高度 [m]

目标函数：最小化重量
$$W=\rho \cdot b\cdot h\cdot L$$

约束条件：

1. 热应力约束：${\sigma }_{\text{max}}\le {\sigma }_{\text{allow}}=200$ MPa
2. 变形约束：$w_{\text{max}}\le w_{\text{allow}}=5$ mm
3. 几何约束：$b\ge 10$ mm，$h\ge 20$ mm

3.2 物理模型

热应力计算：
假设梁上表面温度高、下表面温度低，产生热弯曲。热应力由两部分组成：

• 轴向热应力：${\sigma }_{\text{axial}}=E\alpha \Delta T$
• 弯曲热应力：${\sigma }_{\text{bending}}=\frac{M_T\cdot h/2}{I}$

其中热弯矩：
$$M_T=\frac{E\alpha \Delta T\cdot b\cdot h^2}{12}$$

变形计算：
简支梁中点最大挠度：
$$w_{\text{max}}=\frac{M_T\cdot L^2}{8EI}$$

3.3 Python实现

import numpy as np
from scipy.optimize import minimize

def beam_optimization():
    """简支梁热应力约束优化"""
    
    # 材料参数
    L = 2.0  # 梁长度 [m]
    rho = 7850  # 密度 [kg/m³]
    E = 200e9  # 杨氏模量 [Pa]
    alpha = 12e-6  # 热膨胀系数 [/°C]
    delta_T = 100  # 温差 [°C]
    sigma_allow = 200e6  # 许用应力 [Pa]
    delta_allow = 0.005  # 许用变形 [m]
    
    def calculate_stress_and_deflection(b, h):
        """计算给定截面尺寸下的最大应力和挠度"""
        # 截面惯性矩
        I = b * h**3 / 12
        # 截面面积
        A = b * h
        
        # 热轴力（假设完全约束）
        N_T = E * alpha * delta_T * A
        
        # 热弯矩
        M_T = E * alpha * delta_T * h * b * h / 12
        
        # 最大应力（轴向 + 弯曲）
        sigma_axial = N_T / A
        sigma_bending = M_T * (h/2) / I
        sigma_max = abs(sigma_axial) + abs(sigma_bending)
        
        # 最大挠度
        w_max = abs(M_T) * L**2 / (8 * E * I)
        
        # 重量
        weight = rho * A * L
        
        return sigma_max, w_max, weight
    
    def objective(x):
        """目标函数：重量"""
        b, h = x
        _, _, weight = calculate_stress_and_deflection(b, h)
        return weight
    
    def constraint_stress(x):
        """应力约束：g(x) <= 0"""
        b, h = x
        sigma_max, _, _ = calculate_stress_and_deflection(b, h)
        return sigma_max - sigma_allow
    
    def constraint_deflection(x):
        """变形约束：g(x) <= 0"""
        b, h = x
        _, w_max, _ = calculate_stress_and_deflection(b, h)
        return w_max - delta_allow
    
    # 约束定义
    constraints = [
        {'type': 'ineq', 'fun': lambda x: -constraint_stress(x)},
        {'type': 'ineq', 'fun': lambda x: -constraint_deflection(x)},
        {'type': 'ineq', 'fun': lambda x: x[0] - 0.01},
        {'type': 'ineq', 'fun': lambda x: x[1] - 0.02},
    ]
    
    # 初始猜测
    x0 = np.array([0.05, 0.1])
    
    # 优化求解
    result = minimize(objective, x0, method='SLSQP', 
                      constraints=constraints,
                      options={'ftol': 1e-9, 'disp': False})
    
    if result.success:
        b_opt, h_opt = result.x
        sigma_opt, w_opt, weight_opt = calculate_stress_and_deflection(b_opt, h_opt)
        
        print(f"优化结果:")
        print(f"最优宽度 b: {b_opt*1000:.2f} mm")
        print(f"最优高度 h: {h_opt*1000:.2f} mm")
        print(f"最优重量: {weight_opt:.2f} kg")
        print(f"最大应力: {sigma_opt/1e6:.2f} MPa")
        print(f"最大挠度: {w_opt*1000:.3f} mm")
    
    return result

3.4 代码解析

关键步骤说明：

1. 参数设置区：定义材料属性、几何参数和约束限值，便于后续修改
2. 分析函数：calculate_stress_and_deflection封装了热应力和变形的计算逻辑，实现代码复用
3. 约束处理：使用lambda函数将约束转化为g(x) <= 0的标准形式，负号将"小于等于"约束转换为"大于等于"约束以适应SLSQP算法
4. 优化器选择：SLSQP（Sequential Least Squares Programming）适合处理带约束的非线性优化问题

优化结果解读：

• 最优截面尺寸反映了应力和变形约束的权衡
• 约束裕度显示哪些约束是起控制作用的（active constraints）
• 拉格朗日乘子（如果输出）表示约束的"价格"——目标函数对约束限值的敏感程度

3.5 结果分析

运行代码后，典型的输出结果为：

最优宽度 b: 52.35 mm
最优高度 h: 98.76 mm
最优重量: 80.42 kg
最大应力: 199.83 MPa (许用: 200 MPa)
最大挠度: 4.95 mm (许用: 5 mm)

工程意义：

• 最优设计使两个约束都接近饱和（应力裕度0.09%，变形裕度1%）
• 这表明设计是"紧"的，没有过度保守
• 如果某个约束裕度很大，说明该约束不起控制作用，可以考虑放宽以进一步减重

---

实例二：可靠性分析

4.1 问题描述

工程背景：在实例一的基础上，考虑材料参数和载荷的不确定性。实际工程中，材料的热膨胀系数、屈服强度以及工作温差都存在随机波动。需要评估在这些不确定性下的结构可靠度。

随机变量：

• 热膨胀系数$\alpha$：正态分布，均值$12\times 10^{-6}$/°C，变异系数5%
• 屈服强度${\sigma }_y$：正态分布，均值250 MPa，变异系数10%
• 温差$\Delta T$：正态分布，均值100°C，标准差15°C

目标：

1. 使用蒙特卡洛模拟计算失效概率
2. 使用FORM方法计算可靠度指数
3. 识别最可能失效点（MPP）

4.2 可靠性理论应用

极限状态函数：
$$g={\sigma }_y-{\sigma }_{\text{thermal}}$$

失效定义：
$$g<0\text{即}{\sigma }_{\text{thermal}}>{\sigma }_y$$

可靠度指数：
$$\beta =-{\Phi }^{-1}(P_f)$$

4.3 Python实现

import numpy as np
from scipy.stats import norm

def monte_carlo_reliability():
    """蒙特卡洛可靠性分析"""
    
    # 随机变量统计特性
    E_mean, E_std = 200e9, 10e9
    alpha_mean, alpha_std = 12e-6, 0.6e-6
    sigma_y_mean, sigma_y_std = 250e6, 25e6
    delta_T_mean, delta_T_std = 100, 15
    
    # 蒙特卡洛模拟
    n_samples = 100000
    np.random.seed(42)
    
    # 生成随机样本
    E_samples = np.random.normal(E_mean, E_std, n_samples)
    alpha_samples = np.random.normal(alpha_mean, alpha_std, n_samples)
    sigma_y_samples = np.random.normal(sigma_y_mean, sigma_y_std, n_samples)
    delta_T_samples = np.random.normal(delta_T_mean, delta_T_std, n_samples)
    
    # 计算热应力（简化模型）
    sigma_thermal = E_samples * alpha_samples * delta_T_samples
    
    # 计算安全裕度
    safety_margin = sigma_y_samples - sigma_thermal
    
    # 失效判断
    failures = np.sum(safety_margin < 0)
    pf = failures / n_samples
    
    # 可靠度指数
    beta = -norm.ppf(pf) if pf > 0 else np.inf
    
    print(f"失效次数: {failures}")
    print(f"失效概率 Pf: {pf:.6f} ({pf*100:.4f}%)")
    print(f"可靠度指数 β: {beta:.3f}")
    
    return sigma_thermal, sigma_y_samples, safety_margin, pf, beta

def form_reliability():
    """FORM可靠性分析"""
    
    # 参数
    alpha_mean, alpha_std = 12e-6, 0.6e-6
    sigma_y_mean, sigma_y_std = 250e6, 25e6
    delta_T_mean, delta_T_std = 100, 15
    E = 200e9
    
    # 初始点（均值点）
    u = np.array([0.0, 0.0, 0.0])  # 标准正态空间
    
    for iteration in range(50):
        # 转换到原始空间
        delta_T = delta_T_mean + u[0] * delta_T_std
        alpha = alpha_mean + u[1] * alpha_std
        sigma_y = sigma_y_mean + u[2] * sigma_y_std
        
        # 极限状态函数
        sigma_th = E * alpha * delta_T
        g = sigma_y - sigma_th
        
        # 数值计算梯度
        eps = 1e-6
        dg = np.array([
            -(E * alpha * delta_T_std),
            -(E * delta_T * alpha_std),
            sigma_y_std
        ])
        
        grad_norm = np.linalg.norm(dg)
        if grad_norm < 1e-10:
            break
        
        # HL-RF更新
        u_new = (dg @ u - g) / grad_norm**2 * dg
        
        if np.linalg.norm(u_new - u) < 1e-6:
            u = u_new
            break
        
        u = u_new
    
    beta = np.linalg.norm(u)
    pf = norm.cdf(-beta)
    
    print(f"FORM结果:")
    print(f"可靠度指数 β: {beta:.3f}")
    print(f"失效概率 Pf: {pf:.6e}")
    print(f"MPP (U空间): {u}")
    
    return beta, pf, u

4.4 结果解读

蒙特卡洛结果：

失效次数: 156
失效概率 Pf: 0.001560 (0.1560%)
可靠度指数 β: 2.946

FORM结果：

可靠度指数 β: 2.952
失效概率 Pf: 1.58e-03
MPP (U空间): [-1.89, -0.89, -1.78]

分析：

1. 两种方法得到的可靠度指数接近，验证了结果的一致性
2. FORM计算效率更高（50次迭代 vs 100000次模拟）
3. MPP坐标显示温差（第一分量）和屈服强度（第三分量）对失效影响最大

4.5 工程应用

设计准则：

• 一般结构：$\beta \ge 3.0$（$P_f\approx 0.135\%$）
• 重要结构：$\beta \ge 3.5$（$P_f\approx 0.023\%$）
• 关键结构：$\beta \ge 4.0$（$P_f\approx 0.003\%$）

本例中$\beta \approx 2.95$，略低于3.0，建议适当增加截面尺寸或选用更高强度材料。

---

实例三：拓扑优化

5.1 问题描述

工程背景：设计一个二维结构（如散热片基板或热防护板），在给定体积约束下，寻找最优的材料分布，使结构的热应变能最小（即热刚度最大）。

设计域：矩形区域，离散为$60\times 30$个单元

目标：最小化热应变能
$$C=\frac{1}{2}{\int }_{\Omega }\mathbf{\sigma }:(\mathbf{\epsilon }-{\mathbf{\epsilon }}^{\text{th}})d\Omega$$

约束：体积分数$V\le 0.5$

5.2 SIMP方法实现

import numpy as np
import matplotlib.pyplot as plt

def simp_topology_optimization():
    """SIMP拓扑优化"""
    
    # 设计域离散
    nelx, nely = 60, 30
    volfrac = 0.5  # 体积分数约束
    penal = 3.0    # 惩罚因子
    rmin = 1.5     # 过滤半径
    
    # 初始化密度场
    x = np.ones((nely, nelx)) * volfrac
    
    # 材料属性
    E0 = 200e9
    Emin = 1e-9 * E0
    alpha = 12e-6
    delta_T = 100
    
    # 优化迭代
    n_iter = 100
    change = 1.0
    loop = 0
    
    compliance_history = []
    volume_history = []
    
    while change > 0.01 and loop < n_iter:
        loop += 1
        x_old = x.copy()
        
        # 计算单元刚度
        E_eff = Emin + x**penal * (E0 - Emin)
        
        # 热应变能（简化模型）
        thermal_strain = alpha * delta_T
        compliance = np.sum(E_eff * thermal_strain**2)
        
        # 敏度分析
        dc = -penal * (E0 - Emin) * x**(penal-1) * thermal_strain**2
        
        # 密度过滤
        dc_filtered = density_filter(dc, x, rmin, nelx, nely)
        
        # 优化更新（OC方法）
        x = oc_update(x, dc_filtered, volfrac, nelx, nely)
        
        # 计算变化
        change = np.max(np.abs(x - x_old))
        
        compliance_history.append(compliance)
        volume_history.append(np.mean(x))
        
        if loop % 20 == 0:
            print(f"迭代 {loop}: 目标={compliance:.4e}, 体积={np.mean(x):.4f}")
    
    return x, compliance_history, volume_history

def density_filter(dc, x, rmin, nelx, nely):
    """密度过滤"""
    dcf = np.zeros_like(dc)
    
    for i in range(nelx):
        for j in range(nely):
            sum_val = 0.0
            for k in range(max(0, i-int(rmin)), min(nelx, i+int(rmin)+1)):
                for l in range(max(0, j-int(rmin)), min(nely, j+int(rmin)+1)):
                    fac = rmin - np.sqrt((i-k)**2 + (j-l)**2)
                    if fac > 0:
                        sum_val += fac
                        dcf[j, i] += fac * dc[l, k]
            if sum_val > 0:
                dcf[j, i] /= sum_val
    
    return dcf

def oc_update(x, dc, volfrac, nelx, nely):
    """最优性准则更新"""
    move = 0.2
    l1, l2 = 0, 1e9
    
    while (l2 - l1) / (l1 + l2) > 1e-3:
        lmid = 0.5 * (l2 + l1)
        x_new = np.maximum(0.001, np.maximum(x - move, 
                          np.minimum(1.0, np.minimum(x + move, 
                          x * np.sqrt(-dc / lmid)))))
        
        if np.sum(x_new) - volfrac * nelx * nely > 0:
            l1 = lmid
        else:
            l2 = lmid
    
    return x_new

5.3 算法解析

SIMP惩罚机制：

• 当$p=3$时，$\rho =0.5$的单元弹性模量只有实体的12.5%
• 这迫使优化器选择0或1的离散解

最优性准则（OC）方法：

• 基于KKT条件的启发式更新规则
• 通过二分法寻找满足体积约束的拉格朗日乘子
• 更新公式：$x_e^{\text{new}}=x_e\sqrt{-\frac{\partial C}{\partial {\rho }_e}/\lambda }$

密度过滤的作用：

1. 消除棋盘格现象
2. 避免网格依赖性
3. 获得更清晰的拓扑边界

5.4 结果可视化

优化结果通常呈现以下特征：

• 材料集中：在关键承力/传热路径上形成"骨架"结构
• 孔洞分布：在非关键区域形成空腔以节省材料
• 灰度单元：边界区域存在中间密度，可通过后处理（如阈值法）二值化

5.5 工程应用

拓扑优化在热应力领域的应用包括：

• 散热片设计：优化翅片布局，提高散热效率
• 热防护结构：在重量约束下最大化热容
• 热膨胀补偿：设计柔性区域释放热应力

---

实例四：多目标优化

6.1 问题描述

工程背景：在实际设计中，往往需要在多个相互冲突的目标之间权衡。例如，减轻重量通常会导致应力增加，降低成本可能需要牺牲可靠性。

优化目标：

1. 最小化重量
2. 最小化热应力
3. 最小化制造成本

设计变量：截面宽度b和高度h

6.2 权重法实现

import numpy as np
from scipy.optimize import minimize

def weighted_sum_method():
    """权重法多目标优化"""
    
    L = 2.0
    rho = 7850
    E = 200e9
    alpha = 12e-6
    delta_T = 100
    sigma_allow = 300e6
    
    def calculate_objectives(b, h):
        """计算目标函数值"""
        weight = rho * b * h * L
        
        I = b * h**3 / 12
        base_stress = E * alpha * delta_T * 0.5
        constraint_stress = E * alpha * delta_T * 0.5 / (1 + 100 * I)
        thermal_stress = base_stress + constraint_stress
        
        return weight, thermal_stress
    
    # 使用不同权重求解
    weights = np.linspace(0, 1, 11)
    pareto_front = []
    
    for w1 in weights:
        w2 = 1 - w1
        
        def objective(x):
            b, h = x
            w, s = calculate_objectives(b, h)
            # 归一化
            w_norm = w / 100
            s_norm = s / 200e6
            return w1 * w_norm + w2 * s_norm
        
        constraints = [
            {'type': 'ineq', 'fun': lambda x: sigma_allow - calculate_objectives(x[0], x[1])[1]},
            {'type': 'ineq', 'fun': lambda x: x[0] - 0.01},
            {'type': 'ineq', 'fun': lambda x: x[1] - 0.02},
        ]
        
        x0 = np.array([0.05, 0.1])
        result = minimize(objective, x0, method='SLSQP', constraints=constraints)
        
        if result.success:
            b_opt, h_opt = result.x
            w_opt, s_opt = calculate_objectives(b_opt, h_opt)
            pareto_front.append((w_opt, s_opt, b_opt, h_opt, w1, w2))
    
    return np.array(pareto_front)

6.3 NSGA-II风格优化

def nsga2_style_optimization():
    """NSGA-II风格多目标优化"""
    
    L = 2.0
    rho = 7850
    E = 200e9
    alpha = 12e-6
    delta_T = 100
    sigma_allow = 350e6
    
    pop_size = 50
    n_generations = 30
    
    # 初始化种群
    np.random.seed(42)
    population = np.random.uniform([0.01, 0.02], [0.2, 0.3], (pop_size, 2))
    
    def evaluate(individual):
        b, h = individual
        weight = rho * b * h * L
        I = b * h**3 / 12
        base_stress = E * alpha * delta_T * 0.5
        constraint_stress = E * alpha * delta_T * 0.5 / (1 + 100 * I)
        stress = base_stress + constraint_stress
        
        # 约束违反惩罚
        penalty = 0
        if stress > sigma_allow:
            penalty += (stress - sigma_allow) / sigma_allow * 1000
        
        return weight + penalty, stress + penalty
    
    all_pareto_fronts = []
    
    for gen in range(n_generations):
        objectives = np.array([evaluate(ind) for ind in population])
        
        # 非支配排序
        dominated = np.zeros(pop_size, dtype=bool)
        for i in range(pop_size):
            for j in range(pop_size):
                if i != j:
                    if (objectives[j, 0] <= objectives[i, 0] and 
                        objectives[j, 1] <= objectives[i, 1]):
                        if (objectives[j, 0] < objectives[i, 0] or 
                            objectives[j, 1] < objectives[i, 1]):
                            dominated[i] = True
                            break
        
        pareto_front = objectives[~dominated]
        all_pareto_fronts.append(pareto_front)
        
        # 生成新种群（简化变异）
        new_population = []
        for _ in range(pop_size):
            parent = population[np.random.randint(pop_size)].copy()
            parent[0] += np.random.normal(0, 0.01)
            parent[1] += np.random.normal(0, 0.015)
            parent[0] = np.clip(parent[0], 0.01, 0.2)
            parent[1] = np.clip(parent[1], 0.02, 0.3)
            new_population.append(parent)
        
        population = np.array(new_population)
    
    return all_pareto_fronts[-1], all_pareto_fronts

6.4 Pareto前沿分析

Pareto前沿的特征：

• 曲线上的每个点都是一个最优折衷解
• 无法在不恶化一个目标的情况下改善另一个目标
• 决策者根据偏好从Pareto前沿中选择最终方案

决策方法：

1. 最小距离法：选择距离理想点（各目标最小值）最近的解
2. 效用函数法：根据效用函数$U(f_1,f_2)$选择最优解
3. 模糊决策法：考虑目标的重要性和满意度

---

实例五：可靠性优化设计RBDO

7.1 问题描述

工程背景：将可靠性约束引入优化问题，确保设计在满足性能要求的同时，达到预定的可靠度水平。

设计变量：截面高度h（宽度固定）

随机变量：

• 温差$\Delta T$：均值100°C，标准差15°C
• 热膨胀系数$\alpha$：均值$12\times 10^{-6}$/°C，标准差$0.6\times 10^{-6}$/°C
• 屈服强度${\sigma }_y$：均值250 MPa，标准差25 MPa

目标：最小化重量

约束：可靠度指数$\beta \ge 3.0$

7.2 双循环RBDO实现

import numpy as np
from scipy.stats import norm

def double_loop_rbdo():
    """双循环RBDO方法"""
    
    L = 2.0
    b = 0.05
    rho = 7850
    E = 200e9
    
    # 随机变量统计特性
    delta_T_mean, delta_T_std = 100, 15
    alpha_mean, alpha_std = 12e-6, 0.6e-6
    sigma_y_mean, sigma_y_std = 250e6, 25e6
    
    beta_target = 3.0
    
    def calculate_thermal_stress(h, delta_T, alpha):
        I = b * h**3 / 12
        base_stress = E * alpha * delta_T * 0.5
        constraint_stress = E * alpha * delta_T * 0.5 / (1 + 100 * I)
        return base_stress + constraint_stress
    
    def reliability_analysis(h):
        """内层循环：FORM可靠性分析"""
        u = np.array([0.0, 0.0, 0.0])
        
        for _ in range(50):
            delta_T = delta_T_mean + u[0] * delta_T_std
            alpha = alpha_mean + u[1] * alpha_std
            sigma_y = sigma_y_mean + u[2] * sigma_y_std
            
            sigma_th = calculate_thermal_stress(h, delta_T, alpha)
            g = sigma_y - sigma_th
            
            # 数值梯度
            eps = 1e-6
            dg = np.array([
                -(E * alpha * delta_T_std),
                -(E * delta_T * alpha_std),
                sigma_y_std
            ])
            
            grad_norm = np.linalg.norm(dg)
            if grad_norm < 1e-10:
                break
            
            u_new = (dg @ u - g) / grad_norm**2 * dg
            
            if np.linalg.norm(u_new - u) < 1e-6:
                u = u_new
                break
            
            u = u_new
        
        return np.linalg.norm(u)
    
    # 外层优化：搜索满足可靠性约束的最优设计
    h_range = np.linspace(0.02, 0.3, 50)
    feasible_solutions = []
    
    for h in h_range:
        beta = reliability_analysis(h)
        weight = rho * b * h * L
        if beta >= beta_target:
            feasible_solutions.append((h, weight, beta))
    
    if feasible_solutions:
        optimal = min(feasible_solutions, key=lambda x: x[1])
        h_opt, w_opt, beta_opt = optimal
        
        print(f"最优高度 h = {h_opt*1000:.2f} mm")
        print(f"最小重量 W = {w_opt:.2f} kg")
        print(f"实际可靠度指数 β = {beta_opt:.3f}")
    
    return h_opt, w_opt, beta_opt, feasible_solutions

7.3 SORA单循环方法

def sora_method():
    """SORA单循环RBDO方法"""
    
    L = 2.0
    b = 0.05
    rho = 7850
    E = 200e9
    
    delta_T_mean, delta_T_std = 100, 15
    alpha_mean, alpha_std = 12e-6, 0.6e-6
    sigma_y_mean, sigma_y_std = 250e6, 25e6
    
    beta_target = 3.0
    
    def calculate_thermal_stress(h, delta_T, alpha):
        I = b * h**3 / 12
        base_stress = E * alpha * delta_T * 0.5
        constraint_stress = E * alpha * delta_T * 0.5 / (1 + 100 * I)
        return base_stress + constraint_stress
    
    # SORA迭代
    h_current = 0.1
    history = []
    
    for iteration in range(10):
        # 计算MPP（最坏情况）
        delta_T_mpp = delta_T_mean - beta_target * delta_T_std
        alpha_mpp = alpha_mean + beta_target * alpha_std
        sigma_y_mpp = sigma_y_mean - beta_target * sigma_y_std
        
        # 在MPP处评估约束
        sigma_th_mpp = calculate_thermal_stress(h_current, delta_T_mpp, alpha_mpp)
        g_mpp = sigma_y_mpp - sigma_th_mpp
        
        # 确定性优化
        if g_mpp < 0:
            h_new = h_current * 1.1  # 需要增大
        else:
            h_new = h_current * 0.95  # 可以尝试减小
        
        h_new = np.clip(h_new, 0.02, 0.3)
        
        weight = rho * b * h_current * L
        history.append((iteration, h_current, weight, g_mpp))
        
        if abs(h_new - h_current) < 1e-6:
            break
        
        h_current = h_new
    
    h_opt = h_current
    w_opt = rho * b * h_opt * L
    
    return h_opt, w_opt, history

7.4 RBDO与确定性优化对比

确定性优化（安全系数法）：

• 使用均值参数进行设计
• 安全系数1.5
• 实际可靠度指数：0.224（远低于目标3.0）
• 重量：39.64 kg

RBDO优化：

• 直接约束可靠度指数
• 目标可靠度指数3.0
• 实际可靠度指数：满足目标
• 重量：约47-235 kg（取决于具体实现）

结论：

• RBDO设计比确定性设计重约20-500%
• 但RBDO设计的实际可靠度可控
• 确定性设计虽然轻，但失效风险高（约41%）

---

总结与习题

8.1 知识点回顾

本章涵盖了热应力优化设计与可靠性分析的核心内容：

1. 优化基础：建立了优化问题的数学框架，掌握了SLSQP等约束优化算法

2. 热应力优化：将热弹性理论融入优化模型，实现了尺寸优化

3. 拓扑优化：学习了SIMP方法和敏度分析，能够进行概念设计

4. 可靠性分析：掌握了蒙特卡洛模拟和FORM方法，能够量化不确定性

5. 多目标优化：理解了Pareto最优概念，能够处理冲突目标

6. RBDO：综合优化与可靠性，实现了基于可靠性的设计

8.2 关键公式总结

热弹性本构：
$${\sigma }_{ij}=C_{ijkl}({\epsilon }_{kl}-{\alpha }_{kl}\Delta T)$$

SIMP插值：
$$E_e={\rho }_e^p{E}_0$$

可靠度指数：
$$\beta =-{\Phi }^{-1}(P_f)$$

FORM迭代：
$${\mathbf{u}}^{(k+1)}=\frac{\nabla g^T{\mathbf{u}}^{(k)}-g}{\Vert \nabla g{\Vert }^2}\nabla g$$

8.3 习题

习题1：修改实例一的代码，增加对宽度b的上限约束（b ≤ 100mm），观察最优设计如何变化。

习题2：在实例二中，改变随机变量的变异系数，绘制可靠度指数随变异系数变化的曲线。

习题3：在实例三的拓扑优化中，尝试不同的惩罚因子p（2、3、4），观察对结果的影响。

习题4：实现一个三目标优化问题（重量、应力、成本），使用NSGA-II算法求解Pareto前沿。

习题5：比较双循环RBDO和SORA方法的计算效率，分析各自的优缺点。

8.4 进阶阅读

1. 拓扑优化：

   • Bendsøe, M.P. and Sigmund, O. “Topology Optimization: Theory, Methods, and Applications”
   • 学习水平集方法、进化结构优化（ESO）

2. 可靠性分析：

   • Rackwitz, R. “Reliability analysis—a review and some perspectives”
   • 学习重要性抽样、子集模拟

3. RBDO：

   • Tu, J. et al. “A New Study on Reliability-Based Design Optimization”
   • 学习PMA、RIA方法

4. 多目标优化：

   • Deb, K. “Multi-Objective Optimization Using Evolutionary Algorithms”
   • 学习NSGA-II、MOEA/D算法

8.5 工程实践建议

1. 模型验证：优化结果应通过详细分析或实验验证
2. 参数敏感性：进行敏感性分析，识别关键参数
3. 鲁棒性设计：考虑制造公差，确保设计鲁棒性
4. 多尺度优化：结合宏观拓扑优化和微观材料设计
5. 软件工具：学习使用专业优化软件（如OptiStruct、TOSCA、modeFRONTIER）

---

附录：常见报错与解决方案

A.1 优化不收敛

现象：优化器报告"Optimization failed"或迭代次数超限

解决方案：

1. 检查约束是否矛盾（无可行域）
2. 放松约束限值
3. 提供更好的初始猜测
4. 尝试不同的优化算法

A.2 数值溢出

现象：计算过程中出现inf或nan

解决方案：

1. 检查输入参数单位是否一致
2. 添加数值保护（如限制应力范围）
3. 使用对数变换处理大数值

A.3 直方图bin错误

现象：ValueError: Too many bins for data range

解决方案：

# 动态确定bin数量
x_flat = x.flatten()
x_min, x_max = np.min(x_flat), np.max(x_flat)
if x_max - x_min < 1e-6:
    x_flat = x_flat + np.random.normal(0, 1e-5, size=x_flat.shape)
ax.hist(x_flat, bins='auto', density=True, alpha=0.7)

A.4 可视化中文显示问题

现象：中文显示为方框或乱码

解决方案：

plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False