一种考虑损伤的自由曲面形态创建方法.pdf

本发明提供一种考虑损伤的自由曲面形态创建方法，考虑到随着荷载的施加和变形的积累，材料在自由曲面创建过程中出现的微观损伤：1.考虑损伤的自由曲面创建能够保证自由曲面的受力性能和整体刚度处于局部最优；2.现有创建方法中不考虑自由曲面的损伤会造成的自由曲面创建结果偏于不安全，通过以上技术方案可以改进这一缺点。

1.一种考虑损伤的自由曲面形态创建方法，其特征在于，具体方法步骤如下：
步骤1，根据给定的初始自由曲面的一组控制点坐标初始化自由曲面；
步骤2，对步骤1得到的初始化自由曲面进行网格划分，得到三角形薄壳单元；
步骤3，计算步骤2得到的三角形薄壳单元的单元坐标系下的刚度矩阵；
步骤4，计算整体坐标下的单元刚度矩阵；
步骤5，经过总纲集成得到自由曲面结构的整体刚度矩阵；
步骤6，加入Mazars损伤模型，对结构施加荷载，分N个荷载步进行加载，每个荷载步内
存在M个迭代步，计算自由曲面壳体位移；
步骤7，采用应变能作为自由曲面结构创建的优化目标，节点坐标作为自由曲面结构创
建的优化变量，通过节点荷载和节点位移计算应变能；
步骤8，对考虑Mazars损伤的自由曲面的控制点进行高度调整；
步骤9，设置优化精度ε^*，若当前优化步的应变能对当前优化步的控制点高度的差分的
模小于优化精度，则考虑损伤的自由曲面优化结束，自由曲面创建完成；否则继续迭代进行
自由曲面结构形态优化。
2.根据权利要求1所述的一种考虑损伤的自由曲面形态创建方法，其特征在于，步骤1
中的自由曲面是NURBS曲面。
3.根据权利要求2所述的一种考虑损伤的自由曲面形态创建方法，其特征在于，自由曲
面是一张在u方向p次、v方向q次的NURBS曲面，具有如下形式的双变量分段有理矢值函数：
$S\left(u,v\right)=\underset{i=0}{\overset{n}{\mathrm{\Σ}}}\underset{j=0}{\overset{m}{\mathrm{\Σ}}}{N}_{i,p}\left(u\right)N_{j,q}\left(v\right)w_{i,j}{P}_{i,j}/\underset{i=0}{\overset{n}{\mathrm{\Σ}}}\underset{j=0}{\overset{m}{\mathrm{\Σ}}}{N}_{i,p}\left(u\right)N_{j,q}\left(v\right)P_{i,j}$
其中，n为u方向的控制点个数，且i∈[1,n]；m为v方向的控制点个数，且j∈[1,m]；P_i,j
是u方向编号为i且v方向编号为j的控制点坐标；w_i,j是u方向编号为i且v方向编号为j的控
制点的全因子；N_i,p(u)和N_j,q(v)是分别定义在矢量U和V上的样条基函数，u_i为u方向编号为
i的节点矢量且u_i∈U，v_j为v方向编号为j的
节点矢量且v_j∈V，
N_i,p(u)和N_j,q(v)的表达式如下所示：
$N_{i,p}\left(u\right)=\frac{u-u_i}{u_{i+p}-u_i}{N}_{i,p-1}\left(u\right)+\frac{u_{i+p+1}-u}{u_{i+p+1}-u_{i+1}}{N}_{i+1,p-1}\left(u\right)$
$N_{j,q}\left(v\right)=\frac{v-v_j}{v_{j+q}-v_j}{N}_{j,q-1}\left(v\right)+\frac{v_{j+q+1}-v}{v_{j+q+1}-v_{j+1}}{N}_{j+1,q-1}\left(v\right).$
4.根据权利要求1所述的一种考虑损伤的自由曲面形态创建方法，其特征在于，步骤3
中计算三角形薄壳单元的单元坐标系下的刚度矩阵，具体为：
三角形薄壳单元的刚度矩阵由平面应力状态和弯曲应力状态的刚度矩阵得到，三角形
薄壳单元的三个节点分别记为r、s、t；
在三角形平面应力单元中，节点力和节点位移的关系如下：
$\left\{\begin{array}{c}{F}_r^p\\ F_s^p\\ F_t^p\end{array}\right\}=\[k^p\]\left\{\begin{array}{c}{\stackrel{\‾}{\mathrm{\δ}}}_r^p\\ {\stackrel{\‾}{\mathrm{\δ}}}_s^p\\ {\stackrel{\‾}{\mathrm{\δ}}}_t^p\end{array}\right\}$
其中，[k^p]表示单元坐标系下三角形平面应力单元的刚度矩阵；为单元坐标系下平
面应力单元的r节点的节点力，且由u方向的节点力U_r和v方向的节点力V_r组成，
为单元坐标系下平面应力单元的s节点的节点力，且由u方向的节点力U_s和
v方向的节点力V_s组成，为单元坐标系下平面应力单元的t节点的节点力，
且由u方向的节点力U_t和v方向的节点力V_t组成，表示单元坐标系下三角形
平面应力单元的r节点的位移，且由u方向的位移u_r和v方向的位移v_r组成，
表示单元坐标系下三角形平面应力单元的s节点的位移，且由u方向的位移u_s和v方向的位
移v_s组成，表示单元坐标系下三角形平面应力单元的t节点的位移，且由u
方向的位移u_t和v方向的位移v_t组成，
在三角形弯曲应力单元中，节点力和节点位移的关系如下：
$\left\{\begin{array}{c}{F}_r^b\\ F_s^b\\ F_t^b\end{array}\right\}=\[k^b\]\left\{\begin{array}{c}{\stackrel{\‾}{\mathrm{\δ}}}_r^b\\ {\stackrel{\‾}{\mathrm{\δ}}}_s^b\\ {\stackrel{\‾}{\mathrm{\δ}}}_t^b\end{array}\right\}$
其中，[k^b]表示单元坐标系下三角形弯曲板单元的刚度矩阵；为单元坐标系下弯曲
板单元r节点的节点力，且由垂直于u、v平面方向的节点力W_r、绕u轴旋转的节点力M_ur及绕v
轴旋转的节点力M_vr组成，为单元坐标系下弯曲板单元s节点的节点力，且
由垂直于u、v平面方向的节点力W_s、绕u轴旋转的节点力M_us及绕v轴旋转的节点力M_vs组成，
为单元坐标系下弯曲板单元t节点的节点力，且由垂直于u、v平面方向的
节点力W_t、绕u轴旋转的节点力M_ut及绕v轴旋转的节点力M_vt组成，表示单
元坐标系下三角形弯曲板单元的r节点的位移，垂直于u、v平面方向的线位移w_r、绕u轴旋转
角θ_ur及绕v轴旋转角θ_vr组成，表示单元坐标系下三角形弯曲板单元的s节
点的位移，垂直于u、v平面方向的线位移w_s、绕u轴旋转角θ_us及绕v轴旋转角θ_vs组成，
表示单元坐标系下三角形弯曲板单元的t节点的位移，垂直于u、v平面方
向的线位移w_t、绕u轴旋转角θ_ut及绕v轴旋转角θ_vt组成，
三角形薄壳单元的单元坐标系下的刚度矩阵为：
$\[k\]=\left[\begin{array}{cccccc}\[k^p\]& & 0& 0& 0& 0\\ & & 0& 0& 0& 0\\ 0& 0& & & & 0\\ 0& 0& & \[k^b\]& & 0\\ 0& 0& & & & 0\\ 0& 0& 0& 0& 0& 0\end{array}\right].$
5.根据权利要求4所述的一种考虑损伤的自由曲面形态创建方法，其特征在于，步骤4
中整体坐标下的单元刚度矩阵[k′]为：
[k′]＝[L]^-1[k][L]
其中，[L]为坐标转换矩阵，{δ_r}＝[λ]{δ′_r}，{δ′_r}＝
[u′_r v′_r w′_r θ′_xr θ′_yr θ′_zr]^T，θ_r表示单元坐标系下r节点绕u、v面法线的转角，u′_r、v′_r、w
′_r、θ′_xr、θ′_yr、θ′_zr分别表示整体坐标系下的三角形单元的r节点在整体坐标系下x、y、z方向
的位移以及绕x、y、z方向的旋转角；{F_r}＝[λ]{F′_r}，M_θr表示单元坐标系下r
节点的虚拟弯矩，{F′_r}＝[U′_r V′_r W′_r M′_θxr M′_θyr M′_θzr]^T，U′_r、V′_r、W′_r、M′_θxr、M′_θyr、M′_θzr分
别表示整体坐标系下的三角形单元的r节点在整体坐标系下x、y、z方向的节点力以及绕x、
y、z方向的旋转角。
6.根据权利要求5所述的一种考虑损伤的自由曲面形态创建方法，其特征在于，步骤6
中计算自由曲面壳体位移的具体步骤如下：
(6-1)计算整体坐标系下的节点位移{δ′_I}：
对于第I荷载步，施加荷载{P}/N，{P}表示荷载，N表示荷载步数，I∈[1,N]；根据[K_I]
{δ′_I}＝{P}/N，求第I步的节点位移{δ′_I}；其中I＝1时，结构初始损伤D₁为0，[K₁]＝[K]，[K]
为步骤5中得到的自由曲面结构的整体刚度矩阵；
(6-2)通过坐标转换，将{δ′_I}转换为单元坐标系下的节点位移向量{δ_I}，{δ_I}＝[L]
{δ′_I}；并将{δ_I}分解为平面应力单元对应的位移和弯曲板单元对应的位移和
(6-3)由单元坐标系下的和分别计算平面应力单元和弯曲板单元的应力向量
和并得到第I荷载步的单元总应力向量{σ_I}，具体为：
平面应力单元的应力向量为：
${\left\{{\mathrm{\σ}}_I^p\right\}}^T=\{{\left\{{\mathrm{\σ}}_r^p\right\}}^T,{\left\{{\mathrm{\σ}}_s^p\right\}}^T,{\left\{{\mathrm{\σ}}_t^p\right\}}^T\}$
其中，
[D^p]表示平面应力单元的弹性矩阵，[B^p]表示平面应力单元
的应变矩阵；
弯曲板单元的应力向量为：
${\left\{{\mathrm{\σ}}_I^b\right\}}^T=\{{\left\{{\mathrm{\σ}}_r^b\right\}}^T,{\left\{{\mathrm{\σ}}_s^b\right\}}^T,{\left\{{\mathrm{\σ}}_t^b\right\}}^T\}$
其中，
[D^b]表示弯曲板单元的弹性矩阵，[B^b]表示弯曲板单
元的应变矩阵；
第I荷载步的单元总应力向量为：
{σ_I}^T＝{{σ_Ir}^T，{σ_Is}^T，{σ_It}^T}
其中，{σ_Ir}表示第I荷载步r节点的总应力向量，{σ_Is}表示第I荷载
步s节点的总应力向量，{σ_It}表示第I荷载步t节点的总应力向量，

