CN110096756A - 一种考虑荷载不确定性的自由曲面结构形态创建方法 - Google Patents

Abstract

本发明公开了一种考虑荷载不确定性的自由曲面结构形态创建方法，该方法算法简单有效，考虑在形态创建过程中对荷载不确定性的忽略会导致自由曲面结构的鲁棒性和可靠性无法得到保证的问题。本发明不确定性优化后的自由曲面结构的刚度和鲁棒性都远优于初始结构，具有更好的受力性能；当荷载的大小不变，方向不确定时，不确定性优化结构的鲁棒性和可靠性要远优于确定性优化结构的鲁棒性和可靠性。

Description

一种考虑荷载不确定性的自由曲面结构形态创建方法

技术领域

本发明涉及一种考虑荷载不确定性的自由曲面结构形态创建方法。

背景技术

近年来，自由曲面结构因其新颖独特的造型而得到日益广泛的应用。其设计的核心问题是形态创建，即实现曲面多样性与受力合理性的有机结合。目前主要采用的形态创建方法是基于优化思想的数值法，该方法将自由曲面的几何建模技术与结构优化算法相结合，以曲面形状的几何参数作为优化变量，以结构受力合理作为优化目标，采用合适的优化算法进行求解，获得最优自由曲面形状。

但现有的自由曲面结构的形态创建方法都是建立荷载为确定性量的假设基础上的。但在实际工程中，荷载(大小、方向)不可避免的存在不确定性。而与传统空间结构相比，自由曲面结构几何造型更复杂、设计和施工难度更大，这就使得其作用荷载具有更强的不确定性。在形态创建过程中对荷载不确定性的忽略会导致自由曲面结构的鲁棒性和可靠性无法得到保证。

发明内容

本发明的目的在于解决上述问题，提供一种考虑荷载不确定性的自由曲面结构形态创建方法。该方法算法简单有效，考虑自由曲面创建过程中的考虑荷载不确定性的影响。

为实现上述目的，本发明采用技术方案如下：

步骤1，采用非均匀有理B样条技术NURBS根据给定的初始控制点坐标建立自由曲面初始几何模型；

步骤2，采用三角形壳单元对步骤1得到的自由曲面初始几何模型进行网格划分，得到混凝土三角形壳单元；

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

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

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

步骤6，利用不确定性参数、不确定性荷载大小及方向的期望值和标准差，建立不确定性荷载的概率模型，对自由曲面结构施加不确定性荷载大小和/或方向的期望值；

步骤7，对有限元平衡方程进行求解，通过结点荷载、结点位移以及不确定性荷载概率模型中的不确定性参数计算出应变能标准差作为形态创建优化目标，将结点最大位移作为约束条件，形成有约束的优化问题；

步骤8，采用广义拉格朗日乘子法将步骤7中的有约束优化问题转换为无约束问题，构成增广目标函数；

步骤9，采用梯度法确定优化方向和步长来更新自由曲面控制点坐标，获得新的自由曲面形状；

步骤10，设置优化精度ε^*，若当前优化步的增广函数对当前优化步的控制点高度的差分的模小于优化精度ε^*，则考虑确定性荷载的自由曲面优化结束，自由曲面创建完成；否则继续迭代至步骤2进行自由曲面结构形态优化。

作为本发明的进一步优化方案，步骤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方向的结点矢量且u_i∈U,v_i为v方向的结点矢量且v_i∈V，

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

(4-1)混凝土单元采用三角形薄壳单元模拟，单元含3个结点，每个结点含6个自由度；

(4-2)混凝土单元采用三角形薄壳单元模拟，三角形薄壳单元的刚度矩阵由平面应力状态和弯曲应力状态的刚度矩阵得到，三角形薄壳单元的三个结点分别记为r、s、t；

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

其中，[k^p]表示局部坐标系下混凝土三角形平面应力单元的刚度矩阵；分别为局部坐标系下混凝土平面应力单元的r、s、t结点的结点力，且由u方向的结点力F_u和v方向的结点力F_v组成，分别表示局部坐标系下混凝土三角形平面应力单元的r、s、t结点的位移，且由u方向的位移u_c和v方向的位移v_c组成；

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

