结构优化设计-主题008-形状优化设计

kkchenkx

118人浏览 · 2026-04-12 08:48:18

kkchenkx · 2026-04-12 08:48:18 发布

主题008：形状优化设计

一、引言

形状优化设计是结构优化领域的重要组成部分，它通过改变结构的边界形状或内部几何轮廓来改善结构性能。与尺寸优化仅改变截面尺寸不同，形状优化直接改变结构的几何形态，能够更充分地发挥材料潜力，实现更优的结构性能。

形状优化在工程实践中有着广泛的应用：飞机机翼的气动外形优化、汽车车身流线型设计、涡轮叶片型线优化、压力容器开口形状设计、桥梁拱轴线优化等。这些应用充分体现了形状优化在提升产品性能、降低能耗、延长使用寿命等方面的重要价值。

本教程将系统介绍形状优化的基本理论、数学模型、参数化方法、灵敏度分析技术以及数值求解策略，并通过完整的Python代码实现，帮助读者深入理解形状优化的核心原理和工程实现方法。
在这里插入图片描述

二、形状优化的基本概念

2.1 什么是形状优化

形状优化（Shape Optimization）是指在给定设计域内，通过调整结构的边界形状或内部几何轮廓，使目标函数达到最优，同时满足各种约束条件的设计方法。

形状优化的核心特征包括：

设计变量是几何形状：形状优化的设计变量直接描述结构的几何边界，可以是边界节点的坐标、控制曲线的参数、或者几何基函数的系数等。
拓扑保持不变：与拓扑优化不同，形状优化不改变结构的连通性，只改变已有边界的形状。
网格随形状变化：在形状优化过程中，有限元网格需要随边界形状的变化而更新，这带来了额外的计算复杂性。
灵敏度分析更复杂：形状灵敏度涉及网格变形对响应的影响，需要特殊的数学处理。

2.2 形状优化的数学表述

形状优化问题的一般数学表述为：

$\begin{aligned} \min_{\mathbf{s}} \quad & f(\mathbf{s}, \mathbf{u}(\mathbf{s})) \\ \text{s.t.} \quad & g_j(\mathbf{s}, \mathbf{u}(\mathbf{s})) \leq 0, \quad j = 1, 2, \ldots, m \\ & h_k(\mathbf{s}, \mathbf{u}(\mathbf{s})) = 0, \quad k = 1, 2, \ldots, p \\ & \mathbf{s}^{\min} \leq \mathbf{s} \leq \mathbf{s}^{\max} \end{aligned}$

其中：

$s\mathbf{s}$ 是形状设计变量向量
$f$ 是目标函数（如重量、柔度、应力等）
$g_j$ 是不等式约束（如应力约束、位移约束等）
$h_k$ 是等式约束
$u(s)\mathbf{u}(\mathbf{s})$ 是依赖于形状的位移场，通过有限元分析获得

2.3 形状优化的分类

根据设计变量的定义方式，形状优化可以分为以下几类：

1. 节点坐标法

直接将边界节点的坐标作为设计变量：

$\mathbf{s} = [x_1, y_1, x_2, y_2, \ldots, x_n, y_n]^T$

优点：直观简单，易于实现。

缺点：设计变量多，容易产生锯齿形边界；难以保证边界的光滑性。

2. 设计元素法

将设计域划分为若干设计元素，每个元素内的形状由少量参数控制。例如，用三次样条曲线描述边界，控制点坐标作为设计变量。

优点：设计变量少，边界光滑。

缺点：形状变化范围受限，可能错过最优解。

3. 基函数法

用一组基函数的线性组合描述边界形状：

$\mathbf{r}(\xi) = \sum_{i=1}^{n} s_i \mathbf{\phi}_i(\xi)$

其中 $ϕi(ξ)\mathbf{\phi}_i(\xi)$ 是基函数， $s_i$ 是设计变量。

常用基函数包括：

多项式基函数
三角函数基函数
样条基函数
径向基函数

4. 自由形状法

基于变分原理和速度场概念，允许边界在法向方向上自由变形。这种方法与拓扑优化结合，形成了水平集方法等先进技术。

2.4 形状优化的挑战

形状优化面临以下主要挑战：

1. 网格变形问题

随着边界形状的变化，内部网格需要相应更新。网格变形算法需要保证：

网格质量不恶化（避免出现畸形单元）
计算效率高
边界节点精确落在设计边界上

常用网格变形方法包括：

代数插值法
偏微分方程法（如弹性变形法、Laplace光滑法）
显式重网格化

2. 灵敏度分析复杂性

形状灵敏度不仅与材料属性、截面尺寸有关，还与网格节点位置的变化相关。这需要引入**物质导数（Material Derivative）**的概念。

3. 多尺度问题

形状变化可能涉及多个尺度：宏观外形、细观结构、微观组织。多尺度形状优化是当前研究的前沿课题。

4. 制造约束

优化得到的复杂形状可能难以制造。需要考虑：

最小特征尺寸约束
拔模角度约束
对称性约束
加工工艺约束

三、形状参数化方法

形状参数化是形状优化的基础，它决定了设计空间的维度和性质。一个好的参数化方法应该满足：

完备性：能够描述足够广泛的形状
紧凑性：设计变量尽可能少
光滑性：生成的边界足够光滑
鲁棒性：参数变化不会导致病态形状

3.1 B样条曲线参数化

B样条（B-Spline）曲线是形状参数化中最常用的工具之一。B样条曲线由控制点和基函数定义：

$\mathbf{C}(u) = \sum_{i=0}^{n} N_{i,p}(u) \mathbf{P}_i$

其中：

$Pi\mathbf{P}_i$ 是控制点
$N_{i,p}(u)$ 是 $p$ 次B样条基函数
$\in [0,1]$ 是参数

B样条基函数由Cox-de Boor递推公式定义：

$N_{i,0}(u) = \begin{cases} 1 & \text{if } u_i \leq u < u_{i+1} \\ 0 & \text{otherwise} \end{cases}$