(6-4)由(6-3)得到的第I荷载步的单元总应力向量{σ_I}，计算前I荷载步的累计应力向
量并由{σ^*}求解前I荷载步累计主应力向量和累计主应变向量具体计算方
法如下：
(6-4-1)计算前I荷载步的累计应力向量
${\left\{{\mathrm{\σ}}_I^{*}\right\}}^T=\{{\left\{{\mathrm{\σ}}_{Ir}^{*}\right\}}^T,{\left\{{\mathrm{\σ}}_{Is}^{*}\right\}}^T,{\left\{{\mathrm{\σ}}_{It}^{*}\right\}}^T\}$
其中，r节点的累计应力向量s节点的累计应力向量
t节点的累计应力向量
(6-4-2)由前I荷载步的累计应力向量求前I荷载步累计主应力向量
${\left\{{\mathrm{\σ}}_M^{*}\right\}}^T=\{{\left\{{\mathrm{\σ}}_{Mr}^{*}\right\}}^T,{\left\{{\mathrm{\σ}}_{Ms}^{*}\right\}}^T,{\left\{{\mathrm{\σ}}_{Mt}^{*}\right\}}^T\}$
其中，r节点的节点主应力记为3×1的向量

将r节点的累计应力向量记为3×1的向量；s节点的节点主应
力记为3×1的向量
将s节点的累计应力向量记为3×1的向量；t节点的
节点主应力记为3×1的向量
将t节点的累计应力向量记为3×1的向量；
(6-4-3)由累计单元主应力求累计单元主应变
${\left\{{\mathrm{\ϵ}}_M^{*}\right\}}^T=\{{\left\{{\mathrm{\ϵ}}_{Mr}^{*}\right\}}^T,{\left\{{\mathrm{\ϵ}}_{Ms}^{*}\right\}}^T,{\left\{{\mathrm{\ϵ}}_{Mt}^{*}\right\}}^T\}$
其中，r节点的节点主应变记为3×1的向量
s节点的节点主应变记为3×1
的向量t
节点的节点主应变记为3×1的向量
E为混凝土的弹性模量；
(6-5)求解第I+1步的损伤因子D_I+1:
$D_{I+1}=1/3\·(D_{I+1}^r+D_{I+1}^s+D_{I+1}^t)$
其中，分别为r、s、t节点的损伤因子；
(6-5-1)r节点的损伤因子计算方法包括以下三种情况：
Ⅰ、若均大于0，单元处于拉伸应力状态，则
$D_{I+1}^r=\left\{\begin{array}{cc}0& {\mathrm{\ϵ}}_{Mr1}^{*}\≤{\mathrm{\ϵ}}_f\\ 1-\frac{{\mathrm{\ϵ}}_f(1-A_t)}{{\mathrm{\ϵ}}_{Mr1}^{*}}-\frac{A_t}{\exp \[B_t({\mathrm{\ϵ}}_{Mr1}^{*}-{\mathrm{\ϵ}}_f)\]}& {\mathrm{\ϵ}}_{Mr1}^{*}\>{\mathrm{\ϵ}}_f\end{array}\right.$
其中，A_t和B_t为拉伸时的经验系数，ε_f为峰值应力的应变；
Ⅱ、若均小于0，单元处于压缩应力状态，则：
$D_{I+1}^r=\left\{\begin{array}{cc}0& \left|{\mathrm{\ϵ}}_{Mr2}^{*}\right|\≤{\mathrm{\ϵ}}_f\\ 1-\frac{{\mathrm{\ϵ}}_f(1-A_c)}{\left|{\mathrm{\ϵ}}_{Mr2}^{*}\right|}-\frac{A_c}{\exp \[B_c\left(\right|{\mathrm{\ϵ}}_{Mr2}^{*}|-{\mathrm{\ϵ}}_f)\]}& \left|{\mathrm{\ϵ}}_{Mr2}^{*}\right|\>{\mathrm{\ϵ}}_f\end{array}\right.$
其中，A_c和B_c为压缩时的经验系数；
Ⅲ、若大于0、小于0，则单元处于复杂应力状态，则：
$D_{I+1}^r={\mathrm{\α}}_T{D}_T+{\mathrm{\α}}_c{D}_c$
$D_T=\left\{\begin{array}{cc}0& {\mathrm{\ϵ}}_{Mr1}^{*}\≤{\mathrm{\ϵ}}_f\\ 1-\frac{{\mathrm{\ϵ}}_f(1-A_t)}{{\mathrm{\ϵ}}_{Mr1}^{*}}-\frac{A_t}{\exp \[B_t({\mathrm{\ϵ}}_{Mr1}^{*}-{\mathrm{\ϵ}}_f)\]}& {\mathrm{\ϵ}}_{Mr1}^{*}\>{\mathrm{\ϵ}}_f\end{array}\right.$
$D_c=\left\{\begin{array}{cc}0& \left|{\mathrm{\ϵ}}_{Mr2}^{*}\right|\≤{\mathrm{\ϵ}}_f\\ 1-\frac{{\mathrm{\ϵ}}_f(1-A_c)}{\left|{\mathrm{\ϵ}}_{Mr2}^{*}\right|}-\frac{A_c}{\exp \[B_c\left(\right|{\mathrm{\ϵ}}_{Mr2}^{*}|-{\mathrm{\ϵ}}_f)\]}& \left|{\mathrm{\ϵ}}_{Mr2}^{*}\right|\>{\mathrm{\ϵ}}_f\end{array}\right.$
其中，α_T和α_c为耦合系数，纯压缩时取α_T＝0，纯拉伸时取α_c＝0，组合情况α_T+α_c＝1；
(6-5-2)采用6-5-1的方法计算得到s、t节点的损伤因子
(6-6)求解第I+1步的单元刚度矩阵和自由曲面结构的总体刚矩阵：
第I+1步的单元刚度矩阵计算方法如下：
$\[k_{I+1}\]=\left[\begin{array}{cccccc}\[k_{I+1}^p\]& & 0& 0& 0& 0\\ & & 0& 0& 0& 0\\ 0& 0& & & & 0\\ 0& 0& & \[k_{I+1}^p\]& & 0\\ 0& 0& & & & 0\\ 0& 0& 0& 0& 0& 0\end{array}\right]$
其中，
通过第I+1步的整体坐标系下的单元刚度矩阵，得到第I+1步的整体坐标系下的自由曲
面结构的总体刚度矩阵[K_I+1]。
(6-7)迭代步消除不平衡力：
在整体坐标下，记第J迭代步的位移误差向量为且J∈[1,M]，M为迭代步数，迭代公
式如下：
第I荷载步第1迭代步的公式为：求解{δ′₁}；
第I荷载步第J迭代步公式为：
$\[k_{I+1}\]\·\left\{\stackrel{\‾}{{\mathrm{\δ}}_J^{\′}}\right\}=\frac{1}{N}\·\left\{p\right\}-\[k_{I+1}\]\·\left(\right\{{\mathrm{\δ}}_I^{\′}\}+\{\stackrel{\‾}{{\mathrm{\δ}}_1^{\′}}\}+\{\stackrel{\‾}{{\mathrm{\δ}}_2^{\′}}\}+...+\{\sqrt{{\mathrm{\δ}}_{J-1}^{\′}}\left\}\right)$
第I荷载步第M迭代步的迭代公式为：
$\[k_{I+1}\]\·\left\{\stackrel{\‾}{{\mathrm{\δ}}_M^{\′}}\right\}=\frac{1}{N}\·\left\{p\right\}-\[k_{I+1}\]\·\left(\right\{{\mathrm{\δ}}_I^{\′}\}+\{\stackrel{\‾}{{\mathrm{\δ}}_1^{\′}}\}+\{\stackrel{\‾}{{\mathrm{\δ}}_2^{\′}}\}+...+\{\stackrel{\‾}{{\mathrm{\δ}}_{M-1}^{\′}}\left\}\right)$
(6-8)考虑Mazars损伤模型的{p}荷载作用下的第I荷载步的节点位移{W_I}为：
$\left\{W_I\right\}=\left\{{\mathrm{\δ}}_I^{\′}\right\}+\left\{\stackrel{\‾}{{\mathrm{\δ}}_1^{\′}}\right\}+\left\{\stackrel{\‾}{{\mathrm{\δ}}_2^{\′}}\right\}+...+\left\{\stackrel{\‾}{{\mathrm{\δ}}_M^{\′}}\right\}$
求得第I荷载步的节点位移{W_I}后，返回6-1求解第I+1荷载步的节点位移{W_I+1}，且I∈
[1,N]；
将N个荷载步的所有节点位移累加，得到考虑Mazars损伤模型的{p}荷载作用下的节点
位移{W^*}为：
{W^*}＝{W₁}+{W₂}+…+{W_N}。
7.根据权利要求6所述的一种考虑损伤的自由曲面形态创建方法，其特征在于，考虑
Mazars损伤的自由曲面控制点高度调整，具体为：
记L为自由曲面结构的第L个优化步，L为正整数；第L个优化步和第L+1个优化步的控制
点高度差Δz^(L)为：
${\mathrm{\Δz}}^{\left(L\right)}={\mathrm{\λ}}^{\left(L\right)}\·\underset{i=1,j=1}{\overset{m\×n}{\Σ}}\frac{\∂C\left(z^{\left(L\right)}\right)}{\∂z_{i,j}}{\mathrm{\Δz}}_{i,j}$
其中，表示第L优化步的应变能关于第L优化步的控制点高度的梯度，
λ^(L)表示第L优化步的步长，z_i,j表示u方向的第i个和v方向第j个控制点的高度坐标，Δz_i,j
表示u方向的第i个和v方向第j个控制点的高度坐标微小增量；
进一步，第L+1优化步的控制点高度如下式所示为：
$z_{i,j}^{(L+1)}=z_{i,j}^{\left(L\right)}-{\mathrm{\Δz}}^{\left(L\right)}=z_{i,j}^{\left(L\right)}-{\mathrm{\λ}}^{\left(L\right)}\underset{i=1,j=1}{\overset{m\×n}{\Σ}}\frac{\∂C\left(z^{\left(L\right)}\right)}{\∂z_{i,j}}{\mathrm{\Δz}}_{i,j}$
其中，为第L优化步的控制点高度。
8.根据权利要求7所述的一种考虑损伤的自由曲面形态创建方法，其特征在于，采用应
变能差分算法求解梯度，具体为：
结构应变能C是关于自由曲面控制点高度z^(L)的连续可微函数，将C(z^(L))进行Taylor展
开：
$C(z^{\left(L\right)}+{\mathrm{\Δz}}^{\left(L\right)})=C\left(z^{\left(L\right)}\right)+\underset{i=1,j=1}{\overset{m\×n}{\Σ}}\frac{\∂C\left(z^{\left(L\right)}\right)}{\∂z_{i,j}}{\mathrm{\Δz}}_{i,j}+o\left({\mathrm{\Δz}}_{i,j}\right)$
$C(z^{\left(L\right)}-{\mathrm{\Δz}}^{\left(L\right)})=C\left(z^{\left(L\right)}\right)-\underset{i=1,j=1}{\overset{m\×n}{\Σ}}\frac{\∂C\left(z^{\left(L\right)}\right)}{\∂z_{i,j}}{\mathrm{\Δz}}_{i,j}+o\left({\mathrm{\Δz}}_{i,j}\right)$
得到梯度的表达式为：
$\underset{i=1,j=1}{\overset{m\×n}{\mathrm{\Σ}}}\frac{\∂C\left(z^{\left(L\right)}\right)}{\∂z_{i,j}}{\mathrm{\Δz}}_{i,j}=\frac{C\left(z^{\left(L\right)}+{\mathrm{\Δz}}^{\left(L\right)}\right)-C\left(z^{\left(L\right)}-{\mathrm{\Δz}}^{\left(L\right)}\right)}{2{\mathrm{\Δz}}^{\left(L\right)}}-\frac{{\left({\mathrm{\Δz}}_i\right)}^2}{6}{C}^{\left(3\right)}\left(\mathrm{\ξ}\right)$
其中，o(Δz_i,j)是高阶无穷小量；且ξ∈(0，Δz_i,j)，C⁽³⁾(ξ)表示应变能关于ξ的三阶导
数。
9.根据权利要求8所述的一种考虑损伤的自由曲面形态创建方法，其特征在于，采用黄
金分割法求解步长，具体为：
设定F表示黄金分割法的第F迭代步，λ_F为第F次迭代求得的步长值，λ为黄金分割法最终
要求的步长值。当程序循环至满足step2的条件时将λ_F的值即为第L优化步所求的步长λ^(L)；
step1：置初始区间[a₁,b₁]和精度要求G>0，计算试探点λ₁和μ₁，计算C(λ₁)和C(μ₁)，令F
＝1
λ₁＝a₁+0.382(b₁-a₁)
μ₁＝a₁+0.618(b₁-a₁)
step2：若b_k-a_k<G，则停止计算；否则，当C(λ_k)>C(μ_k)时，转step3；当C(λ_F)<C(μ_F)时，转
step4；
step3：置a_F+1＝λ_F，b_F+1＝b_F，λ_F+1＝μ_F
μ_F+1＝a_F+1+0.618(b_F+1-a_F+1)
计算应变能C(μ_F+1)，转step5；
step4：置a_F+1＝a_F，b_F+1＝μ_F，μ_F+1＝λ_F
λ_F+1＝a_F+1+0.382(b_F+1-a_F+1)
计算应变能C(λ_F+1)，转step5；
step5：置F＝F+1，返回step2。

一种考虑损伤的自由曲面形态创建方法

技术领域

本发明涉及一种自由曲面形态创建方法，特别涉及一种考虑损伤的自由曲面形态
创建方法。

背景技术

国内外目前自由曲面结构的形态创建方法可以概括为以下两种：

1.基于试验的物理方法，即以工程经验或物理实验为主要手段，有意识地确定合
理结构形状；

2.基于数值的优化方法，即基于计算机技术，以计算机图形学和结构优化等为手
段创建合理结构形态。

在现有的自由曲面形态结构创建方法中，自由曲面的模型都采用线弹性的方法进
行分析和优化。此类方法没有考虑到随着荷载的施加和变形的积累，材料在自由曲面创建
过程中出现的微观损伤。而材料的损伤将直接导致结构刚度下降，因此如果不考虑材料的
削弱对结构设计的影响，会导致结构设计偏于不安全。

发明内容

本发明所要解决的技术问题是提供一种考虑损伤的自由曲面形态创建方法。该方
法算法简单有效，考虑自由曲面创建过程中的材料损伤。

本发明为解决上述技术问题采用以下技术方案：

本发明提供一种考虑损伤的自由曲面形态创建方法，具体方法步骤如下：

步骤1，根据给定的初始自由曲面的一组控制点坐标初始化自由曲面；

步骤2，对步骤1得到的初始化自由曲面进行网格划分，得到三角形薄壳单元；

步骤3，计算步骤2得到的三角形薄壳单元的单元坐标系下的刚度矩阵；

步骤4，计算整体坐标下的单元刚度矩阵；

步骤5，经过总纲集成得到自由曲面结构的整体刚度矩阵；

步骤6，加入Mazars损伤模型，对结构施加荷载，分N个荷载步进行加载，每个荷载
步内存在M个迭代步，计算自由曲面壳体位移；

步骤7，采用应变能作为自由曲面结构创建的优化目标，节点坐标作为自由曲面结
构创建的优化变量，通过节点荷载和节点位移计算应变能；

步骤8，对考虑Mazars损伤的自由曲面的控制点进行高度调整；

步骤9，设置优化精度ε^*，若当前优化步的应变能对当前优化步的控制点高度的差
分的模小于优化精度，则考虑损伤的自由曲面优化结束，自由曲面创建完成；否则继续迭代
进行自由曲面结构形态优化。

作为本发明的进一步优化方案，步骤1中的自由曲面是NURBS曲面。

作为本发明的进一步优化方案，自由曲面是一张在u方向p次、v方向q次的NURBS曲
面，具有如下形式的双变量分段有理矢值函数：

其中，n为u方向的控制点个数，且i∈[1,n]；m为v方向的控制点个数，且j∈[1,m]；
P_i,j是u方向编号为i且v方向编号为j的控制点坐标；w_i,j是u方向编号为i且v方向编号为j的
控制点的全因子；N_i,p(u)和N_j,q(v)是分别定义在矢量U和V上的样条基函数，u_i为u方向编号
为i的节点矢量且u_i∈U，v_j为v方向编号为j
的节点矢量且

N_i,p(u)和N_j,q(v)的表达式如下所示：

作为本发明的进一步优化方案，步骤3中计算三角形薄壳单元的单元坐标系下的
刚度矩阵，具体为：

三角形薄壳单元的刚度矩阵由平面应力状态和弯曲应力状态的刚度矩阵得到，三
角形薄壳单元的三个节点分别记为r、s、t；

在三角形平面应力单元中，节点力和节点位移的关系如下：

其中，[k^p]表示单元坐标系下三角形平面应力单元的刚度矩阵；为单元坐标系
下平面应力单元的r节点的节点力，且由u方向的节点力U_r和v方向的节点力V_r组成，
为单元坐标系下平面应力单元的s节点的节点力，且由u方向的节点力U_s
和v方向的节点力V_s组成，为单元坐标系下平面应力单元的t节点的节点
力，且由u方向的节点力U_t和v方向的节点力V_t组成，表示单元坐标系下三
角形平面应力单元的r节点的位移，且由u方向的位移u_r和v方向的位移v_r组成，
表示单元坐标系下三角形平面应力单元的s节点的位移，且由u方向的位
移u_s和v方向的位移v_s组成，表示单元坐标系下三角形平面应力单元的t
节点的位移，且由u方向的位移u_t和v方向的位移v_t组成，

在三角形弯曲应力单元中，节点力和节点位移的关系如下：

其中，[k^b]表示单元坐标系下三角形弯曲板单元的刚度矩阵；为单元坐标系下
弯曲板单元r节点的节点力，且由垂直于u、v平面方向的节点力W_r、绕u轴旋转的节点力M_ur及
绕v轴旋转的节点力M_vr组成，为单元坐标系下弯曲板单元s节点的节点
力，且由垂直于u、v平面方向的节点力W_s、绕u轴旋转的节点力M_us及绕v轴旋转的节点力M_vs组
成，为单元坐标系下弯曲板单元t节点的节点力，且由垂直于u、v平面方
向的节点力W_t、绕u轴旋转的节点力M_ut及绕v轴旋转的节点力M_vt组成，表
示单元坐标系下三角形弯曲板单元的r节点的位移，垂直于u、v平面方向的线位移w_r、绕u轴
旋转角θ_ur及绕v轴旋转角θ_vr组成，表示单元坐标系下三角形弯曲板单元
的s节点的位移，垂直于u、v平面方向的线位移w_s、绕u轴旋转角θ_us及绕v轴旋转角θ_vs组成，
表示单元坐标系下三角形弯曲板单元的t节点的位移，垂直于u、v平面方
向的线位移w_t、绕u轴旋转角θ_ut及绕v轴旋转角θ_vt组成，

三角形薄壳单元的单元坐标系下的刚度矩阵为：

作为本发明的进一步优化方案，步骤4中整体坐标下的单元刚度矩阵[k′]为：

[k′]＝[L]^-1[k][L]

其中，[L]为坐标转换矩阵，{δ′_r}
＝[u′_rv′_rw′_rθ′_xrθ′_yrθ′_zr]^T，θ_r表示单元坐标系下r节点绕u、v面法线的转角，u′_r、v′_r、w′_r、
θ′_xr、θ′_yr、θ′_zr分别表示整体坐标系下的三角形单元的r节点在整体坐标系下x、y、z方向的
位移以及绕x、y、z方向的旋转角；{F_r}＝[λ]{F′_r}，M_θr表示单元坐标系下r
节点的虚拟弯矩，{F′_r}＝[U′_r V′_r W′_rM′_θxr M′_θyr M′_θzr]^T，U′_r、V′_r、W′_r、M′_θxr、M′_θyr、M′_θzr分
别表示整体坐标系下的三角形单元的r节点在整体坐标系下x、y、z方向的节点力以及绕x、
y、z方向的旋转角。

作为本发明的进一步优化方案，步骤6中计算自由曲面壳体位移的具体步骤如下：

(6-1)计算整体坐标系下的节点位移{δ′_I}：

对于第I荷载步，施加荷载{P}/N，{P}表示荷载，N表示荷载步数，I∈[1,N]；根据
[K_I]{δ′_I}＝{P}/N，求第I步的节点位移{δ′_I}；其中I＝1时，结构初始损伤D₁为0，[K₁]＝[K]，
[K]为步骤5中得到的自由曲面结构的整体刚度矩阵；

(6-2)通过坐标转换，将{δ′_I}转换为单元坐标系下的节点位移向量{δ_I}，{δ_I}＝
[L]{δ′_I}；并将{δ_I}分解为平面应力单元对应的位移和弯曲板单元对应的位移和

(6-3)由单元坐标系下的和分别计算平面应力单元和弯曲板单元的应力
向量和并得到第I荷载步的单元总应力向量{σ_I}，具体为：

平面应力单元的应力向量为：

其中，
[D^p]表示平面应力单元的弹性矩阵，[B^p]表示平面应力单元
的应变矩阵；

弯曲板单元的应力向量为：

其中，
[D^b]表示弯曲板单元的弹性矩阵，[B^b]表示弯曲板单
元的应变矩阵；

第I荷载步的单元总应力向量为：

{σ_I}^T＝{{σ_Ir}^T，{σ_Is}^T，{σ_It}^T}

其中，{σ_Ir}表示第I荷载步r节点的总应力向量，表示第I
荷载步s节点的总应力向量，{σ_It}表示第I荷载步t节点的总应力向量，

(6-4)由(6-3)得到的第I荷载步的单元总应力向量{σ_I}，计算前I荷载步的累计应
力向量并由{σ^*}求解前I荷载步累计主应力向量和累计主应变向量具体计
算方法如下：

(6-4-1)计算前I荷载步的累计应力向量

其中，r节点的累计应力向量s节点的累计应力向
量t节点的累计应力向量

(6-4-2)由前I荷载步的累计应力向量求前I荷载步累计主应力向量

其中，r节点的节点主应力记为3×1的向量

将r节点的累计应力向量记为3×1的向
量；s节点的节点主应力记为3×1的向量

将s节点的累计应力向量记为3×1
的向量；t节点的节点主应力记为3×1的向量

将t节点的累计应力向量记为3×1的向量；

(6-4-3)由累计单元主应力求累计单元主应变

其中，r节点的节点主应变记为3×1的向量
s节点的节点主应变记为3×1
的向量t
节点的节点主应变记为3×1的向量
E为混凝土的弹性模量；

(6-5)求解第I+1步的损伤因子D_I+1:

其中，分别为r、s、t节点的损伤因子；

(6-5-1)r节点的损伤因子计算方法包括以下三种情况：

Ⅰ、若均大于0，单元处于拉伸应力状态，则

其中，A_t和B_t为拉伸时的经验系数，ε_f为峰值应力的应变；

Ⅱ、若均小于0，单元处于压缩应力状态，则：

其中，A_c和B_c为压缩时的经验系数；

Ⅲ、若大于0、小于0，则单元处于复杂应力状态，则：

其中，α_T和α_c为耦合系数，纯压缩时取α_T＝0，纯拉伸时取α_c＝0，组合情况α_T+α_c＝1；

(6-5-2)采用6-5-1的方法计算得到s、t节点的损伤因子

(6-6)求解第I+1步的单元刚度矩阵和自由曲面结构的总体刚矩阵：

第I+1步的单元刚度矩阵计算方法如下：

其中，

通过第I+1步的整体坐标系下的单元刚度矩阵，得到第I+1步的整体坐标系下的自
由曲面结构的总体刚度矩阵[K_I+1]。

(6-7)迭代步消除不平衡力：

在整体坐标下，记第J迭代步的位移误差向量为且J∈[1,M]，M为迭代步数，迭
代公式如下：

第I荷载步第1迭代步的公式为：求解{δ′₁}；

第I荷载步第J迭代步公式为：

第I荷载步第M迭代步的迭代公式为：

(6-8)考虑Mazars损伤模型的{p}荷载作用下的第I荷载步的节点位移{W_I}为：

求得第I荷载步的节点位移{W_I}后，返回6-1求解第I+1荷载步的节点位移{W_I+1}，
且I∈[1,N]；

将N个荷载步的所有节点位移累加，得到考虑Mazars损伤模型的{p}荷载作用下的
节点位移{W^*}为：

作为本发明的进一步优化方案，考虑Mazars损伤的自由曲面控制点高度调整，具
体为：

记L为自由曲面结构的第L个优化步，L为正整数；第L个优化步和第L+1个优化步的
控制点高度差Δz^(L)为：

其中，表示第L优化步的应变能关于第L优化步的控制点高度的
梯度，λ^(L)表示第L优化步的步长，z_i,j表示u方向的第i个和v方向第j个控制点的高度坐标，
Δz_i,j表示u方向的第i个和v方向第j个控制点的高度坐标微小增量；

进一步，第L+1优化步的控制点高度如下式所示为：

其中，为第L优化步的控制点高度。

作为本发明的进一步优化方案，采用应变能差分算法求解梯度，具体为：

结构应变能C是关于自由曲面控制点高度z^(L)的连续可微函数，将C(z^(L))进行
Taylor展开：

得到梯度的表达式为：

其中，o(Δz_i,j)是高阶无穷小量；且ξ∈(0，Δz_i,j)，C⁽³⁾(ξ)表示应变能关于ξ的三
阶导数。

作为本发明的进一步优化方案，采用黄金分割法求解步长，具体为：

设定F表示黄金分割法的第F迭代步，λ_F为第F次迭代求得的步长值，λ为黄金分割
法最终要求的步长值。当程序循环至满足step2的条件时将λ_F的值即为第L优化步所求的步
长λ^(L)；

step1：置初始区间[a₁,b₁]和精度要求G>0，计算试探点λ₁和μ₁，计算C(λ₁)和C(μ₁)，
令F＝1

λ₁＝a₁+0.382(b₁-a₁)

μ₁＝a₁+0.618(b₁-a₁)

step2：若b_k-a_k<G，则停止计算；否则，当C(λ_k)>C(μ_k)时，转step3；当C(λ_F)<C(μ_F)
时，转step4；

step3：置a_F+1＝λ_F，b_F+1＝b_F，λ_F+1＝μ_F

μ_F+1＝a_F+1+0.618(b_F+1-a_F+1)

计算应变能C(μ_F+1)，转step5；

step4：置a_F+1＝a_F，b_F+1＝μ_F，μ_F+1＝λ_F

λ_F+1＝a_F+1+0.382(b_F+1-a_F+1)

计算应变能C(λ_F+1)，转step5；

step5：置F＝F+1，返回step2。

本发明采用以上技术方案与现有技术相比，具有以下技术效果：本发明一种考虑
损伤的自由曲面形态创建方法，考虑到随着荷载的施加和变形的积累，材料在自由曲面创
建过程中出现的微观损伤：

1.考虑损伤的自由曲面创建能够保证自由曲面的受力性能和整体刚度处于局部
最优；

2.现有创建方法中不考虑自由曲面的损伤会造成的自由曲面创建结果偏于不安
全，通过以上技术方案可以改进这一缺点。

附图说明

图1是本发明的方法流程图。

图2是本发明实施例的自由结构的正视图和俯视图，其中，(a)为正视图，(b)为俯
视图。

图3是本发明实施例的考虑mazars损伤模型的优化结果和不考虑mazars损伤模型
的优化结果，其中，(a)为考虑mazars损伤模型的优化结果，(b)为不考虑mazars损伤模型的
优化结果。

图4是本发明实施例的考虑损伤的结构损伤和不考虑损伤的结构损伤示意图，其
中，(a)为考虑损伤的结构损伤，(b)为不考虑损伤的结构损伤。

图5是本发明实施例的考虑损伤的结构位移和不考虑损伤的结构位移，其中，(a)
为考虑损伤的结构位移，(b)为不考虑损伤的结构位移。

图6是本发明实施例的考虑损伤的结构Mises应力和不考虑损伤的结构Mises应
力，其中，(a)为考虑损伤的结构Mises应力，(b)为不考虑损伤的结构Mises应力。

具体实施方式

下面结合附图以及具体实施例，对本发明的技术方案做进一步的详细说明：

本发明提供一种考虑损伤的自由曲面形态创建方法，如图1所示，具体方法步骤如
下：

步骤1，根据给定的初始自由曲面的一组控制点坐标{P_i,j}初始化自由曲面，其中
{P_i,j}表示由控制点坐标P_i,j组成的集合。

本发明中自由曲面为是一张在u方向p次、v方向q次的NURBS曲面，具有如下形式的
双变量分段有理矢值函数：

其中，n为u方向的控制点个数，且i∈[1,n]；m为v方向的控制点个数，且j∈[1,m]；
P_i,j是u方向编号为i且v方向编号为j的控制点坐标；w_i,j是u方向编号为i且v方向编号为j的
控制点的全因子；N_i,p(u)和N_j,q(v)是分别定义在矢量U和V上的样条基函数，u_i为u方向编号
为i的节点矢量且u_i∈U，v_j为v方向编号为j
的节点矢量且

N_i,p(u)和N_j,q(v)的表达式如下所示：

步骤2，对步骤1得到的初始化自由曲面进行网格划分，得到三角形薄壳单元。

步骤3，计算步骤2得到的三角形薄壳单元的单元坐标系下的刚度矩阵。

弹性薄壳的应力状态由平面应力状态和弯曲应力状态组合而成。三角形薄壳单元
的刚度矩阵也由这两种应力状态的刚度矩阵加以组合得到。

对于一个三角形薄壳单元，存在三个节点分别记为r，s，t。

在三角形平面应力单元中，节点力和节点位移的关系如下：

其中，[k^p]表示单元坐标系下三角形平面应力单元的刚度矩阵；为单元坐标
系下平面应力单元的r节点的节点力，且由u方向的节点力U_r和v方向的节点力V_r组成，
为单元坐标系下平面应力单元的s节点的节点力，且由u方向的节点力U_s
和v方向的节点力V_s组成，为单元坐标系下平面应力单元的t节点的节点
力，且由u方向的节点力U_t和v方向的节点力V_t组成，表示单元坐标系下三
角形平面应力单元的r节点的位移，且由u方向的位移u_r和v方向的位移v_r组成，
表示单元坐标系下三角形平面应力单元的s节点的位移，且由u方向的位
移u_s和v方向的位移v_s组成，表示单元坐标系下三角形平面应力单元的t
节点的位移，且由u方向的位移u_t和v方向的位移v_t组成，

在三角形弯曲应力单元中，节点力和节点位移的关系如下：

其中，[k^b]表示单元坐标系下三角形弯曲板单元的刚度矩阵；为单元坐标系
下弯曲板单元r节点的节点力，且由垂直于u、v平面方向的节点力W_r、绕u轴旋转的节点力M_ur
及绕v轴旋转的节点力M_vr组成，为单元坐标系下弯曲板单元s节点的节
点力，且由垂直于u、v平面方向的节点力W_s、绕u轴旋转的节点力M_us及绕v轴旋转的节点力M_vs
组成，为单元坐标系下弯曲板单元t节点的节点力，且由垂直于u、v平面
方向的节点力W_t、绕u轴旋转的节点力M_ut及绕v轴旋转的节点力M_vt组成，
表示单元坐标系下三角形弯曲板单元的r节点的位移，垂直于u、v平面方向的线位移w_r、绕u
轴旋转角θ_ur及绕v轴旋转角θ_vr组成，表示单元坐标系下三角形弯曲板单
元的s节点的位移，垂直于u、v平面方向的线位移w_s、绕u轴旋转角θ_us及绕v轴旋转角θ_vs组成，
表示单元坐标系下三角形弯曲板单元的t节点的位移，垂直于u、v平面方
向的线位移w_t、绕u轴旋转角θ_ut及绕v轴旋转角θ_vt组成，

把平面应力和弯曲应力加以组合。其中θ_r表示绕u,v面法线的转角，虽然转角θ_r不
影响单元应力，为了便于以后把局部坐标系的刚度矩阵转化为整体坐标系下的刚度矩阵并
集合成整体刚度矩阵，把θ_r也包括在节点位移中，并在节点力中相应的包括一个虚拟弯矩
M_θr。

由于平面应力状态下的节点力和弯曲应力状态下的节点位移互不影响；弯曲应力
状态下的节点力和平面应力状态下的节点位移也互不影响。所以，组合应力状态下的单元
坐标系下的刚度矩阵[k]可以写成如下形式：

步骤4，计算整体坐标下的单元刚度矩阵。

整体坐标系下的单元刚度矩阵[k′]由局部坐标系下的单元刚度矩阵[k]乘以坐标
转换矩阵[L]得到。

整体坐标下的单元刚度矩阵[k′]为：

[k′]＝[L]^-1[k][L]

其中，[L]为坐标转换矩阵，{δ′_r}
＝[u′_r v′_r w′_r θ′_xr θ′_yr θ′_zr]^T，θ_r表示单元坐标系下r节点绕u、v面法线的转角，u′_r、v′_r、
w′_r、θ′_xr、θ′_yr、θ′_zr分别表示整体坐标系下的三角形单元的r节点在整体坐标系下x、y、z方
向的位移以及绕x、y、z方向的旋转角；{F_r}＝[λ]{F′_r}，M_θr表示单元坐标系
下r节点的虚拟弯矩，{F′_r}＝[U′_r V′_r W′_r M′_θxr M′_θyr M′_θzr]^T，U′_r、V′_r、W′_r、M′_θxr、M′_θyr、
M′_θzr分别表示整体坐标系下的三角形单元的r节点在整体坐标系下x、y、z方向的节点力以
及绕x、y、z方向的旋转角。

把整个三角形单元(含r，s，t三个节点)在单元中的节点位移和节点力分别用{δ_i}^e
和{F_i}^e表示，在整体坐标中的节点位移和节点力分别用{δ′_i}^e和{F′_i}^e表示可得到：

{δ_i}^e＝[L]{δ′_i}^e

{F_i}^e＝[L]{F′_i}^e

由上式子可得：

{F′_i}^e＝[L]^-1{F_i}^e＝[L]^-1[k]{δ_i}^e＝[L]^-1[k][L]{δ′_i}^e

因此，求得整体坐标下的刚度矩阵[k′]为：

[k′]＝[L]^-1[k][L]

步骤5，经过总纲集成得到自由曲面结构的整体刚度矩阵。

步骤6，加入Mazars损伤模型，对结构施加荷载{P}，分N个荷载步进行加载，每个荷
载步内存在M个迭代步，采用以下方法计算自由曲面壳体位移{W^*}。

每次荷载步加载{P}/N，记I为荷载步的第I步，且I∈[1,N]。记第I荷载步内整体坐
标系下的刚度矩阵为[K_I]，局部坐标系下的单元刚度矩阵记为[k_I]，与之对应的平面应力状
态下的刚度矩阵记为弯曲应力状态下的刚度矩阵记为记第I荷载步的损伤因子
为D_I。整体坐标下的节点位移向量为{δ′_I}，单元坐标系的节点位移向量为{δ_I}，其中单元坐
标系下的节点位移分平面应力单元对应的位移和弯曲板单元对应的位移分别记为和
同样的单元坐标系下的平面应力单元和弯曲板单元的应力向量分别记为和
第I荷载步的单元总应力向量记为{σ_I}。前I荷载步的累计单元应力记为

(6-1)计算整体坐标系下的节点位移{δ′_I}：

对于第I荷载步，施加荷载{P}/N，{P}表示荷载，N表示荷载步数，I∈[1,N]；根据
[K_I]{δ′_I}＝{P}/N，求第I步的节点位移{δ′_I}。特别的，对于第一荷载步，即I＝1时，认为结
构初始损伤D₁为0，[K₁]＝[K]。

(6-2)通过坐标转换，将{δ′_I}转换为单元坐标系下的节点位移向量{δ_I}，{δ_I}＝
[L]{δ′_I}；并将{δ_I}分解为平面应力单元对应的位移和弯曲板单元对应的位移和

(6-3)由单元坐标系下的和分别计算平面应力单元和弯曲板单元的应力
向量和并得到第I荷载步的单元总应力向量{σ_I}，具体为：

平面应力单元的应力向量为：

其中，
[D^p]表示平面应力单元的弹性矩阵，[B^p]表示平面应力单元
的应变矩阵；

弯曲板单元的应力向量为：

其中，
[D^b]表示弯曲板单元的弹性矩阵，[B^b]表示弯曲板单
元的应变矩阵；

第I荷载步的单元总应力向量为：

{σ_I}^T＝{{σ_Ir}^T，{σ_Is}^T，{σ_It}^T}

其中，{σ_Ir}表示第I荷载步r节点的总应力向量，表示第I
荷载步s节点的总应力向量，表示第I荷载步t节点的总应力向
量，

(6-4)由(6-3)得到的第I荷载步的单元总应力向量{σ_I}，计算前I荷载步的累计应
力向量并由{σ^*}求解前I荷载步累计主应力向量和累计主应变向量具体计
算方法如下：

(6-4-1)计算前I荷载步的累计应力向量

其中，r节点的累计应力向量s节点的累计应力向
量t节点的累计应力向量

(6-4-2)由前I荷载步的累计应力向量求前I荷载步累计主应力向量

其中，r节点的节点主应力记为3×1的向量

将r节点的累计应力向量记为3×1的向
量；s节点的节点主应力记为3×1的向量

将s节点的累计应力向量记为3×1的向量；t
节点的节点主应力记为3×1的向量

将t节点的累计应力向量记为3×1的向量；

(6-4-3)由累计单元主应力求累计单元主应变

其中，r节点的节点主应变记为3×1的向量
s节点的节点主应变记为3×1的向
量t
节点的节点主应变记为3×1的向量
E为混凝土的弹性模量；

(6-5)求解第I+1步的损伤因子D_I+1。

记第I+1荷载步节点r的损伤因子为(同理可得s，t节点的损伤因子
其中，其计算公式如下：

Mazars采用下式来描述峰值应力以后的软化曲线。其中的ε是指r节点主应变，也
就是(6-4)步骤中求的累计节点主应变中的

情况1：拉伸情况下，峰值应力之后其应力-应变曲线可以用下式来描述，式中A_t和
B_t为拉伸时的经验系数，ε_f为峰值应力的应变。对于混凝土材料取0.5<A_t<1.0，10⁴<B_t<10⁵。

由上式可以得到拉伸状态下损伤演化方程，则

情况2：压缩情况下，峰值应力之后其应力-应变曲线可以用下式来描述，式中A_c和
B_c为压缩时的经验系数，ε_f为峰值应力的应变。对于混凝土材料取1.0<A_c<1.5，1000<B_c<
2000。

由上式可以得到单轴压缩损伤演化方程，则

情况3：复杂应力状态情况下，拉、压损伤之间存在耦合效应。

Mazars建议用下时描述单元总损伤

式中：α_T和α_c为耦合系数，纯压缩时取α_T＝0，纯拉伸时取α_c＝0，组合情况α_T+α_c＝1。

当均大于0时，则单元处于拉伸应力状态，取中的主拉应变
绝对值较大值代入以下情况1中的ε值。

当均小于0时，则单元处于压缩应力状态，取中的主压应变
绝对值较大值代入以下情况2中的ε值。

当均大于0，小于0时，则单元处于复杂应力状态，取代入D_T公式中的
ε，取代入D_c公式中的ε值。

(6-5-2)采用6-5-1的方法计算得到s、t节点的损伤因子

(6-5-3)事实上当单元划分较细时，三角形单元三个节点(r，s，t)的损伤因子相
近。取第I荷载三角形壳单元的三个节点损伤因子的均值，作为该单元第I+1步的损伤因子
D_I+1，即：

其中，分别为r、s、t节点的损伤因子。

(6-6)求解第I+1步的单元刚度矩阵和自由曲面结构的总体刚矩阵：

第I+1步的单元刚度矩阵计算方法如下：

其中，

第I+1步的整体坐标系的单元刚度矩阵形成方法同上述(4)-(5)步骤所述，由此，
通过第I+1步的整体坐标系下的单元刚度矩阵，得到第I+1步的整体坐标系下的自由曲面结
构的总体刚度矩阵[K_I+1]。

(6-7)迭代步消除不平衡力：

在整体坐标下，记第J迭代步的位移误差向量为且J∈[1,M]，M为迭代步数，迭
代公式如下：

第I荷载步第1迭代步的公式为：求解{δ′₁}；

第I荷载步第J迭代步公式为：

第I荷载步第M迭代步的迭代公式为：

(6-8)考虑Mazars损伤模型的{p}荷载作用下的第I荷载步的节点位移{W_I}为：

求得第I荷载步的节点位移{W_I}后，返回6-1求解第I+1荷载步的节点位移{W_I+1}，
且I∈[1,N]；

将N个荷载步的所有节点位移累加，得到考虑Mazars损伤模型的{p}荷载作用下的
节点位移{W^*}为：

{W^*}＝{W₁}+{W₂}+…+{W_N}。

步骤7，采用应变能作为自由曲面结构创建的优化目标，节点坐标作为自由曲面结
构创建的优化变量，通过节点荷载和节点位移计算应变能。

结构的平衡方程为：

[K]·{W^*}＝{p}

式中[K]为结构的整体刚度，{W^*}为结构节点位移向量，{p}为结构荷载向量。结构
的应变能表达式为：

根据应变能的变化调整自由曲面形状，以实现应变能的最小化，改善曲面结构的
整体刚度。若将应变能看成曲面坐标的函数C(z)，合理的自由曲面结构形态的问题可用下
式表达：

步骤8，对考虑Mazars损伤的自由曲面的控制点进行高度调整。

记L为自由曲面结构的第L个优化步，且L属于正整数集。Δz^(L)表示第L个优化步和
第L+1个优化步的控制点高度差，其求解方法如下所示。其中，表示第L优
化步的应变能关于第L优化步的控制点高度的梯度。λ^(L)表示第L优化步的步长，z_i,j表示u方
向的第i和v方向第j个控制点的高度坐标，Δz_i,j表示该控制点的高度坐标微小增量。

由上式得到控制点高度差Δz^(L)后，可以求得第L+1优化步的控制点高度如
下式所示：

(8-1)应变能差分算法求解梯度：

结构应变能C是自由曲面控制点高度z^(L)的函数。假定应变能C是关于控制点高度z
^(L)的连续可微函数，将C(z^(L))进行Taylor展开：

可以得到梯度的表达式如下所示，其中m和n分别为u和v方向的
控制点个数：

式中：o(Δz_i,j)是高阶无穷小量；且ξ∈(0，Δz_i,j)，C⁽³⁾(ξ)表示应变能关于ξ的三
阶导数。

(8-2)采用黄金分割法求解步长：

采用黄金分割法，即通过取试探点使包含极小值点的区间不断缩减其计算步骤如
下，其中F表示黄金分割法的第F迭代步，λ_F为第F次迭代求得的步长值，λ为黄金分割法最终
要求的步长值。当程序循环至满足step2的条件时将λ_F的值赋给λ，该λ值即为第L优化步所
求的步长λ^(L)：

step1：置初始区间[a₁,b₁]和精度要求G>0，计算试探点λ₁和μ₁，计算C(λ₁)和C(μ₁)，
令F＝1

λ₁＝a₁+0.382(b₁-a₁)

μ₁＝a₁+0.618(b₁-a₁)

step2:若b_k-a_k<G，则停止计算。否则，当C(λ_k)>C(μ_k)时，转step3；当C(λ_F)<C(μ_F)
时，转step4。

step3：置a_F+1＝λ_F，b_F+1＝b_F，λ_F+1＝μ_F

μ_F+1＝a_F+1+0.618(b_F+1-a_F+1)

计算应变能C(μ_F+1)，转step5。

step4：置a_F+1＝a_F，b_F+1＝μ_F，μ_F+1＝λ_F

λ_F+1＝a_F+1+0.382(b_F+1-a_F+1)

计算应变能C(λ_F+1)，转step5。

step5：置F＝F+1，返回step2。

步骤9，设置优化精度ε^*，若当第L优化步的应变能C(z^(L))对第L优化步的控制点高
度z^(L)的差分的模小于优化精度ε^*，则考虑损伤的自由曲面优化结束，自由曲面创建完成；
否则将第L+1步的控制点高度z^(L+1)带入第I+1优化步的控制点坐标，回到(1)，继续进行自由
曲面结构形态优化。

具体实施例

设一混凝土自由曲面壳结构，纵向最大跨度L为20m，横向最大跨度B为16m。壳体混
凝土弹性模量E＝3×10¹⁰N·m^-2，泊松比μ＝0.2。混凝土壳厚度为t＝0.1m，壳体受竖向均布
荷载作用，两纵边简支约束。结构初始曲面形状如图2所示，图2中(a)和(b)分别为该自由结
构的正视图和俯视图。

自由曲面创建的控制点位置如图2中(b)的空心圆点所示。优化目标为结构应变能
最小。自由曲面结构取均布荷载q＝120KN/m²，对该自由曲面结构进行优化。考虑mazars损
伤模型的优化结果和不考虑mazars损伤模型的优化结果如图3所示，其中，(a)和(b)分别为
考虑损伤的优化结果和不考虑损伤的优化结果。

表1显示了考虑损伤的自由曲面创建结果的结构最大损伤、位移和Mises应力以及
不考虑损伤的自由曲面创建结果的结构最大损伤、位移和Mises应力。

表1自由曲面结构的最大损伤、位移、mises应力的比较

图4至图6分别为考虑损伤以及不考虑损伤的自由曲面优化结果的损伤、位移和
Mises应力对比。

以上所述，仅为本发明中的具体实施方式，但本发明的保护范围并不局限于此，任
何熟悉该技术的人在本发明所揭露的技术范围内，可理解想到的变换或替换，都应涵盖在
本发明的包含范围之内，因此，本发明的保护范围应该以权利要求书的保护范围为准。