其中，[k^b]表示局部坐标系下混凝土三角形弯曲应力单元的刚度矩阵；分别为局部坐标系下混凝土弯曲应力单元r、s、t结点的结点力，且由垂直于u、v平面w方向的结点力F_w、绕u轴旋转的结点力M_u及绕v轴旋转的结点力M_v组成；分别表示局部坐标系下混凝土三角形弯曲应力单元的r、s、t结点的位移，由垂直于u、v平面w 方向的线位移w_c，绕u轴旋转角θ_u及绕v轴旋转角θ_v组成；

混凝土三角形薄壳单元的局部坐标系下的刚度矩阵[k_c]为：

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

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

其中，[L]为坐标转换矩阵，[k]表示局部坐标系下的单元刚度矩阵，L₁表示局部坐标系各坐标轴与整体坐标系坐标轴之间的方向余弦。

作为本发明的进一步优化方案，步骤5的具体步骤为，由整体坐标系下的单元刚度矩阵 [k′]拼装得到整体坐标系下的整体刚度矩阵K。

作为本发明的进一步优化方案，步骤6中利用不确定性参数、不确定性荷载大小及方向的期望值和标准差，建立不确定性荷载的概率模型，对结构施加不确定性荷载大小和/或方向的期望值，具体步骤如下

(7-1)采用概率模型描述不确定性荷载UF：

假定不确定性荷载UF₁、UF₂……UF_NC均为正态分布随机变量，且相互独立，记为

其中，UF_β表示第β个不确定性荷载，N(·)表示正态分布符号，σ_β分别表示第β个不确定性荷载大小方向的期望值和标准差，NC表示不确定性荷载的个数；

(7-2)在结构上施加NC个荷载

作为本发明的进一步优化方案，步骤7对有限元平衡方程进行求解，通过结点荷载、结点位移以及概率模型中的不确定性参数计算出应变能标准差作为形态创建优化目标，将结点最大位移作为约束条件，具体步骤如下

(8-1)求解有限元平衡方程

根据有限元平衡方程Kd＝F，求解自由曲面壳体结点位移d，并找出结点最大位移，记为d_max，将其作为约束条件；

其中，K表示自由曲面整体刚度矩阵，d表示自由曲面结点位移向量，F表示荷载向量；

(8-2)计算应变能

应变能等于曲面结点力乘以结点位移的1/2，计算公式为

(8-3)将应变能标准差σ{C}作为优化目标：

其中，var{C}为应变能的方差，NC为不确定性荷载个数，x_β表示第β个不确定性荷载，σ_β为第β个不确定荷载x_β的标准差。

作为本发明的进一步优化方案，步骤8中增广目标函数为：

其中，fun为增广Lagrange函数、是控制点高度z的函数，σ{C}表示应变能标准差，p表示惩罚因子，ω表示乘子，d_max表示结点最大位移，g(d_max)表示约束函数，即 g(d_max)＝d_max-d^*，d^*表示最大位移约束。

作为本发明的进一步优化方案，步骤9采用梯度法确定优化方向和步长来更新NURBS 曲面控制点坐标，获得新的自由曲面形状，具体步骤如下

记L为自由曲面结构的第L个优化步，且L属于正整数集，Δz^(L)表示第L个优化步和第 L+1个优化步的控制点高度差，其求解方法如下所示：

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

其中，m和n分别为u和v方向的控制点个数；z^(L)表示第L优化步控制点高度；表示第L优化步的增广目标函数fun关于第L优化步的控制点高度的梯度，求解方法如(10-1)所示；λ^(L)表示第L优化步的步长，其求解方法如(10-2)所示；

(10-1)增广函数fun差分算法求解梯度

增广目标函数fun是z^(L)的函数，假定增广目标函数fun是关于z^(L)的连续可微函数，将 fun(z^(L))行Taylor展开：

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

式中：Δz_i自由曲面控制点上纵坐标增量；o(Δz_i)是高阶无穷小量；且ξ∈(0,Δz_i)，fun⁽³⁾(ξ) 表示增广函数fun关于ξ的三阶导数；

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

其计算步骤如下，其中H表示黄金分割法的第H迭代步，且H属于正整数集，λ_H为第H次迭代求得的步长值，当程序循环至满足step2的条件时将λ_H的值赋给λ^(L):

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

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

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

step2:若b_k-a_k＜G，则停止计算；否则，当fun(λ_H)＞fun(μ_H)时，转step3；当 fun(λ_H)＜fun(μ_H)时，转step4；