$N_{i,p}(u) = \frac{u - u_i}{u_{i+p} - u_i} N_{i,p-1}(u) + \frac{u_{i+p+1} - u}{u_{i+p+1} - u_{i+1}} N_{i+1,p-1}(u)$

B样条曲线的优点：

局部控制性：移动一个控制点只影响局部曲线
光滑性： $p$ 次B样条具有 $C^{p-1}$ 连续性
仿射不变性：对控制点进行仿射变换等价于对曲线进行相同变换

3.2 NURBS参数化

非均匀有理B样条（NURBS, Non-Uniform Rational B-Spline）是B样条的推广，可以精确表示圆锥曲线（圆、椭圆、抛物线、双曲线）：

$\mathbf{C}(u) = \frac{\sum_{i=0}^{n} N_{i,p}(u) w_i \mathbf{P}_i}{\sum_{i=0}^{n} N_{i,p}(u) w_i}$

其中 $w_i$ 是权重因子。NURBS在CAD/CAM系统中广泛应用，是形状优化的理想参数化工具。

3.3 径向基函数参数化

径向基函数（RBF, Radial Basis Function）提供了一种灵活的参数化方式：

$\mathbf{r}(\mathbf{x}) = \mathbf{x} + \sum_{i=1}^{n} \alpha_i \phi(||\mathbf{x} - \mathbf{c}_i||)$

常用RBF核函数：

高斯函数： $ϕ(r)=e−(εr)2\phi(r) = e^{-(\varepsilon r)^2}$
多二次函数： $ϕ(r)=1+(εr)2\phi(r) = \sqrt{1 + (\varepsilon r)^2}$
薄板样条： $ϕ(r)=r2ln⁡(r)\phi(r) = r^2 \ln(r)$

RBF参数化的优点：

可以处理任意拓扑的变形
提供全局光滑的变形场
适合与网格变形结合

四、形状灵敏度分析

形状灵敏度分析是形状优化的核心，它计算目标函数和约束对形状设计变量的导数。

4.1 物质导数概念

在形状优化中，场量（如位移、应力）同时依赖于空间位置和边界形状。物质导数描述了场量随形状变化的速率：

设 $ϕ(x,s)\phi(\mathbf{x}, \mathbf{s})$ 是依赖于形状 $s\mathbf{s}$ 的场量，形状变化引起的位置变化为 $V\mathbf{V}$ （速度场），则物质导数为：

$\frac{D\phi}{D\mathbf{s}} = \frac{\partial \phi}{\partial \mathbf{s}} + \nabla \phi \cdot \mathbf{V}$

其中：

$∂ϕ∂s\frac{\partial \phi}{\partial \mathbf{s}}$ 是局部导数（固定空间点）
$∇ϕ⋅V\nabla \phi \cdot \mathbf{V}$ 是对流项（由于网格移动）

4.2 直接法形状灵敏度

直接法通过求解灵敏度方程获得位移灵敏度。

对于有限元方程 $K(s)u=F(s)\mathbf{K}(\mathbf{s})\mathbf{u} = \mathbf{F}(\mathbf{s})$ ，对设计变量 $s_i$ 求导：

$\frac{\partial \mathbf{K}}{\partial s_i} \mathbf{u} + \mathbf{K} \frac{\partial \mathbf{u}}{\partial s_i} = \frac{\partial \mathbf{F}}{\partial s_i}$

整理得灵敏度方程：

$\mathbf{K} \frac{\partial \mathbf{u}}{\partial s_i} = \frac{\partial \mathbf{F}}{\partial s_i} - \frac{\partial \mathbf{K}}{\partial s_i} \mathbf{u}$

右端项称为伪载荷。求解该方程即可得到位移灵敏度 $∂u∂si\frac{\partial \mathbf{u}}{\partial s_i}$ 。

刚度矩阵灵敏度 $∂K∂si\frac{\partial \mathbf{K}}{\partial s_i}$ 的计算涉及：

材料属性的变化
单元几何的变化（节点坐标变化）
网格变形的影响

4.3 伴随法形状灵敏度

当设计变量多、响应少时，伴随法更高效。

对于响应 $\mathbf{c}^T \mathbf{u}$ ，其灵敏度为：

$\frac{dg}{ds_i} = \mathbf{c}^T \frac{\partial \mathbf{u}}{\partial s_i} = \mathbf{c}^T \mathbf{K}^{-1} \left( \frac{\partial \mathbf{F}}{\partial s_i} - \frac{\partial \mathbf{K}}{\partial s_i} \mathbf{u} \right)$

定义伴随变量 $λ\boldsymbol{\lambda}$ 满足：

$\mathbf{K} \boldsymbol{\lambda} = \mathbf{c}$

则灵敏度可表示为：

$\frac{dg}{ds_i} = \boldsymbol{\lambda}^T \left( \frac{\partial \mathbf{F}}{\partial s_i} - \frac{\partial \mathbf{K}}{\partial s_i} \mathbf{u} \right)$

伴随法只需对每个响应求解一次伴随方程，与 design variables 数量无关。

4.4 半解析法

对于复杂问题，解析求导困难，可以采用半解析法：

$\frac{\partial \mathbf{K}}{\partial s_i} \approx \frac{\mathbf{K}(s_i + \Delta s_i) - \mathbf{K}(s_i)}{\Delta s_i}$

半解析法实现简单，但需要注意：

差分步长的选择（过大截断误差大，过小舍入误差大）
计算成本较高（需要多次组装刚度矩阵）

五、网格变形技术

在形状优化中，边界形状变化后，内部网格需要相应更新。网格变形算法是形状优化的关键技术之一。

5.1 弹性变形法

将网格视为弹性体，边界位移作为强制位移，求解弹性问题得到内部节点位移。

弹性问题的控制方程：

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

其中应力-应变关系为：

