基于元胞自动机的结构拓扑优化方法

# 基于元胞自动机的结构拓扑优化方法 ## 1. 复合层合板连续域 二维连续体的一个特殊情况是正交各向异性纤维增强复合层合板，其中纤维在特定方向上提供刚度和强度。因此，确定纤维取向角是设计中的重要环节，这通常被称为定制。对于拓扑和纤维角的组合设计，各向同性连续域的基本元胞自动机（CA）元素几乎保持不变，但第 $i$ 个单元的状态修改如下： \[ \Sigma(i) = \left\{ u_{i}^{(1...m)}, f_{i}^{(1...m)}, (\rho_i, \theta_i), Q_i \right\} \] 其中，$\theta_i$ 是第 $i$ 个单元的纤维角度，$Q_i$ 是简化的变换刚度矩阵，由于对称性，仅存储矩阵的上半对角线。 ## 2. 分析更新规则 ### 2.1 桁架结构 局部分析规则源自单元与其相邻单元的平衡条件。在一个单元内，邻域结构的每个桁架成员（$k = 1, ..., 8$）具有杨氏模量 $E$、变形前的长度 $L_{k}^{i}$ 和横截面积 $A_{k}^{i}$。与该单元相关的总势能是邻域结构中八个桁架成员的应变能之和，以及施加在该单元上的外力势能： \[ \Pi_i = \sum_{k = 1}^{8} \frac{E A_{k}^{i} L_{k}^{i} (\varepsilon_{k}^{i})^2}{2} - f_{x}^{i} u_i - f_{y}^{i} v_i \] 其中，$\varepsilon_{k}^{i}$ 是桁架成员的应变，它取决于相邻单元的相对位移。应变使用桁架成员的格林应变定义进行评估： \[ \varepsilon_{k}^{i} = \frac{(u_i - u_{k}^{i}) \cos \theta_{k} - (v_i - v_{k}^{i}) \sin \theta_{k}}{L_{k}^{i}} \] 其中，$(u_{k}^{i}, v_{k}^{i})$ 是相邻单元的位移，$\theta_{k}$ 是第 $k$ 个桁架元素成员相对于单元中心的取向角。 通过相对于单元位移 $u_i$ 和 $v_i$ 最小化总势能，可得到平衡方程： \[ \frac{\partial \Pi_i}{\partial u_i} = 0, \frac{\partial \Pi_i}{\partial v_i} = 0 \] ### 2.2 各向同性连续体结构 邻域结构的平衡（如图所示）再次用于制定局部分析更新规则。与一个单元相关的总势能是邻域结构中每个元素的应变能之和，以及直接施加在该单元上的外力势能。在这种离散化方案中，只有少数单元会受到外力作用。大多数单元的平衡仅涉及局部单元通过它们之间的固体域进行相互作用，在本章中，这些相互作用的区域被称为元素： \[ \Pi_i = \sum_{k = 1}^{N_{element}} U_{k}^{i} - \mathbf{f}_i \cdot \mathbf{u}_i \] 其中，$N_{element}$ 是围绕一个单元的元素数量，$U_{k}^{i}$ 是第 $k$ 个元素的应变能，$\mathbf{f}_i$ 是施加的力向量，$\mathbf{u}_i$ 是包括该单元及其邻域所有单元的位移向量。 元素的应变能表示为基础材料应变能的形式： \[ U_{k} = \bar{\rho}^p \tilde{U} \] 其中， \[ \tilde{U} = \frac{1}{2} \int_{element} \boldsymbol{\epsilon}^T \cdot \mathbf{Q} \cdot \boldsymbol{\epsilon} dx dy dz \] 是基础材料的应变能，$\boldsymbol{\epsilon}$ 是小应变张量，$\mathbf{Q}$ 是简化的面内刚度矩阵。方程中的符号 $p$ 称为惩罚参数，用于设计目的，其作用将在后续设计规则中解释。 元素密度 $\bar{\rho}$ 通过平均密度插值获得： \[ \frac{1}{\bar{\rho}^p} = \frac{1}{N_{cell}} \sum_{i = 1}^{N_{cell}} \frac{1}{\rho_{i}^p} \] 其中，$\rho_i$ 是围绕元素的单元的密度度量，$N_{cell}$ 是定义元素的单元数量。对于二维邻域结构，$N_{cell} = 4$；对于三维邻域结构，$N_{cell} = 8$。 该密度插值方案的选择使得任何密度度量低于阈值的节点都会关闭其参与的所有四个元素。使用此方案，在优化过程中会自动抑制棋盘格图案。 通过相对于单元位移最小化总势能，可得到平衡方程： \[ \min_{\mathbf{u}_C} \Pi_i \] 每个单元的平衡方程以残差形式表示： \[ \mathbf{R}_C(\mathbf{u}_C, \mathbf{u}_N) = \begin{bmatrix} \mathbf{G}_C(\mathbf{u}_C, \mathbf{u}_N) \\ \mathbf{G}_N(\mathbf{u}_C, \mathbf{u}_N) \end{bmatrix} + \begin{bmatrix} \mathbf{f}_C \\ \mathbf{f}_N \end{bmatrix} = 0 \] 其中，$\mathbf{u}_C$ 和 $\mathbf{u}_N$ 分别是单元和邻域的位移向量，$\mathbf{G}_C$ 和 $\mathbf{G}_N$ 是内力向量，$\mathbf{f}_C$ 和 $\mathbf{f}_N$ 分别是相对于单元的施加力向量和相对于邻域的内力向量。 对向量 $\mathbf{R}_C$ 相对于 $\mathbf{u}_C$ 的分量求导，线性刚度矩阵可写为： \[ \mathbf{K} = - \frac{\partial \mathbf{R}_C}{\partial \mathbf{u}_C} (\mathbf{u}_C, \mathbf{u}_N) \] 刚度矩阵 $\mathbf{K}$ 也可以表示为总势能的海森矩阵： \[ K_{pq} = \frac{\partial^2 \Pi_i}{\partial u_p \partial u_q} \] 因此，单元位移更新如下： \[ \mathbf{u}_{C}^{(t + 1)} = \mathbf{u}_{C}^{(t)} + \Delta \mathbf{u}_C \] \[ \Delta \mathbf{u}_C = (\mathbf{K}_C)^{-1} \cdot \left( \mathbf{G}_C(\mathbf{u}_{N}^{(t + 1)}) + \mathbf{f}_C \right) \] 其中，$\mathbf{K}_C$ 是二维或三维情况下的 $(2 \times 2)$ 或 $(3 \times 3)$ 单元刚度矩阵。 ### 2.3 复合层合板连续体结构 在拓扑优