step3：置a_H+1＝λ_H，b_H+1＝b_H，λ_H+1＝μ_H，μ_H+1＝a_H+1+0.618(b_H+1-a_H+1)；计算增广函数fun(μ_H+1)转step5；

step4：置a_H+1＝a_H，b_H+1＝μ_H，μ_H+1＝λ_H，λ_H+1＝a_H+1+0.382(b_H+1-a_H+1)置a_H+1＝ a_H，b_H+1＝μ_H，μ_H+1＝λ_H，λ_H+1＝a_H+1+0.382(b_H+1-a_H+1)；计算增广函数fun(λ_H+1)，转step5；

step5：置H＝H+1，返回step2；

(10-3)修正乘子法中的乘子ω

乘子ω的修正公式如下

ω^L+1＝max(0,ω^L-pg(d_max))

其中，L表示第L个优化步。

本发明提供一种考虑荷载不确定性的自由曲面结构形态创建方法，具有以下技术效果：

1.该方法算法简单有效，考虑自由曲面创建过程中的考虑荷载不确定性的影响；

2.不确定性优化后的自由曲面结构的刚度和鲁棒性都远优于初始结构，具有更好的受力性能；

3.当荷载的大小不变，方向不确定时，不确定性优化结构的鲁棒性和可靠性要远优于确定性优化结构的鲁棒性和可靠性。

附图说明

图1为考虑荷载不确定性的自由曲面结构形态创建方法的步骤流程图；

图2为算例计算模型侧视图；

图3为不确定性形状优化过程中自由曲面的应变能和鲁棒性指标的变化曲线；

图4为不确定性形状优化过程中自由曲面的形状变化情况，其中，(a)为第1(初始曲面) 优化步时的曲面形状，(b)为第10优化步时的曲面形状，(c)为第100优化步时的曲面形状， (d)为第200优化步时的曲面形状，(e)为第400优化步时的曲面形状，(f)为第527(最终曲面)优化步时的曲面形状；

图5为确定性形状优化过程中自由曲面的应变能和鲁棒性指标的变化曲线；

图6为不同方向的荷载作用下确定性优化结构和不确定性优化结构的应变能变化曲线对比。

具体实施方法

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

本发明提供一种考虑荷载不确定性的自由曲面结构形态创建方法，如图1所示，具体步骤如下：

步骤1，采用非均匀有理B样条技术NURBS根据给定的初始控制点坐标建立自由曲面初始几何模型。

步骤1所述的自由曲面是一张在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方向的结点矢量且u_i∈U,v_i为v方向的结点矢量且v_i∈V，

步骤2，采用三角形壳单元对步骤1得到的自由曲面初始几何模型进行网格划分，得到混凝土三角形壳单元。

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

所述步骤3具体步骤为

(3-1)混凝土单元采用三角形薄壳单元模拟，单元含3个结点，每个结点含6个自由度。

(3-2)混凝土单元采用三角形薄壳单元模拟，三角形薄壳单元的刚度矩阵由平面应力状态和弯曲应力状态的刚度矩阵得到，三角形薄壳单元的三个结点分别记为r、s、t。

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

其中，[k^p]表示局部坐标系下混凝土三角形平面应力单元的刚度矩阵；分别为局部坐标系下混凝土平面应力单元的r、s、t结点的结点力，且由u方向的结点力F_u和 v方向的结点力F_v组成，分别表示局部坐标系下混凝土三角形平面应力单元的r、s、t结点的位移，且由u方向的位移u_c和v方向的位移v_c组成；

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

其中，[k^b]表示局部坐标系下混凝土三角形弯曲应力单元的刚度矩阵；分别为局部坐标系下混凝土弯曲应力单元r、s、t结点的结点力，且由垂直于u、v平面w方向的结点力F_w、绕u轴旋转的结点力M_u及绕v轴旋转的结点力M_v组成；分别表示局部坐标系下混凝土三角形弯曲应力单元的r、s、t结点的位移，由垂直于u、v平面w 方向的线位移w_c，绕u轴旋转角θ_u及绕v轴旋转角θ_v组成。

混凝土三角形薄壳单元的局部坐标系下的刚度矩阵[k_c]为：

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

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

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

其中，[L]为坐标转换矩阵，[k]为局部坐标系下的单元刚度矩阵，L₁表示局部坐标系各坐标轴与整体坐标系坐标轴之间的方向余弦。

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