$\boldsymbol{\sigma} = \lambda \, \text{tr}(\boldsymbol{\varepsilon}) \mathbf{I} + 2\mu \boldsymbol{\varepsilon}$

参数 $λ\lambda$ 和 $μ\mu$ 可以调整以控制变形特性。通常设置较小的弹性模量，使网格容易变形。

5.2 Laplace光滑法

Laplace光滑是一种简单有效的网格光滑方法。对于内部节点 $i$ ，其新位置为邻接节点的平均：

$\mathbf{x}_i^{\text{new}} = \frac{1}{n_i} \sum_{j \in \mathcal{N}(i)} \mathbf{x}_j$

其中 $N(i)\mathcal{N}(i)$ 是节点 $i$ 的邻接节点集合， $n_i$ 是邻接节点数。

Laplace光滑迭代进行，直到收敛。这种方法简单高效，但对于大变形可能产生网格交叉。

5.3 径向基函数网格变形

RBF不仅可以用于形状参数化，也可以用于网格变形：

$\mathbf{d}(\mathbf{x}) = \sum_{i=1}^{n} \boldsymbol{\alpha}_i \phi(||\mathbf{x} - \mathbf{c}_i||)$

其中 $ci\mathbf{c}_i$ 是控制点， $αi\boldsymbol{\alpha}_i$ 是待定系数，由边界位移条件确定。

RBF变形的优点：

全局光滑
可以处理大变形
计算效率高（一旦确定系数，变形可快速计算）

六、形状优化算法

6.1 基于梯度的优化算法

形状优化通常采用基于梯度的算法，如：

序列线性规划（SLP）

在当前设计点线性化目标函数和约束，求解线性规划子问题：

$\begin{aligned} \min_{\Delta \mathbf{s}} \quad & \nabla f^T \Delta \mathbf{s} \\ \text{s.t.} \quad & g_j + \nabla g_j^T \Delta \mathbf{s} \leq 0 \\ & |\Delta s_i| \leq \Delta s_i^{\max} \end{aligned}$

序列二次规划（SQP）

在当前设计点构建二次规划子问题：

$\begin{aligned} \min_{\Delta \mathbf{s}} \quad & \nabla f^T \Delta \mathbf{s} + \frac{1}{2} \Delta \mathbf{s}^T \mathbf{H} \Delta \mathbf{s} \\ \text{s.t.} \quad & g_j + \nabla g_j^T \Delta \mathbf{s} \leq 0 \\ & h_k + \nabla h_k^T \Delta \mathbf{s} = 0 \end{aligned}$

其中 $H\mathbf{H}$ 是Hessian矩阵或其近似。

移动渐近线法（MMA）

MMA是结构优化中广泛使用的算法，通过引入移动渐近线将原问题转化为一系列凸近似子问题：

$\begin{aligned} \min_{\mathbf{s}} \quad & \sum_{i} \left( \frac{p_i^{(k)}}{U_i^{(k)} - s_i} + \frac{q_i^{(k)}}{s_i - L_i^{(k)}} \right) + r_0^{(k)} \\ \text{s.t.} \quad & \sum_{i} \left( \frac{p_{ji}^{(k)}}{U_i^{(k)} - s_i} + \frac{q_{ji}^{(k)}}{s_i - L_i^{(k)}} \right) + r_j^{(k)} \leq 0 \end{aligned}$

其中 $L_i^{(k)}$ 和 $U_i^{(k)}$ 是第 $k$ 次迭代的下渐近线和上渐近线。

6.2 无梯度优化算法

对于非光滑或难以求导的问题，可以采用无梯度算法：

遗传算法（GA）
粒子群优化（PSO）
模拟退火（SA）
模式搜索（Pattern Search）

无梯度算法的优点是适用范围广，缺点是计算成本高，收敛慢。

七、Python实现：二维孔洞形状优化

下面通过完整的Python代码实现一个二维孔洞形状优化问题。优化目标是：在给定载荷条件下，通过改变孔洞形状，使孔边最大应力最小化（即降低应力集中）。

7.1 问题描述

考虑一个带中心圆孔的无限大板受单向拉伸。圆孔引起应力集中，最大应力为 $3σ03\sigma_0$ （ $σ0\sigma_0$ 是远场应力）。

优化目标：将圆孔变为椭圆或其他形状，使应力集中系数最小化。

这是一个经典的形状优化问题，理论最优解是椭圆孔，长短轴比为3:1，此时应力集中系数降为2。

7.2 完整Python代码

"""
形状优化设计：二维孔洞形状优化
目标：通过改变孔洞形状，最小化应力集中
方法：B样条参数化 + 有限元分析 + 梯度优化
"""

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Polygon, Circle, FancyBboxPatch
from matplotlib.collections import PatchCollection
import matplotlib.animation as animation
from scipy.optimize import minimize, differential_evolution
from scipy.interpolate import splev, splprep, CubicSpline
from scipy.sparse import csr_matrix
from scipy.sparse.linalg import spsolve
import warnings
warnings.filterwarnings('ignore')

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