所述步骤5的具体步骤为，由整体坐标系下的单元刚度矩阵[k′]拼装得到整体坐标系下的整体刚度矩阵K。

步骤6，利用利用不确定性参数—不确定性荷载大小及方向的期望值和标准差，建立不确定性荷载的概率模型，对自由曲面结构施加不确定性荷载大小和/或方向的期望值。

所述步骤6具体步骤如下：

(6-1)采用不确定性荷载概率模型描述不确定性荷载UF：

假定不确定性荷载UF₁、UF₂……UF_NC均为正态分布随机变量，且相互独立，记为

其中，UF_β表示第β个不确定性荷载，N(·)表示正态分布符号，σ_β分别表示第β个不确定性荷载大小方向的期望值和标准差，NC表示不确定性荷载的个数。

(6-2)在结构上施加NC个荷载

步骤7，对有限元平衡方程进行求解，通过结点荷载、结点位移以及不确定性荷载概率模型中的不确定性参数计算出应变能标准差作为形态创建优化目标，将结点最大位移作为约束条件，形成有约束的优化问题。

所述步骤7具体步骤如下：

(7-1)求解有限元平衡方程

根据有限元平衡方程Kd＝F，求解自由曲面壳体结点位移d，并找出结点最大位移，记为d_max，将其作为约束条件。

其中，K表示自由曲面整体刚度矩阵，d表示自由曲面结点位移向量，F表示荷载向量。

(7-2)计算应变能

应变能等于曲面结点力乘以结点位移的1/2，计算公式为

(7-3)推导应变能标准差

将结构总应变能方程在荷载不确定参数期望值处进行一阶Taylor展开：

静荷载作用下的线性有限元平衡方程为Kd＝F，平衡方程两端同时关于各荷载不确定参数求导，如下式所示

将上式代入结构应变能的一阶Taylor展开式中得：

其中，x_i表示第i个不确定性荷载，NC为不确定性荷载个数。

根据应变能的一阶Taylor展开式，结构应变能的方差可表示为

最终，就可以得到结构的应变能标准差，并将其作为优化目标：

其中，σ_β为第β个不确定荷载x_β的标准差。

步骤8，采用广义拉格朗日乘子法将步骤7中的有约束优化问题转换为无约束问题，构成增广目标函数。

所述步骤8具体步骤为：

(8-1)建立增广目标函数

定义增广Lagrange函数

其中，fun为增广Lagrange函数，是控制点高度z的函数，σ{C}表示应变能标准差，p表示惩罚因子，ω表示乘子，d_max表示结点最大位移，g(d_max)表示约束函数，即 g(d_max)＝d_max-d^*，d^*表示最大位移约束。

(8-2)有约束优化问题转换为无约束问题

以增广Lagrange函数fun替换应变能标准差σ{C}作为目标优化函数。

步骤9，采用梯度法确定优化方向和步长来更新自由曲面控制点坐标，获得新的自由曲面形状。

所述步骤9具体步骤如下：

记L为自由曲面结构的第L个优化步，且L属于正整数集，Δz^(L)表示第L个优化步和第 L+1个优化步的控制点高度差，其求解方法如下所示：

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

其中，m和n分别为u和v方向的控制点个数；z^(L)表示第L优化步控制点高度；表示第L优化步的增广函数fun关于第L优化步的控制点高度的梯度，求解公式如(9-1)所示；λ^(L)表示第L优化步的步长，其求解公式如(9-2)所示。

(9-1)增广函数fun差分算法求解梯度

增广函数fun是z^(L)的函数，假定增广函数fun是关于z^(L)的连续可微函数，将fun(z^(L))行 Taylor展开：

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

式中：Δz_i自由曲面控制点上纵坐标增量；o(Δz_i)是高阶无穷小量；且ξ∈(0,Δz_i)，fun⁽³⁾(ξ) 表示增广函数fun关于ξ的三阶导数。

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

其计算步骤如下，其中H表示黄金分割法的第H迭代步，且H属于正整数集，λ_H为第H次迭代求得的步长值，当程序循环至满足step2的条件时将λ_H的值赋给λ^(L):

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

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

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

step2:若b_k-a_k＜G，则停止计算；否则，当fun(λ_H)＞fun(μ_H)时，转step3；当 fun(λ_H)＜fun(μ_H)时，转step4；