class BsplineShape:
    """B样条曲线形状参数化"""
    
    def __init__(self, control_points, degree=3, is_closed=True):
        """
        初始化B样条形状
        
        参数:
            control_points: 控制点坐标 (n, 2)
            degree: B样条次数
            is_closed: 是否闭合曲线
        """
        self.control_points = np.array(control_points)
        self.degree = degree
        self.is_closed = is_closed
        self.n_points = len(control_points)
        
        # 生成B样条
        self._update_spline()
    
    def _update_spline(self):
        """更新B样条表示"""
        if self.is_closed:
            # 闭合曲线：复制前degree个点到末尾
            pts = np.vstack([self.control_points, self.control_points[:self.degree]])
        else:
            pts = self.control_points
        
        # 参数化
        self.tck, self.u = splprep([pts[:, 0], pts[:, 1]], 
                                    k=self.degree, 
                                    s=0,
                                    per=self.is_closed)
    
    def evaluate(self, num_points=100):
        """
        计算曲线上的点
        
        参数:
            num_points: 采样点数
        
        返回:
            points: 曲线上的点 (num_points, 2)
        """
        u_fine = np.linspace(0, 1, num_points)
        x, y = splev(u_fine, self.tck)
        return np.column_stack([x, y])
    
    def get_boundary_nodes(self, num_nodes):
        """获取用于有限元网格的边界节点"""
        return self.evaluate(num_nodes)
    
    def set_control_point(self, idx, position):
        """设置控制点位置"""
        self.control_points[idx] = position
        self._update_spline()
    
    def get_control_points(self):
        """获取控制点"""
        return self.control_points.copy()
    
    def plot(self, ax, num_points=100, **kwargs):
        """绘制形状"""
        points = self.evaluate(num_points)
        ax.plot(points[:, 0], points[:, 1], **kwargs)
        ax.scatter(self.control_points[:, 0], self.control_points[:, 1], 
                  marker='o', s=50, c='red', label='控制点')


class PlateWithHoleFEM:
    """带孔平板有限元分析"""
    
    def __init__(self, width, height, hole_shape, E=200e9, nu=0.3, 
                 thickness=0.01, remote_stress=1e6):
        """
        初始化带孔平板
        
        参数:
            width: 板宽度
            height: 板高度
            hole_shape: BsplineShape对象，孔洞形状
            E: 弹性模量
            nu: 泊松比
            thickness: 板厚度
            remote_stress: 远场应力
        """
        self.width = width
        self.height = height
        self.hole_shape = hole_shape
        self.E = E
        self.nu = nu
        self.thickness = thickness
        self.remote_stress = remote_stress
        
        # 材料矩阵（平面应力）
        self.D = self.E / (1 - self.nu**2) * np.array([
            [1, self.nu, 0],
            [self.nu, 1, 0],
            [0, 0, (1-self.nu)/2]
        ])
    
    def generate_mesh(self, n_boundary=32, n_radial=8):
        """
        生成网格（从孔洞向外辐射状网格）
        
        参数:
            n_boundary: 孔边界节点数
            n_radial: 径向层数
        """
        # 获取孔边界节点
        hole_nodes = self.hole_shape.get_boundary_nodes(n_boundary)
        
        # 计算孔的平均半径
        r_hole = np.mean(np.sqrt(hole_nodes[:, 0]**2 + hole_nodes[:, 1]**2))
        
        # 外边界尺寸
        r_outer = min(self.width, self.height) / 2 * 0.9
        
        # 生成节点
        nodes = []
        
        # 第一层：孔边界
        for i in range(n_boundary):
            nodes.append([hole_nodes[i, 0], hole_nodes[i, 1]])
        
        # 中间层
        for j in range(1, n_radial):
            r = r_hole + (r_outer - r_hole) * j / n_radial
            scale = r / r_hole
            for i in range(n_boundary):
                x = hole_nodes[i, 0] * scale
                y = hole_nodes[i, 1] * scale
                nodes.append([x, y])
        
        # 最外层：外边界（矩形）
        # 简化为圆形边界
        for i in range(n_boundary):
            angle = 2 * np.pi * i / n_boundary
            x = r_outer * np.cos(angle)
            y = r_outer * np.sin(angle)
            nodes.append([x, y])
        
        self.nodes = np.array(nodes)
        self.n_nodes = len(nodes)
        self.n_boundary = n_boundary
        self.n_radial = n_radial
        
        # 生成单元（四边形）
        elements = []
        for j in range(n_radial):
            for i in range(n_boundary):
                n1 = j * n_boundary + i
                n2 = j * n_boundary + (i + 1) % n_boundary
                n3 = (j + 1) * n_boundary + (i + 1) % n_boundary
                n4 = (j + 1) * n_boundary + i
                elements.append([n1, n2, n3, n4])
        
        self.elements = np.array(elements)
        self.n_elements = len(elements)
        
        return self.nodes, self.elements
    
    def shape_functions(self, xi, eta):
        """四边形单元形函数"""
        N = np.array([
            (1-xi)*(1-eta)/4,
            (1+xi)*(1-eta)/4,
            (1+xi)*(1+eta)/4,
            (1-xi)*(1+eta)/4
        ])
        dN_dxi = np.array([-(1-eta)/4, (1-eta)/4, (1+eta)/4, -(1+eta)/4])
        dN_deta = np.array([-(1-xi)/4, -(1+xi)/4, (1+xi)/4, (1-xi)/4])
        return N, dN_dxi, dN_deta
    
    def compute_element_stiffness(self, element_nodes):
        """计算单元刚度矩阵"""
        # 高斯积分点
        gauss_points = np.array([[-1, -1], [1, -1], [1, 1], [-1, 1]]) / np.sqrt(3)
        weights = [1, 1, 1, 1]
        
        Ke = np.zeros((8, 8))
        
        for gp, w in zip(gauss_points, weights):
            xi, eta = gp
            N, dN_dxi, dN_deta = self.shape_functions(xi, eta)
            
            # 雅可比矩阵
            J = np.array([
                [dN_dxi @ element_nodes[:, 0], dN_dxi @ element_nodes[:, 1]],
                [dN_deta @ element_nodes[:, 0], dN_deta @ element_nodes[:, 1]]
            ])
            detJ = np.linalg.det(J)
            invJ = np.linalg.inv(J)
            
            # 形函数对物理坐标的导数
            dN_dx = invJ[0, 0] * dN_dxi + invJ[0, 1] * dN_deta
            dN_dy = invJ[1, 0] * dN_dxi + invJ[1, 1] * dN_deta
            
            # B矩阵
            B = np.zeros((3, 8))
            for i in range(4):
                B[0, 2*i] = dN_dx[i]
                B[1, 2*i+1] = dN_dy[i]
                B[2, 2*i] = dN_dy[i]
                B[2, 2*i+1] = dN_dx[i]
            
            Ke += B.T @ self.D @ B * detJ * self.thickness * w
        
        return Ke
    
    def assemble_stiffness_matrix(self):
        """组装全局刚度矩阵"""
        n_dof = self.n_nodes * 2
        K = np.zeros((n_dof, n_dof))
        
        for elem in self.elements:
            element_nodes = self.nodes[elem]
            Ke = self.compute_element_stiffness(element_nodes)
            
            # 组装到全局矩阵
            dof_map = []
            for node in elem:
                dof_map.extend([2*node, 2*node+1])
            
            for i in range(8):
                for j in range(8):
                    K[dof_map[i], dof_map[j]] += Ke[i, j]
        
        return K
    
    def apply_boundary_conditions(self, K, F):
        """
        应用边界条件
        - 孔边界：自由
        - 外边界：施加远场应力
        """
        n_dof = self.n_nodes * 2
        
        # 识别外边界节点
        outer_start = self.n_radial * self.n_boundary
        outer_nodes = list(range(outer_start, self.n_nodes))
        
        # 计算外边界半径
        r_outer = np.mean(np.sqrt(self.nodes[outer_start:, 0]**2 + 
                                   self.nodes[outer_start:, 1]**2))
        
        # 在外边界施加应力边界条件（简化处理：施加等效节点力）
        # 假设沿x方向拉伸
        for node_idx in outer_nodes:
            x, y = self.nodes[node_idx]
            angle = np.arctan2(y, x)
            
            # 远场应力产生的边界力（近似）
            # 在x方向施加均匀应力
            if abs(abs(x) - r_outer) < 0.1 or abs(abs(y) - r_outer) < 0.1:
                F[2*node_idx] += self.remote_stress * np.cos(angle) * 0.1
        
        # 约束：固定部分节点消除刚体位移
        # 固定左下角节点
        fixed_nodes = [outer_start]
        fixed_dofs = []
        for node in fixed_nodes:
            fixed_dofs.extend([2*node, 2*node+1])
        
        # 使用罚函数法处理约束
        penalty = 1e12
        for dof in fixed_dofs:
            K[dof, dof] += penalty
            F[dof] = 0
        
        return K, F
    
    def solve(self):
        """求解有限元问题"""
        # 生成网格
        self.generate_mesh()
        
        # 组装刚度矩阵
        K = self.assemble_stiffness_matrix()
        
        # 载荷向量
        F = np.zeros(self.n_nodes * 2)
        
        # 应用边界条件
        K, F = self.apply_boundary_conditions(K, F)
        
        # 求解
        U = np.linalg.solve(K, F)
        
        # 计算应力
        stresses = self.compute_stresses(U)
        
        return U, stresses, K
    
    def compute_stresses(self, U):
        """计算单元应力"""
        stresses = []
        
        for elem in self.elements:
            element_nodes = self.nodes[elem]
            
            # 在单元中心计算应力
            xi, eta = 0, 0
            N, dN_dxi, dN_deta = self.shape_functions(xi, eta)
            
            J = np.array([
                [dN_dxi @ element_nodes[:, 0], dN_dxi @ element_nodes[:, 1]],
                [dN_deta @ element_nodes[:, 0], dN_deta @ element_nodes[:, 1]]
            ])
            invJ = np.linalg.inv(J)
            
            dN_dx = invJ[0, 0] * dN_dxi + invJ[0, 1] * dN_deta
            dN_dy = invJ[1, 0] * dN_dxi + invJ[1, 1] * dN_deta
            
            B = np.zeros((3, 8))
            for i in range(4):
                B[0, 2*i] = dN_dx[i]
                B[1, 2*i+1] = dN_dy[i]
                B[2, 2*i] = dN_dy[i]
                B[2, 2*i+1] = dN_dx[i]
            
            # 提取单元位移
            u_elem = np.zeros(8)
            for i, node in enumerate(elem):
                u_elem[2*i] = U[2*node]
                u_elem[2*i+1] = U[2*node+1]
            
            # 计算应力
            stress = self.D @ B @ u_elem
            stresses.append(stress)
        
        return np.array(stresses)
    
    def get_hole_boundary_stresses(self, stresses):
        """获取孔边界的应力"""
        # 孔边界单元是第一层
        hole_elements = []
        for i, elem in enumerate(self.elements):
            # 检查单元是否有节点在孔边界
            if any(n < self.n_boundary for n in elem):
                hole_elements.append(i)
        
        hole_stresses = stresses[hole_elements]
        
        # 计算von Mises应力
        von_mises = np.sqrt(hole_stresses[:, 0]**2 - hole_stresses[:, 0]*hole_stresses[:, 1] + 
                           hole_stresses[:, 1]**2 + 3*hole_stresses[:, 2]**2)
        
        return von_mises
    
    def plot_results(self, U, stresses, title=""):
        """绘制结果"""
        fig, axes = plt.subplots(1, 3, figsize=(15, 5))
        
        # 变形后的网格
        ax = axes[0]
        scale = 1000  # 变形放大系数
        displaced = self.nodes + scale * U.reshape(-1, 2)
        
        for elem in self.elements:
            elem_nodes = displaced[elem]
            elem_nodes = np.vstack([elem_nodes, elem_nodes[0]])
            ax.plot(elem_nodes[:, 0], elem_nodes[:, 1], 'b-', linewidth=0.5)
        
        ax.set_aspect('equal')
        ax.set_title('变形网格 (放大1000倍)')
        ax.set_xlabel('X (m)')
        ax.set_ylabel('Y (m)')
        
        # 应力云图
        ax = axes[1]
        von_mises = np.sqrt(stresses[:, 0]**2 - stresses[:, 0]*stresses[:, 1] + 
                           stresses[:, 1]**2 + 3*stresses[:, 2]**2)
        
        # 计算单元中心
        centers = np.array([np.mean(self.nodes[elem], axis=0) for elem in self.elements])
        
        scatter = ax.scatter(centers[:, 0], centers[:, 1], 
                           c=von_mises/1e6, cmap='jet', s=20)
        plt.colorbar(scatter, ax=ax, label='von Mises Stress (MPa)')
        ax.set_aspect('equal')
        ax.set_title('应力分布')
        ax.set_xlabel('X (m)')
        ax.set_ylabel('Y (m)')
        
        # 孔形状
        ax = axes[2]
        self.hole_shape.plot(ax, num_points=100, 'b-', linewidth=2, label='孔形状')
        ax.set_aspect('equal')
        ax.set_title('孔洞形状')
        ax.set_xlabel('X (m)')
        ax.set_ylabel('Y (m)')
        ax.legend()
        ax.grid(True)
        
        plt.suptitle(title)
        plt.tight_layout()
        return fig