step3：置a_H+1＝λ_H，b_H+1＝b_H，λ_H+1＝μ_H，μ_H+1＝a_H+1+0.618(b_H+1-a_H+1)；计算增广函数fun(μ_H+1)转step5；

step4：置a_H+1＝a_H，b_H+1＝μ_H，μ_H+1＝λ_H，λ_H+1＝a_H+1+0.382(b_H+1-a_H+1)置a_H+1＝ a_H，b_H+1＝μ_H，μ_H+1＝λ_H，λ_H+1＝a_H+1+0.382(b_H+1-a_H+1)；计算增广函数fun(λ_H+1)，转step5；

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

(9-3)修正乘子法中的乘子ω

乘子ω的修正公式如下

ω^L+1＝max(0,ω^L-pg(d_max))

其中，L表示第L个优化步。

步骤10，判断本优化是否收敛的具体步骤为：设置优化精度ε^*，若当前优化步的增广目标函数对当前优化步的控制点高度的差分的模小于优化精度ε^*，则考虑确定性荷载的自由曲面优化结束，自由曲面创建完成；否则继续迭代至步骤2进行自由曲面结构形态优化。

具体实施例

下面结合附图和具体实施例，进一步阐明本发明。应理解下述具体实施方式仅用于说明本发明而不用于限制本发明的范围。

如图2所示一平面投影为正方形(12m×12m)的自由曲面结构，在四个角点简支。曲面厚度t＝0.1m，的弹性模量E＝3×1010MPa，泊松比ν＝0.2。控制点设置如图2所示。结构的控制点6、7、10、11处分别作用竖直向下的不确定集中荷载：荷载值f＝200kN为确定量；荷载方向(以荷载方向与Z轴正方向的夹角θ表示)为不确定量，且服从均值为π、标准差为0.5 的正态分布，记为

θ～N(π,0.5²)

在形态创建过程中，角点处控制点1、4、13、16位置不变，将其他控制点竖向坐标Z作为优化变量，以结点最大竖向位移≤0.015m作为约束条件，以鲁棒性指标(应变能标准差) 作为目标优化函数，对算例计算模型进行不确定性形状优化。有限元计算采用三角形薄壳单元，单元数为4608，结点数为2401。

图3给出了不确定性形状优化过程中结构的鲁棒性指标和应变能的变化曲线，其中横坐标为优化步，纵坐标为结构性能函数值；表1则给出了确定性形状优化前后结构的应变能和鲁棒性指标的变化情况；图4中的(a)至(f)则分别给出了自由曲面结构在第1(初始曲面)、 10、100、200、400、527(最终曲面)优化步时的曲面形状，显示了不确定性形状优化过程中的曲面形状变化。

表1不确定性形状优化前后结构性能函数变化

性能函数	初始结构(N·m)	最终结构(N·m)	下降率
				应变能	35106.43	3905.41	88.87％
鲁棒性指标	107.10	67.15	37.30％

通过图3和表1可以发现，不确定性优化后的自由曲面结构的刚度和鲁棒性都远优于初始结构，具有更好的受力性能。进一步地，为了说明本发明相比于确定性形状优化的优势，下面将以应变能作为目标优化函数，其他条件同上述不确定性形状优化，对图2所示算例模型进行确定性形状优化，并进行对比分析。

图5给出了确定性形状优化过程中结构的应变能和鲁棒性指标的变化曲线，其中横坐标为优化步，纵坐标为结构性能函数值；表2则给出了确定性形状优化前后结构的应变能和不确定性指标的对比。

表2确定性形状优化前后结构性能函数对比

性能函数	初始结构(N·m)	最终结构(N·m)	下降率
				应变能	35106.43	1531.48	95.64％
鲁棒性指标	107.10	185.35	-73.06％

图6则给出了在大小一定，方向不同的荷载作用下确定性优化结构和不确定性优化结构的应变能变化曲线对比。对比分析后可以得出结论：当荷载的大小不变，方向不确定时，不确定性优化结构的鲁棒性和可靠性要远优于确定性优化结构的鲁棒性和可靠性。

本发明方案所公开的技术手段不仅限于上述技术手段所公开的技术手段，还包括由以上技术特征任意组合所组成的技术方案。

以上述依据本发明的理想实施例为启示，通过上述的说明内容，相关工作人员完全可以在不偏离本项发明技术思想的范围内，进行多样的变更以及修改。本项发明的技术性范围并不局限于说明书上的内容，必须要根据权利要求范围来确定其技术性范围。