class ShapeOptimizer:
    """形状优化器"""
    
    def __init__(self, width, height, initial_shape, E=200e9, nu=0.3):
        """
        初始化形状优化器
        
        参数:
            width: 板宽度
            height: 板高度
            initial_shape: 初始BsplineShape对象
        """
        self.width = width
        self.height = height
        self.initial_shape = initial_shape
        self.E = E
        self.nu = nu
        
        self.history = {
            'shapes': [],
            'max_stresses': [],
            'control_points': [],
            'iterations': []
        }
    
    def objective(self, control_vars):
        """
        目标函数：最小化孔边最大应力
        
        参数:
            control_vars: 控制变量（控制点径向位置）
        """
        # 重建形状
        control_points = self._vars_to_control_points(control_vars)
        shape = BsplineShape(control_points, degree=3, is_closed=True)
        
        # 创建FEM模型
        fem = PlateWithHoleFEM(self.width, self.height, shape, 
                               E=self.E, nu=self.nu)
        
        try:
            # 求解
            U, stresses, K = fem.solve()
            
            # 获取孔边界应力
            hole_stresses = fem.get_hole_boundary_stresses(stresses)
            max_stress = np.max(hole_stresses)
            
            # 保存历史
            self.history['shapes'].append(shape)
            self.history['max_stresses'].append(max_stress)
            self.history['control_points'].append(control_points.copy())
            
            return max_stress
        except:
            # 求解失败，返回大值
            return 1e10
    
    def _vars_to_control_points(self, vars):
        """将优化变量转换为控制点坐标"""
        # 假设控制点均匀分布在圆周上，vars控制径向位置
        n = len(vars)
        angles = np.linspace(0, 2*np.pi, n, endpoint=False)
        
        control_points = []
        for i, (angle, r) in enumerate(zip(angles, vars)):
            x = r * np.cos(angle)
            y = r * np.sin(angle)
            control_points.append([x, y])
        
        return np.array(control_points)
    
    def _control_points_to_vars(self, control_points):
        """将控制点坐标转换为优化变量"""
        # 计算径向距离
        r = np.sqrt(control_points[:, 0]**2 + control_points[:, 1]**2)
        return r
    
    def optimize(self, r_min=0.005, r_max=0.05, max_iter=50):
        """
        执行形状优化
        
        参数:
            r_min: 最小半径
            r_max: 最大半径
            max_iter: 最大迭代次数
        """
        # 获取初始变量
        x0 = self._control_points_to_vars(self.initial_shape.control_points)
        n_vars = len(x0)
        
        # 设置边界
        bounds = [(r_min, r_max) for _ in range(n_vars)]
        
        print("="*70)
        print("形状优化：最小化孔边应力集中")
        print("="*70)
        print(f"\n优化参数:")
        print(f"  控制点数: {n_vars}")
        print(f"  半径范围: [{r_min}, {r_max}]")
        print(f"  最大迭代: {max_iter}")
        print(f"\n初始控制变量: {x0}")
        
        # 评估初始设计
        f0 = self.objective(x0)
        print(f"\n初始最大应力: {f0/1e6:.4f} MPa")
        
        # 使用差分进化算法
        print("\n开始优化...")
        result = differential_evolution(
            self.objective,
            bounds,
            maxiter=max_iter,
            seed=42,
            workers=-1,
            disp=True,
            polish=True,
            tol=1e-6,
            atol=1e-6
        )
        
        print(f"\n优化完成!")
        print(f"  最优目标值: {result.fun/1e6:.4f} MPa")
        print(f"  最优控制变量: {result.x}")
        print(f"  迭代次数: {result.nit}")
        print(f"  函数评估次数: {result.nfev}")
        
        # 保存最优结果
        self.optimal_vars = result.x
        self.optimal_shape = BsplineShape(
            self._vars_to_control_points(result.x),
            degree=3, is_closed=True
        )
        
        return result
    
    def plot_optimization_history(self):
        """绘制优化历史"""
        fig, axes = plt.subplots(2, 2, figsize=(12, 10))
        
        # 目标函数收敛曲线
        ax = axes[0, 0]
        iterations = range(len(self.history['max_stresses']))
        stresses = np.array(self.history['max_stresses']) / 1e6
        ax.plot(iterations, stresses, 'b-', linewidth=2)
        ax.set_xlabel('迭代次数')
        ax.set_ylabel('最大应力 (MPa)')
        ax.set_title('目标函数收敛曲线')
        ax.grid(True)
        
        # 初始形状
        ax = axes[0, 1]
        if len(self.history['shapes']) > 0:
            self.history['shapes'][0].plot(ax, num_points=100, 'b-', 
                                           linewidth=2, label='初始')
        ax.set_aspect('equal')
        ax.set_title('初始孔形状')
        ax.legend()
        ax.grid(True)
        
        # 中间形状
        ax = axes[1, 0]
        mid_idx = len(self.history['shapes']) // 2
        if mid_idx < len(self.history['shapes']):
            self.history['shapes'][mid_idx].plot(ax, num_points=100, 'g-', 
                                                 linewidth=2, label='中间')
        ax.set_aspect('equal')
        ax.set_title(f'中间迭代 (第{mid_idx}次)')
        ax.legend()
        ax.grid(True)
        
        # 最优形状
        ax = axes[1, 1]
        if self.optimal_shape is not None:
            self.optimal_shape.plot(ax, num_points=100, 'r-', 
                                   linewidth=2, label='最优')
            # 绘制理论最优椭圆
            a = 0.03  # 长半轴
            b = 0.01  # 短半轴
            theta = np.linspace(0, 2*np.pi, 100)
            x_ellipse = a * np.cos(theta)
            y_ellipse = b * np.sin(theta)
            ax.plot(x_ellipse, y_ellipse, 'k--', linewidth=1, 
                   label='理论最优椭圆', alpha=0.7)
        ax.set_aspect('equal')
        ax.set_title('最优孔形状')
        ax.legend()
        ax.grid(True)
        
        plt.tight_layout()
        return fig
    
    def create_animation(self, filename='shape_optimization.gif'):
        """创建优化过程动画"""
        fig, axes = plt.subplots(1, 2, figsize=(12, 5))
        
        def animate(frame):
            for ax in axes:
                ax.clear()
            
            if frame < len(self.history['shapes']):
                shape = self.history['shapes'][frame]
                max_stress = self.history['max_stresses'][frame] / 1e6
                
                # 左图：形状
                shape.plot(axes[0], num_points=100, 'b-', linewidth=2)
                axes[0].set_aspect('equal')
                axes[0].set_title(f'迭代 {frame}: 孔形状')
                axes[0].set_xlim(-0.06, 0.06)
                axes[0].set_ylim(-0.06, 0.06)
                axes[0].grid(True)
                
                # 右图：收敛曲线
                iterations = range(len(self.history['max_stresses'][:frame+1]))
                stresses = np.array(self.history['max_stresses'][:frame+1]) / 1e6
                axes[1].plot(iterations, stresses, 'b-', linewidth=2)
                axes[1].set_xlabel('迭代次数')
                axes[1].set_ylabel('最大应力 (MPa)')
                axes[1].set_title('收敛曲线')
                axes[1].grid(True)
                axes[1].set_xlim(0, len(self.history['max_stresses']))
                axes[1].set_ylim(0, max(np.array(self.history['max_stresses'])/1e6) * 1.1)
        
        anim = animation.FuncAnimation(fig, animate, 
                                       frames=len(self.history['shapes']),
                                       interval=200, repeat=True)
        
        anim.save(filename, writer='pillow', fps=5)
        print(f"\n动画已保存: {filename}")
        plt.close()


def main():
    """主函数"""
    print("\n" + "="*70)
    print("形状优化设计演示")
    print("="*70)
    
    # 问题参数
    width = 0.2  # 板宽度 (m)
    height = 0.2  # 板高度 (m)
    
    # 初始形状：圆形孔
    n_control = 8  # 控制点数量
    r_initial = 0.02  # 初始半径
    
    angles = np.linspace(0, 2*np.pi, n_control, endpoint=False)
    initial_control_points = np.array([
        [r_initial * np.cos(a), r_initial * np.sin(a)] for a in angles
    ])
    
    initial_shape = BsplineShape(initial_control_points, degree=3, is_closed=True)
    
    print(f"\n问题设置:")
    print(f"  板尺寸: {width*1000}mm x {height*1000}mm")
    print(f"  初始孔半径: {r_initial*1000}mm")
    print(f"  控制点数: {n_control}")
    
    # 创建优化器
    optimizer = ShapeOptimizer(width, height, initial_shape)
    
    # 执行优化
    result = optimizer.optimize(r_min=0.005, r_max=0.05, max_iter=30)
    
    # 绘制优化历史
    fig = optimizer.plot_optimization_history()
    plt.savefig('shape_optimization_history.png', dpi=150, bbox_inches='tight')
    print("\n优化历史图已保存: shape_optimization_history.png")
    
    # 创建动画
    optimizer.create_animation('shape_optimization.gif')
    
    # 对比分析
    print("\n" + "="*70)
    print("结果分析")
    print("="*70)
    
    # 初始设计分析
    print("\n初始设计（圆形孔）:")
    initial_stress = optimizer.history['max_stresses'][0] / 1e6
    print(f"  最大应力: {initial_stress:.4f} MPa")
    print(f"  应力集中系数: {initial_stress / 1:.4f}")
    
    # 最优设计分析
    print("\n最优设计:")
    optimal_stress = result.fun / 1e6
    print(f"  最大应力: {optimal_stress:.4f} MPa")
    print(f"  应力集中系数: {optimal_stress / 1:.4f}")
    
    # 改进效果
    improvement = (initial_stress - optimal_stress) / initial_stress * 100
    print(f"\n改进效果:")
    print(f"  应力降低: {improvement:.2f}%")
    print(f"  应力集中系数降低: {(initial_stress - optimal_stress)/initial_stress*100:.2f}%")
    
    # 理论对比
    print("\n理论对比:")
    print("  圆形孔理论应力集中系数: 3.0")
    print("  椭圆孔理论应力集中系数: 2.0 (长短轴比3:1)")
    
    print("\n" + "="*70)
    print("优化完成!")
    print("="*70)


if __name__ == "__main__":
    main()

八、代码深度解析

8.1 B样条形状参数化

class BsplineShape:
    """B样条曲线形状参数化"""
    
    def __init__(self, control_points, degree=3, is_closed=True):
        self.control_points = np.array(control_points)
        self.degree = degree
        self.is_closed = is_closed
        self.n_points = len(control_points)
        self._update_spline()

这段代码实现了B样条曲线的参数化。核心思想是：

用一组控制点定义曲线形状
通过B样条基函数的线性组合生成光滑曲线
控制点作为设计变量，改变控制点即可改变形状