Claims (10)

1.一种考虑荷载不确定性的自由曲面结构形态创建方法，其特征在于，具体方法步骤如下：

步骤1，采用非均匀有理B样条技术NURBS根据给定的初始控制点坐标建立自由曲面初始几何模型；

步骤2，采用三角形壳单元对步骤1得到的自由曲面初始几何模型进行网格划分，得到混凝土三角形壳单元；

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

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

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

步骤6，利用不确定性参数、不确定性荷载大小及方向的期望值和标准差，建立不确定性荷载的概率模型，对自由曲面结构施加不确定性荷载大小和/或方向的期望值；

步骤7，对有限元平衡方程进行求解，通过结点荷载、结点位移以及不确定性荷载概率模型中的不确定性参数计算出应变能标准差作为形态创建优化目标，将结点最大位移作为约束条件，形成有约束的优化问题；

步骤8，采用广义拉格朗日乘子法将步骤7中的有约束优化问题转换为无约束问题，构成增广目标函数；

步骤9，采用梯度法确定优化方向和步长来更新自由曲面控制点坐标，获得新的自由曲面形状；

步骤10，设置优化精度ε^*，若当前优化步的增广目标函数对当前优化步的控制点高度的差分的模小于优化精度ε^*，则考虑确定性荷载的自由曲面优化结束，自由曲面创建完成；否则继续迭代至步骤2进行自由曲面结构形态优化。

2.根据权利要求1所述的一种考虑荷载不确定性的自由曲面结构形态创建方法，其特征在于，步骤1所述的自由曲面是NURBS曲面。

3.根据权利要求1所述的一种考虑荷载不确定性的自由曲面结构形态创建方法，其特征在于，自由曲面是一张在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方向的结点矢量且u_i∈U,v_i为v方向的结点矢量且v_i∈V，

4.根据权利要求1所述的一种考虑荷载不确定性的自由曲面结构形态创建方法，其特征在于，步骤3计算混凝土三角形壳单元局部坐标系下的单元刚度矩阵具体步骤为

(4-1)混凝土单元采用三角形薄壳单元模拟，单元含3个结点，每个结点含6个自由度；

(4-2)混凝土单元采用三角形薄壳单元模拟，三角形薄壳单元的刚度矩阵由平面应力状态和弯曲应力状态的刚度矩阵得到，三角形薄壳单元的三个结点分别记为r、s、t；

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

其中，[k^p]表示局部坐标系下混凝土三角形平面应力单元的刚度矩阵；分别为局部坐标系下混凝土平面应力单元的r、s、t结点的结点力，且由u方向的结点力F_u和v方向的结点力F_v组成，分别表示局部坐标系下混凝土三角形平面应力单元的r、s、t结点的位移，且由u方向的位移u_c和v方向的位移v_c组成；

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

其中，[k^b]表示局部坐标系下混凝土三角形弯曲应力单元的刚度矩阵；分别为局部坐标系下混凝土弯曲应力单元r、s、t结点的结点力，且由垂直于平面w方向的结点力F_w、绕u轴旋转的结点力M_u及绕轴旋转的结点力M_v组成；分别表示局部坐标系下混凝土三角形弯曲应力单元的r、s、t结点的位移，由垂直于平面w方向的线位移w_c，绕u轴旋转角θ_u及绕v轴旋转角θ_v组成；

混凝土三角形薄壳单元的局部坐标系下的刚度矩阵[k_c]为：

5.根据权利要求1所述的一种考虑荷载不确定性的自由曲面结构形态创建方法，其特征在于，步骤4中整体坐标下的单元刚度矩阵为[k′]：

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

其中，[L]为坐标转换矩阵，[k]表示局部坐标系下的单元刚度矩阵，L₁表示局部坐标系各坐标轴与整体坐标系坐标轴之间的方向余弦。

6.根据权利要求1所述的一种考虑荷载不确定性的自由曲面结构形态创建方法，其特征在于，步骤5的具体步骤为，由整体坐标系下的单元刚度矩阵[k′]拼装得到整体坐标系下的整体刚度矩阵K。

7.根据权利要求1所述的一种考虑荷载不确定性的自由曲面结构形态创建方法，其特征在于，步骤6中利用不确定性参数、不确定性荷载大小及方向的期望值和标准差，建立不确定性荷载的概率模型，对结构施加不确定性荷载大小和/或方向的期望值，具体步骤如下

(7-1)采用概率模型描述不确定性荷载UF：

假定不确定性荷载UF₁、UF₂……UF_NC均为正态分布随机变量，且相互独立，记为

其中，UF_β表示第β个不确定性荷载，N(·)表示正态分布符号，分别表示第β个不确定性荷载大小方向的期望值和标准差，NC表示不确定性荷载的个数；

(7-2)在结构上施加NC个荷载

8.根据权利要求1所述的一种考虑荷载不确定性的自由曲面结构形态创建方法，其特征在于，步骤7对有限元平衡方程进行求解，通过结点荷载、结点位移以及概率模型中的不确定性参数计算出应变能标准差作为形态创建优化目标，将结点最大位移作为约束条件，具体步骤如下

(8-1)求解有限元平衡方程

根据有限元平衡方程Kd＝F，求解自由曲面壳体结点位移d，并找出结点最大位移，记为d_max，将其作为约束条件；

其中，K表示自由曲面整体刚度矩阵，d表示自由曲面结点位移向量，F表示荷载向量；

(8-2)计算应变能

应变能等于曲面结点力乘以结点位移的1/2，计算公式为

(8-3)将应变能标准差σ{C}作为优化目标：

其中，var{C}为应变能的方差，NC为不确定性荷载个数，x_β表示第β个不确定性荷载，σ_β为第β个不确定荷载x_β的标准差。

9.根据权利要求1所述的一种考虑荷载不确定性的自由曲面结构形态创建方法，其特征在于，步骤8中增广目标函数为：

其中，fun为增广Lagrange函数、是控制点高度z的函数，σ{C}表示应变能标准差，p表示惩罚因子，ω表示乘子，d_max表示结点最大位移，g(d_max)表示约束函数，即g(d_max)＝d_max-d^*，d^*表示最大位移约束。

10.根据权利要求1所述的一种考虑荷载不确定性的自由曲面结构形态创建方法，其特征在于，步骤9采用梯度法确定优化方向和步长来更新NURBS曲面控制点坐标，获得新的自由曲面形状，具体步骤如下

记L为自由曲面结构的第L个优化步，且L属于正整数集，Δz^(L)表示第L个优化步和第L+1个优化步的控制点高度差，其求解方法如下所示：

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

其中，m和n分别为u和v方向的控制点个数；z^(L)表示第L优化步控制点高度；表示第L优化步的增广目标函数fun关于第L优化步的控制点高度的梯度，求解方法如(10-1)所示；λ^(L)表示第L优化步的步长，其求解方法如(10-2)所示；

(10-1)增广函数fun差分算法求解梯度

增广目标函数fun是z^(L)的函数，假定增广目标函数fun是关于z^(L)的连续可微函数，将fun(z^(L))行Taylor展开：

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

式中：Δz_i自由曲面控制点上纵坐标增量；o(Δz_i)是高阶无穷小量；且ξ∈(0,Δz_i)，fun⁽³⁾(ξ)表示增广函数fun关于ξ的三阶导数；

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

其计算步骤如下，其中H表示黄金分割法的第H迭代步，且H属于正整数集，λ_H为第H次迭代求得的步长值，当程序循环至满足step2的条件时将λ_H的值赋给λ^(L):

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

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

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

step2:若b_k-a_k＜G，则停止计算；否则，当fun(λ_H)＞fun(μ_H)时，转step3；当fun(λ_H)＜fun(μ_H)时，转step4；

step3：置a_H+1＝λ_H，b_H+1＝b_H，λ_H+1＝μ_H，μ_H+1＝a_H+1+0.618(b_H+1-a_H+1)；计算增广函数fun(μ_H+1)转step5；

step4：置a_H+1＝a_H，b_H+1＝μ_H，μ_H+1＝λ_H，λ_H+1＝a_H+1+0.382(b_H+1-a_H+1)置a_H+1＝a_H，b_H+1＝μ_H，μ_H+1＝λ_H，λ_H+1＝a_H+1+0.382(b_H+1-a_H+1)；计算增广函数fun(λ_H+1)，转step5；

step5：置H＝H+1，返回step2；

(10-3)修正乘子法中的乘子ω

乘子ω的修正公式如下

ω^L+1＝max(0,ω^L-pg(d_max))

其中，L表示第L个优化步。