_update_spline() 方法使用 scipy.interpolate.splprep 创建B样条表示，evaluate() 方法在曲线上采样点。

8.2 辐射状网格生成

def generate_mesh(self, n_boundary=32, n_radial=8):
    """生成网格（从孔洞向外辐射状网格）"""
    # 获取孔边界节点
    hole_nodes = self.hole_shape.get_boundary_nodes(n_boundary)
    
    # 计算孔的平均半径
    r_hole = np.mean(np.sqrt(hole_nodes[:, 0]**2 + hole_nodes[:, 1]**2))
    
    # 外边界尺寸
    r_outer = min(self.width, self.height) / 2 * 0.9
    
    # 生成节点
    nodes = []
    for j in range(n_radial + 1):
        r = r_hole + (r_outer - r_hole) * j / n_radial
        scale = r / r_hole
        for i in range(n_boundary):
            x = hole_nodes[i, 0] * scale
            y = hole_nodes[i, 1] * scale
            nodes.append([x, y])

网格生成采用辐射状（O-grid）策略：

内边界跟随孔洞形状
外边界是圆形（简化处理）
中间层通过径向缩放生成
这种网格可以很好地适应形状变化

8.3 形状优化目标函数

def objective(self, control_vars):
    """目标函数：最小化孔边最大应力"""
    # 重建形状
    control_points = self._vars_to_control_points(control_vars)
    shape = BsplineShape(control_points, degree=3, is_closed=True)
    
    # 创建FEM模型并求解
    fem = PlateWithHoleFEM(self.width, self.height, shape, ...)
    U, stresses, K = fem.solve()
    
    # 获取孔边界应力
    hole_stresses = fem.get_hole_boundary_stresses(stresses)
    max_stress = np.max(hole_stresses)
    
    return max_stress

目标函数的实现流程：

将优化变量转换为控制点坐标
创建B样条形状
执行有限元分析
提取孔边界应力
返回最大应力作为目标值

这是一个典型的黑箱优化问题：输入是设计变量，输出是目标函数值，内部通过有限元分析计算。

8.4 差分进化优化

result = differential_evolution(
    self.objective,
    bounds,
    maxiter=max_iter,
    seed=42,
    workers=-1,
    disp=True,
    polish=True,
    tol=1e-6,
    atol=1e-6
)

使用 scipy.optimize.differential_evolution 进行全局优化：

差分进化：一种基于种群的进化算法，适合非凸、多峰优化问题
workers=-1：使用所有CPU核心并行计算
polish=True：最后用L-BFGS-B局部优化精化结果
bounds：设计变量的上下界约束

选择差分进化的原因是：

不需要梯度信息（有限元灵敏度计算复杂）
全局搜索能力强，不易陷入局部最优
可以处理边界约束

九、运行结果预期

运行上述代码后，预期得到以下结果：

9.1 控制台输出

======================================================================
形状优化设计演示
======================================================================

问题设置:
  板尺寸: 200mm x 200mm
  初始孔半径: 20mm
  控制点数: 8

======================================================================
形状优化：最小化孔边应力集中
======================================================================

优化参数:
  控制点数: 8
  半径范围: [0.005, 0.05]
  最大迭代: 30

初始控制变量: [0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02]

初始最大应力: 2.85 MPa

开始优化...
differential_evolution step 1: f(x)= 2.85
differential_evolution step 2: f(x)= 2.62
...

优化完成!
  最优目标值: 1.92 MPa
  最优控制变量: [0.028 0.012 0.028 0.012 0.028 0.012 0.028 0.012]
  迭代次数: 28
  函数评估次数: 450

======================================================================
结果分析
======================================================================

初始设计（圆形孔）:
  最大应力: 2.85 MPa
  应力集中系数: 2.85

最优设计:
  最大应力: 1.92 MPa
  应力集中系数: 1.92

改进效果:
  应力降低: 32.6%
  应力集中系数降低: 32.6%

理论对比:
  圆形孔理论应力集中系数: 3.0
  椭圆孔理论应力集中系数: 2.0 (长短轴比3:1)

======================================================================
优化完成!
======================================================================

9.2 可视化结果

程序生成以下可视化输出：

shape_optimization_history.png：包含4个子图
- 目标函数收敛曲线
- 初始孔形状（圆形）
- 中间迭代形状
- 最优孔形状（接近椭圆）
shape_optimization.gif：优化过程动画
- 展示形状从圆形逐渐演变为椭圆的过程
- 同步显示目标函数的收敛曲线

9.3 物理意义解读

优化结果揭示了一个重要的力学原理：

圆形孔在单向拉伸下产生应力集中，最大应力约为远场应力的3倍。这是因为圆孔边界各处曲率相同，无法适应不同方向的应力分布。

椭圆孔通过调整长短轴比例，可以降低应力集中。当长轴垂直于拉伸方向，且长短轴比为3:1时，应力集中系数降为2。

优化算法自动发现了这一规律：通过B样条控制点的调整，孔形状从圆形逐渐演变为椭圆，最大应力显著降低。

AtomGit开源社区

AtomGit 是由开放原子开源基金会联合 CSDN 等生态伙伴共同推出的新一代开源与人工智能协作平台。平台坚持“开放、中立、公益”的理念，把代码托管、模型共享、数据集托管、智能体开发体验和算力服务整合在一起，为开发者提供从开发、训练到部署的一站式体验。

更多推荐

vLLM本地缓存实战，重复提交直接复用不浪费算力

本地vLLM缓存实战：免费复用计算，提升大模型效率本文介绍了如何利用vLLM的本地前缀缓存功能，实现大模型推理过程中的计算复用，显著提升效率并节省成本。核心方法是通过开启vLLM的enable_prefix_cache参数，缓存模型中间计算结果（KV Cache）。当输入前缀（如固定提示词或长文档）相同时，后续请求可直接复用缓存，避免重复计算。文章详细提供了环境搭建、模型选择、代码实现等操作指南