CN107274490A - 一种考虑材料损伤和钢筋作用的自由曲面结构形态创建方法 - Google Patents

Abstract

本发明提供一种考虑材料损伤和钢筋作用的自由曲面形态创建方法，考虑到随着荷载的施加和变形的积累，材料在自由曲面创建过程中出现的微观损伤以及钢筋作用对自由曲面形态的影响，具体步骤如下：1，得到初始化自由曲面；2，对初始化自由曲面进行网格划分，得到钢筋混凝土三角形壳单元；3，计算钢筋混凝土三角形壳单元的局部坐标系下的刚度矩阵；4，计算整体坐标下的单元刚度矩阵；5，得到自由曲面结构的整体刚度矩阵；6，将Mazars损伤模型引入混凝土单元、将双线性模型引入钢筋单元，对结构施加荷载，步骤7，计算应变能；8，对考虑Mazars损伤模型和钢筋作用的自由曲面的控制点进行高度调整；步骤9，设置优化精度ε^*。本发明简单有效。

Description

一种考虑材料损伤和钢筋作用的自由曲面结构形态创建方法

技术领域

本发明涉及一种考虑材料损伤和钢筋作用的自由曲面结构形态创建方法。

背景技术

国内外目前自由曲面结构的形态创建方法可以概括为以下两种。1.基于试验的物理方法，即以工程经验或物理实验为主要手段，有意识地确定合理结构形状。2.基于数值的优化方法，即基于计算机技术，以计算机图形学和结构优化等为手段创建合理结构形态。

在现有的自由曲面形态结构创建方法中，自由曲面的模型都采用纯混凝土线弹性的方法进行分析和优化。此类方法没有考虑到随着荷载的施加和变形的积累，材料在自由曲面创建过程中出现的微观损伤，也没有考虑钢筋的作用对自由曲面形态创建结果的影响。基于此，提出考虑材料损伤和钢筋作用的自由曲面形态创建办法。

发明内容

本发明的目是解决上述问题，提供一种考虑材料损伤和钢筋作用的自由曲面形态创建方法。该方法算法简单有效，考虑自由曲面创建过程中的材料损伤和钢筋的作用。为实现上述目的，本发明采用如下技术方案。

一种考虑材料损伤和钢筋作用的自由曲面创建方法，具体步骤为：

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

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

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

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

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

步骤6，将Mazars损伤模型引入混凝土单元、将双线性模型引入钢筋单元，对结构施加荷载，分N个荷载步进行加载，每个荷载步内存在M个迭代步，计算自由曲面壳体位移；

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

步骤8，对考虑Mazars损伤模型和钢筋作用的自由曲面的控制点进行高度调整；

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

步骤1中的自由曲面是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，

所述步骤3计算步骤2得到的钢筋混凝土三角形壳单元的局部坐标系下的刚度矩阵具体步骤为：

钢筋混凝土三角形单元分为混凝土单元和钢筋单元，混凝土单元和钢筋单元共结点，两者刚度矩阵叠加得到钢筋混凝土三角形单元的刚度矩阵。混凝土单元采用三角形薄壳单元模拟，单元含3个结点，每个结点含6个自由度。钢筋单元采用三角形平面应力单元模拟，单元含3个结点，每个结点含2个自由度。

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

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

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

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

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

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

(3-2)钢筋单元采用三角形平面应力单元模拟，在钢筋三角形平面应力单元中，结点力和结点位移的关系如下：

其中，表示局部坐标系下钢筋三角形平面应力单元的刚度矩阵； 分别为局部坐标系下钢筋平面应力单元的r、s、t结点的结点力，且由u方向的结点力和v方向的结点力组成；分别表示局部坐标系下钢筋平面应力单元的r、s、t结点的位移，且由u方向的位移和v方向的位移组成。

(3-3)得到钢筋混凝土三角形单元刚度矩阵的局部坐标系下的刚度矩阵[k]为：

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

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

其中，[L]为坐标转换矩阵；

由整体坐标系下的单元刚度矩阵[k′]拼装得到整体坐标系下的整体刚度矩阵[K]。

步骤6中计算自由曲面壳体位移的具体步骤如下：

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

对于第I荷载步，施加荷载{P}/N，{P}表示荷载，N表示荷载步数，I∈[1，N]且I属于正整数集；根据[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}，计算前I荷载步混凝土单元的累计应力向量并由求解前I荷载步混凝土单元累计主应力向量和累计主应变向量求解第I+1步的损伤因子D_I+1，最后求解第I+1荷载步混凝土单元刚度矩阵，具体如下：

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

混凝土平面应力单元的应力向量为：

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

混凝土弯曲应力单元的应力向量为：

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

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

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

其中，{σ_Ir}、{σ_Is}、{σ_It}分别表示第I荷载步混凝土单元r、s、t结点的总应力向量，

(6-3-2)根据(6-3-1)得到的第I荷载步的混凝土单元总应力向量{σ_I}，计算前I荷载步混凝土单元的累计应力向量具体计算方法如下：

计算前I荷载步混凝土单元的累计应力向量

其中，混凝土单元r结点的累计应力向量s结点的累计应力向量t结点的累计应力向量

(6-3-3)根据(6-3-2)得到的前I荷载步混凝土单元的累计应力向量求前I荷载步混凝土单元累计主应力向量

前I荷载步混凝土单元r、s、t结点累计应力向量和r、s、t结点主应力向量如下：

r结点：

第一主应力：

第二主应力：

第三主应力：

其中，r结点的累计应力向量用u方向累计应力v方向累计应力和垂直u、v平面的w方向累计应力表达为

s结点：

第一主应力：

第二主应力：

第三主应力：

其中，s结点的累计应力向量用u方向累计应力v方向累计应力和垂直u、v平面的w方向累计应力表达为

t结点：

第一主应力：

第二主应力：

第三主应力：

其中，t结点的累计应力向量用u方向累计应力v方向累计应力和垂直u、v平面的w方向累计应力表达为

(6-3-4)根据(6-3-3)得到的混凝土单元累计单元主应力求混凝土单元累计单元主应变

其中，混凝土单元r结点的结点主应变记为 分别表示r结点第一、第二、第三主应变， s结点的结点主应变记为 分别表示s结点第一、第二、第三主应变， t结点的结点主应变记为分别表示t结点第一、第二、第三主应变，E_c为混凝土的弹性模量；

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

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

结点损伤因子是结点最大主应变ε的函数，其中r结点的损伤因子计算方法如下，由于只讨论

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

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

当大于0，小于0时，则单元处于复杂应力状态，分别取和代入以下情况3的ε值；

情况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.⁵，1000＜B_c＜2000。

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

情况3：复杂应力状态情况下，拉、压损伤之间存在耦合效应，Mazars建议用下时描述单元总损伤

其中，α_T和α_C为耦合系数，纯压缩时取α_T＝0，纯拉伸时取α_C＝0，组合情况α_T+α_C=1，D_T和D_C为耦合损伤因子；

采用相同方法计算s、t结点得损伤因子

(6-3-6)求解第I+1步的局部坐标和整体坐标系下混凝土单元刚度矩阵：

第I+1步的局部坐系下标混凝土单元刚度矩阵计算方法如下：

其中，

第I+1步的整体坐标系的混凝土单元刚度矩阵形成方法同上述步骤4所述。

(6-4)由局部坐标系下的计算钢筋平面应力单元的应力向量根据计算前I荷载步的钢筋平面应力单元累计应力向量并由求解前I荷载步累计主应力向量和累计主应变向量由累计主应变向量判断钢筋是否屈服并求解第I+1荷载步钢筋单元刚度矩阵，具体如下：

(6-4-1)由局部坐标系下的计算钢筋单元平面应力单元的应力向量钢筋平面应力单元的应力向量为：

其中，表示钢筋平面应力单元的弹性矩阵，表示钢筋平面应力单元的应变矩阵；

(6-4-2)根据(6-4-1)得到的第I荷载步钢筋单元平面应力单元的应力向量计算前I荷载步钢筋单元的累计应力向量具体计算方法如下：

计算前I荷载步钢筋单元的累计应力向量

其中，钢筋单元r结点的累计应力向量s结点的累计应力向量t结点的累计应力向量

(6-4-3)根据(6-4-2)得到的前I荷载步钢筋单元的累计应力向量求前I荷载步钢筋单元累计主应力向量

前I荷载步钢筋单元r、s、t结点累计应力向量和r、s、t结点主应力向量如下：

r结点：

第一主应力：

第二主应力：

第三主应力：

其中，r结点的累计应力向量用u方向累计应力v方向累计应力和垂直u、v平面的w方向累计应力表达为

s结点：

第一主应力：

第二主应力：

第三主应力：

其中，s结点的累计应力向量用u方向累计应力v方向累计应力和垂直u、v平面的w方向累计应力表达为

t结点：

第一主应力：

第二主应力：

第三主应力：

其中，t结点的累计应力向量用u方向累计应力v方向累计应力和垂直u、v平面的w方向累计应力表达为

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

其中，钢筋单元r结点的结点主应变记为向量 分别表示r结点第一、第二、第三主应变， s结点的结点主应变记为向量 分别表示s结点第一、第二、第三主应变， t结点的结点主应变记为向量分别表示t结点第一、第二、第三主应变，E_s为钢筋的弹性模量；

(6-4-5)求解第I+1步钢筋单元的弹性模量和刚度矩阵：

取钢筋单元主累计应变中绝对值最大值为钢筋单元最大主应变为ε_s；

若ε_s≤f_y/E_s，则钢筋单元未屈服，第I+1步钢筋单元的弹性模量为E_s，钢筋单元的刚度矩阵为：

其中，表示钢筋平面应力单元的弹性矩阵；

若ε_s＞f_y/E_s，则钢筋单元屈服，第I+1步钢筋单元的弹性模量为E_s2，钢筋单元的刚度矩阵为：

其中，表示钢筋平面应力单元的弹性矩阵。

(6-5)求解第I+1步钢筋混凝土单元刚度矩阵[k_I+1]：

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

在整体坐标下，记第J迭代步的位移误差向量为且J∈[1，M]；M为迭代次数，在程序运行时设定；

对于第I荷载步第1迭代步的公式如下：

可以求解{δ′_I}，对于第J迭代步公式如下：

对于第M迭代步的迭代公式如下：

(6-7)由上述步骤可以得到同时考虑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_N}。

步骤8所述的对考虑Mazars损伤模型和钢筋作用的自由曲面的控制点进行高度调整，具体步骤为：

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

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

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

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

结构应变能C是z^(L)的函数，假定应变能C是关于z^(L)的连续可微函数，将C(z^(L))进行Taylor展开：

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

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

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

其计算步骤如下，其中F表示黄金分割法的第F迭代步，且F属于正整数集，λ_F为第F次迭代求得的步长值，当程序循环至满足step2的条件时将λ_F的值赋给λ^(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为根据给定的初始自由曲面的一组控制点坐标得到初始化自由曲面，

其中：图2(a)为该自由结构的正视图，

图2(b)为该自由结构的俯视图；

图3(a)为考虑和不考虑材料损伤和钢筋作用的曲面创建结果的曲面形态，

图3(b)为不考虑材料损伤和钢筋作用的曲面创建结果的曲面形态，

图4(a)考虑材料损伤和钢筋作用的创建结果的曲面损伤因子，

图4(b)不考虑材料损伤和钢筋作用的创建结果的曲面损伤因子，

图5(a)考虑材料损伤和钢筋作用的创建结果的曲面位移，

图5(b)不考虑材料损伤和钢筋作用的创建结果的曲面位移，

图6(a)考虑材料损伤和钢筋作用的创建结果的曲面Mises应力，

图6(b)不考虑材料损伤和钢筋作用的创建结果的曲面Mises应力。

具体实施方式

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

附图1为考虑材料损伤和钢筋作用的自由曲面结构形态创建方法的步骤流程图；结合附图可见，本方法具体步骤为：

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

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

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

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

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

步骤6，将Mazars损伤模型引入混凝土单元、将双线性模型引入钢筋单元，对结构施加荷载，分N个荷载步进行加载，每个荷载步内存在M个迭代步，计算自由曲面壳体位移；

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

步骤8，对考虑Mazars损伤模型和钢筋作用的自由曲面的控制点进行高度调整；

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

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

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

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

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

附图4-图6分别为考虑不考虑材料损伤和钢筋作用的自由曲面创建结果的损伤因子、结点位移和Mises应力对比，可以发现采用本方法创建的自由曲面在结构损伤因子、整体刚度和力学响应方面表现更优。

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

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

Claims (7)

1.一种考虑材料损伤和钢筋作用的自由曲面形态创建方法，其特征在于，具体方法步骤如下：

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

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

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

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

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

步骤6，将Mazars损伤模型引入混凝土单元、将双线性模型引入钢筋单元，对结构施加荷载，分N个荷载步进行加载，每个荷载步内存在M个迭代步，计算自由曲面壳体位移；

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

步骤8，对考虑Mazars损伤模型和钢筋作用的自由曲面的控制点进行高度调整；

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

2.根据权利要求1所述的一种考虑材料损伤和钢筋作用的自由曲面形态创建方法，其特征在于，步骤1中的自由曲面是NURBS曲面。

3.根据权利要求1所述的一种考虑材料损伤和钢筋作用的自由曲面形态创建方法，其特征在于，自由曲面是一张在u方向p次、v方向q次的NURBS曲面，具有如下形式的双变量分段有理矢值函数：

<mrow> <mi>S</mi> <mrow> <mo>(</mo> <mi>u</mi> <mo>,</mo> <mi>v</mi> <mo>)</mo> </mrow> <mo>=</mo> <munderover> <mi>&amp;Sigma;</mi> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>n</mi> </munderover> <munderover> <mi>&amp;Sigma;</mi> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>m</mi> </munderover> <msub> <mi>N</mi> <mrow> <mi>i</mi> <mo>,</mo> <mi>p</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>u</mi> <mo>)</mo> </mrow> <msub> <mi>N</mi> <mrow> <mi>j</mi> <mo>,</mo> <mi>q</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>v</mi> <mo>)</mo> </mrow> <msub> <mi>w</mi> <mrow> <mi>i</mi> <mo>,</mo> <mi>j</mi> </mrow> </msub> <msub> <mi>P</mi> <mrow> <mi>i</mi> <mo>,</mo> <mi>j</mi> </mrow> </msub> <mo>/</mo> <munderover> <mi>&amp;Sigma;</mi> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>n</mi> </munderover> <munderover> <mi>&amp;Sigma;</mi> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>m</mi> </munderover> <msub> <mi>N</mi> <mrow> <mi>i</mi> <mo>,</mo> <mi>p</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>u</mi> <mo>)</mo> </mrow> <msub> <mi>N</mi> <mrow> <mi>j</mi> <mo>,</mo> <mi>q</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>v</mi> <mo>)</mo> </mrow> <msub> <mi>P</mi> <mrow> <mi>i</mi> <mo>,</mo> <mi>j</mi> </mrow> </msub> </mrow>

其中，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，

<mrow> <msub> <mi>N</mi> <mrow> <mi>i</mi> <mo>,</mo> <mi>p</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>u</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mrow> <mi>u</mi> <mo>-</mo> <msub> <mi>u</mi> <mi>i</mi> </msub> </mrow> <mrow> <msub> <mi>u</mi> <mrow> <mi>i</mi> <mo>+</mo> <mi>p</mi> </mrow> </msub> <mo>-</mo> <msub> <mi>u</mi> <mi>i</mi> </msub> </mrow> </mfrac> <msub> <mi>N</mi> <mrow> <mi>i</mi> <mo>,</mo> <mi>p</mi> <mo>-</mo> <mn>1</mn> </mrow> </msub> <mrow> <mo>(</mo> <mi>u</mi> <mo>)</mo> </mrow> <mo>+</mo> <mfrac> <mrow> <msub> <mi>u</mi> <mrow> <mi>i</mi> <mo>+</mo> <mi>p</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>-</mo> <mi>u</mi> </mrow> <mrow> <msub> <mi>u</mi> <mrow> <mi>i</mi> <mo>+</mo> <mi>p</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>-</mo> <msub> <mi>u</mi> <mrow> <mi>i</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> </mrow> </mfrac> <msub> <mi>N</mi> <mrow> <mi>i</mi> <mo>+</mo> <mn>1</mn> <mo>,</mo> <mi>p</mi> <mo>-</mo> <mn>1</mn> </mrow> </msub> <mrow> <mo>(</mo> <mi>u</mi> <mo>)</mo> </mrow> </mrow>

<mrow> <msub> <mi>N</mi> <mrow> <mi>j</mi> <mo>,</mo> <mi>q</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>v</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mrow> <mi>v</mi> <mo>-</mo> <msub> <mi>v</mi> <mi>j</mi> </msub> </mrow> <mrow> <msub> <mi>v</mi> <mrow> <mi>j</mi> <mo>+</mo> <mi>q</mi> </mrow> </msub> <mo>-</mo> <msub> <mi>v</mi> <mi>j</mi> </msub> </mrow> </mfrac> <msub> <mi>N</mi> <mrow> <mi>j</mi> <mo>,</mo> <mi>q</mi> <mo>-</mo> <mn>1</mn> </mrow> </msub> <mrow> <mo>(</mo> <mi>v</mi> <mo>)</mo> </mrow> <mo>+</mo> <mfrac> <mrow> <msub> <mi>v</mi> <mrow> <mi>j</mi> <mo>+</mo> <mi>q</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>-</mo> <mi>v</mi> </mrow> <mrow> <msub> <mi>v</mi> <mrow> <mi>j</mi> <mo>+</mo> <mi>q</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>-</mo> <msub> <mi>v</mi> <mrow> <mi>j</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> </mrow> </mfrac> <msub> <mi>N</mi> <mrow> <mi>j</mi> <mo>+</mo> <mn>1</mn> <mo>,</mo> <mi>q</mi> <mo>-</mo> <mn>1</mn> </mrow> </msub> <mrow> <mo>(</mo> <mi>v</mi> <mo>)</mo> </mrow> </mrow>

4.根据权利要求1所述的一种考虑材料损伤和钢筋作用的自由曲面形态创建方法，其特征在于，步骤3中计算钢筋混凝土三角形单元的局部坐标系下的刚度矩阵，具体为：

钢筋混凝土三角形单元分为混凝土单元和钢筋单元；混凝土单元和钢筋单元共结点，两者刚度矩阵叠加得到钢筋混凝土三角形单元的刚度矩阵；混凝土单元采用三角形薄壳单元模拟，单元含3个结点，每个结点含6个自由度；钢筋单元采用三角形平面应力单元模拟，单元含3个结点，每个结点含2个自由度；

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

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

<mrow> <mfenced open = "{" close = "}"> <mtable> <mtr> <mtd> <msubsup> <mi>F</mi> <mi>r</mi> <mi>p</mi> </msubsup> </mtd> </mtr> <mtr> <mtd> <msubsup> <mi>F</mi> <mi>s</mi> <mi>p</mi> </msubsup> </mtd> </mtr> <mtr> <mtd> <msubsup> <mi>F</mi> <mi>t</mi> <mi>p</mi> </msubsup> </mtd> </mtr> </mtable> </mfenced> <mo>=</mo> <mo>&amp;lsqb;</mo> <msubsup> <mi>k</mi> <mi>c</mi> <mi>p</mi> </msubsup> <mo>&amp;rsqb;</mo> <mfenced open = "{" close = "}"> <mtable> <mtr> <mtd> <msubsup> <mi>&amp;delta;</mi> <mi>r</mi> <mi>p</mi> </msubsup> </mtd> </mtr> <mtr> <mtd> <msubsup> <mi>&amp;delta;</mi> <mi>s</mi> <mi>p</mi> </msubsup> </mtd> </mtr> <mtr> <mtd> <msubsup> <mi>&amp;delta;</mi> <mi>t</mi> <mi>p</mi> </msubsup> </mtd> </mtr> </mtable> </mfenced> </mrow>

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

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

<mrow> <mfenced open = "{" close = "}"> <mtable> <mtr> <mtd> <msubsup> <mi>F</mi> <mi>r</mi> <mi>b</mi> </msubsup> </mtd> </mtr> <mtr> <mtd> <msubsup> <mi>F</mi> <mi>s</mi> <mi>b</mi> </msubsup> </mtd> </mtr> <mtr> <mtd> <msubsup> <mi>F</mi> <mi>t</mi> <mi>b</mi> </msubsup> </mtd> </mtr> </mtable> </mfenced> <mo>=</mo> <mo>&amp;lsqb;</mo> <msubsup> <mi>k</mi> <mi>c</mi> <mi>b</mi> </msubsup> <mo>&amp;rsqb;</mo> <mfenced open = "{" close = "}"> <mtable> <mtr> <mtd> <msubsup> <mi>&amp;delta;</mi> <mi>r</mi> <mi>b</mi> </msubsup> </mtd> </mtr> <mtr> <mtd> <msubsup> <mi>&amp;delta;</mi> <mi>s</mi> <mi>b</mi> </msubsup> </mtd> </mtr> <mtr> <mtd> <msubsup> <mi>&amp;delta;</mi> <mi>t</mi> <mi>b</mi> </msubsup> </mtd> </mtr> </mtable> </mfenced> </mrow>

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

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

<mrow> <mo>&amp;lsqb;</mo> <msub> <mi>k</mi> <mi>c</mi> </msub> <mo>&amp;rsqb;</mo> <mo>=</mo> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <mrow> <mo>&amp;lsqb;</mo> <msubsup> <mi>k</mi> <mi>c</mi> <mi>p</mi> </msubsup> <mo>&amp;rsqb;</mo> </mrow> </mtd> <mtd> <mrow></mrow> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> </mtr> <mtr> <mtd> <mrow></mrow> </mtd> <mtd> <mrow></mrow> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mrow></mrow> </mtd> <mtd> <mrow></mrow> </mtd> <mtd> <mrow></mrow> </mtd> <mtd> <mn>0</mn> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mrow></mrow> </mtd> <mtd> <mrow> <mo>&amp;lsqb;</mo> <msubsup> <mi>k</mi> <mi>c</mi> <mi>b</mi> </msubsup> <mo>&amp;rsqb;</mo> </mrow> </mtd> <mtd> <mrow></mrow> </mtd> <mtd> <mn>0</mn> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mrow></mrow> </mtd> <mtd> <mrow></mrow> </mtd> <mtd> <mrow></mrow> </mtd> <mtd> <mn>0</mn> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> </mtr> </mtable> </mfenced> </mrow>

(3-2)钢筋单元采用三角形平面应力单元模拟，在钢筋三角形平面应力单元中，结点力和结点位移的关系如下：

<mrow> <mfenced open = "{" close = "}"> <mtable> <mtr> <mtd> <msubsup> <mover> <mi>F</mi> <mo>&amp;OverBar;</mo> </mover> <mi>r</mi> <mi>p</mi> </msubsup> </mtd> </mtr> <mtr> <mtd> <msubsup> <mover> <mi>F</mi> <mo>&amp;OverBar;</mo> </mover> <mi>s</mi> <mi>p</mi> </msubsup> </mtd> </mtr> <mtr> <mtd> <msubsup> <mover> <mi>F</mi> <mo>&amp;OverBar;</mo> </mover> <mi>t</mi> <mi>p</mi> </msubsup> </mtd> </mtr> </mtable> </mfenced> <mo>=</mo> <mo>&amp;lsqb;</mo> <msubsup> <mi>k</mi> <mi>s</mi> <mi>p</mi> </msubsup> <mo>&amp;rsqb;</mo> <mfenced open = "{" close = "}"> <mtable> <mtr> <mtd> <msubsup> <mover> <mi>&amp;delta;</mi> <mo>&amp;OverBar;</mo> </mover> <mi>r</mi> <mi>p</mi> </msubsup> </mtd> </mtr> <mtr> <mtd> <msubsup> <mover> <mi>&amp;delta;</mi> <mo>&amp;OverBar;</mo> </mover> <mi>s</mi> <mi>p</mi> </msubsup> </mtd> </mtr> <mtr> <mtd> <msubsup> <mover> <mi>&amp;delta;</mi> <mo>&amp;OverBar;</mo> </mover> <mi>t</mi> <mi>p</mi> </msubsup> </mtd> </mtr> </mtable> </mfenced> </mrow>

其中，表示局部坐标系下钢筋三角形平面应力单元的刚度矩阵； 分别为局部坐标系下钢筋平面应力单元的r、s、t结点的结点力，且由u方向的结点力和v方向的结点力组成；分别表示局部坐标系下钢筋平面应力单元的r、s、t结点的位移，且由u方向的位移和v方向的位移组成；

(3-3)得到钢筋混凝土三角形单元刚度矩阵的局部坐标系下的刚度矩阵[k]为：

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

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

其中，[L]为坐标转换矩阵；

由整体坐标系下的单元刚度矩阵[k’]拼装得到整体坐标系下的整体刚度矩阵[K]。

6.根据权利要求1所述的一种考虑材料损伤的自由曲面形态创建方法，其特征在于，步骤6中计算自由曲面壳体位移的具体步骤如下：

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

对于第I荷载步，施加荷载{P}/N，{P}表示荷载，N表示荷载步数，I∈[1，N]且I属于正整数集；根据[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}，计算前I荷载步混凝土单元的累计应力向量并由求解前I荷载步混凝土单元累计主应力向量和累计主应变向量求解第I+1步的损伤因子D_I+1，最后求解第I+1荷载步混凝土单元刚度矩阵，具体如下：

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

混凝土平面应力单元的应力向量为：

<mrow> <mo>{</mo> <msubsup> <mi>&amp;sigma;</mi> <mi>I</mi> <mi>p</mi> </msubsup> <mo>}</mo> <mo>=</mo> <msup> <mrow> <mo>{</mo> <msup> <mrow> <mo>{</mo> <msubsup> <mi>&amp;sigma;</mi> <mi>r</mi> <mi>p</mi> </msubsup> <mo>}</mo> </mrow> <mi>T</mi> </msup> <mo>,</mo> <msup> <mrow> <mo>{</mo> <msubsup> <mi>&amp;sigma;</mi> <mi>s</mi> <mi>p</mi> </msubsup> <mo>}</mo> </mrow> <mi>T</mi> </msup> <mo>,</mo> <msup> <mrow> <mo>{</mo> <msubsup> <mi>&amp;sigma;</mi> <mi>t</mi> <mi>p</mi> </msubsup> <mo>}</mo> </mrow> <mi>T</mi> </msup> <mo>}</mo> </mrow> <mi>T</mi> </msup> </mrow>

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

混凝土弯曲应力单元的应力向量为：

<mrow> <mo>{</mo> <msubsup> <mi>&amp;sigma;</mi> <mi>I</mi> <mi>b</mi> </msubsup> <mo>}</mo> <mo>=</mo> <msup> <mrow> <mo>{</mo> <msup> <mrow> <mo>{</mo> <msubsup> <mi>&amp;sigma;</mi> <mi>r</mi> <mi>b</mi> </msubsup> <mo>}</mo> </mrow> <mi>T</mi> </msup> <mo>,</mo> <msup> <mrow> <mo>{</mo> <msubsup> <mi>&amp;sigma;</mi> <mi>s</mi> <mi>b</mi> </msubsup> <mo>}</mo> </mrow> <mi>T</mi> </msup> <mo>,</mo> <msup> <mrow> <mo>{</mo> <msubsup> <mi>&amp;epsiv;</mi> <mi>t</mi> <mi>b</mi> </msubsup> <mo>}</mo> </mrow> <mi>T</mi> </msup> <mo>}</mo> </mrow> <mi>T</mi> </msup> </mrow>

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

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

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

其中，{σ_Is}、{σ_Is}、{σ_k}分别表示第I荷载步混凝土单元r、s、t结点的总应力向量，

(6-3-2)根据(6-3-1)得到的第I荷载步的混凝土单元总应力向量{σ_I}，计算前I荷载步混凝土单元的累计应力向量具体计算方法如下：

计算前I荷载步混凝土单元的累计应力向量

<mrow> <mo>{</mo> <msubsup> <mi>&amp;sigma;</mi> <mi>I</mi> <mo>*</mo> </msubsup> <mo>}</mo> <mo>=</mo> <msup> <mrow> <mo>{</mo> <msup> <mrow> <mo>{</mo> <msubsup> <mi>&amp;sigma;</mi> <mrow> <mi>I</mi> <mi>r</mi> </mrow> <mo>*</mo> </msubsup> <mo>}</mo> </mrow> <mi>T</mi> </msup> <mo>,</mo> <msup> <mrow> <mo>{</mo> <msubsup> <mi>&amp;sigma;</mi> <mrow> <mi>I</mi> <mi>s</mi> </mrow> <mo>*</mo> </msubsup> <mo>}</mo> </mrow> <mi>T</mi> </msup> <mo>,</mo> <msup> <mrow> <mo>{</mo> <msubsup> <mi>&amp;sigma;</mi> <mrow> <mi>I</mi> <mi>t</mi> </mrow> <mo>*</mo> </msubsup> <mo>}</mo> </mrow> <mi>T</mi> </msup> <mo>}</mo> </mrow> <mi>T</mi> </msup> </mrow>

其中，混凝土单元r结点的累计应力向量s结点的累计应力向量t结点的累计应力向量

(6-3-3)根据(6-3-2)得到的前I荷载步混凝土单元的累计应力向量求前I荷载步混凝土单元累计主应力向量

<mrow> <mo>{</mo> <msubsup> <mi>&amp;sigma;</mi> <mi>M</mi> <mo>*</mo> </msubsup> <mo>}</mo> <mo>=</mo> <msup> <mrow> <mo>{</mo> <msup> <mrow> <mo>{</mo> <msubsup> <mi>&amp;sigma;</mi> <mrow> <mi>M</mi> <mi>r</mi> </mrow> <mo>*</mo> </msubsup> <mo>}</mo> </mrow> <mi>T</mi> </msup> <mo>,</mo> <msup> <mrow> <mo>{</mo> <msubsup> <mi>&amp;sigma;</mi> <mrow> <mi>M</mi> <mi>s</mi> </mrow> <mo>*</mo> </msubsup> <mo>}</mo> </mrow> <mi>T</mi> </msup> <mo>,</mo> <msup> <mrow> <mo>{</mo> <msubsup> <mi>&amp;sigma;</mi> <mrow> <mi>M</mi> <mi>t</mi> </mrow> <mo>*</mo> </msubsup> <mo>}</mo> </mrow> <mi>T</mi> </msup> <mo>}</mo> </mrow> <mi>T</mi> </msup> </mrow>

前I荷载步混凝土单元r、s、t结点累计应力向量和r、s、t结点主应力向量如下：

r结点：

第一主应力：

第二主应力：

第三主应力：

其中，r结点的累计应力向量用u方向累计应力v方向累计应力和垂直u、v平面的w方向累计应力表达为

s结点：

第一主应力：

第二主应力：

第三主应力：

其中，s结点的累计应力向量用u方向累计应力v方向累计应力和垂直u、v平面的w方向累计应力表达为

t结点：

第一主应力：

第二主应力：

第三主应力：

其中，t结点的累计应力向量用u方向累计应力v方向累计应力和垂直u、v平面的w方向累计应力表达为

(6-3-4)根据(6-3-3)得到的混凝土单元累计单元主应力求混凝土单元累计单元主应变

<mrow> <mo>{</mo> <msubsup> <mi>&amp;epsiv;</mi> <mi>M</mi> <mo>*</mo> </msubsup> <mo>}</mo> <mo>=</mo> <msup> <mrow> <mo>{</mo> <msup> <mrow> <mo>{</mo> <msubsup> <mi>&amp;epsiv;</mi> <mrow> <mi>M</mi> <mi>r</mi> </mrow> <mo>*</mo> </msubsup> <mo>}</mo> </mrow> <mi>T</mi> </msup> <mo>,</mo> <msup> <mrow> <mo>{</mo> <msubsup> <mi>&amp;epsiv;</mi> <mrow> <mi>M</mi> <mi>s</mi> </mrow> <mo>*</mo> </msubsup> <mo>}</mo> </mrow> <mi>T</mi> </msup> <mo>,</mo> <msup> <mrow> <mo>{</mo> <msubsup> <mi>&amp;epsiv;</mi> <mrow> <mi>M</mi> <mi>t</mi> </mrow> <mo>*</mo> </msubsup> <mo>}</mo> </mrow> <mi>T</mi> </msup> <mo>}</mo> </mrow> <mi>T</mi> </msup> </mrow>

其中，混凝土单元r结点的结点主应变记为 分别表示r结点第一、第二、第三主应变， s结点的结点主应变记为 分别表示s结点第一、第二、第三主应变，t结点的结点主应变记为 分别表示t结点第一、第二、第三主应变，E_c为混凝土的弹性模量；

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

<mrow> <msub> <mi>D</mi> <mrow> <mi>I</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>=</mo> <mfrac> <mn>1</mn> <mn>3</mn> </mfrac> <mo>&amp;CenterDot;</mo> <mrow> <mo>(</mo> <msubsup> <mi>D</mi> <mrow> <mi>I</mi> <mo>+</mo> <mn>1</mn> </mrow> <mi>r</mi> </msubsup> <mo>+</mo> <msubsup> <mi>D</mi> <mrow> <mi>I</mi> <mo>+</mo> <mn>1</mn> </mrow> <mi>s</mi> </msubsup> <mo>+</mo> <msubsup> <mi>D</mi> <mrow> <mi>I</mi> <mo>+</mo> <mn>1</mn> </mrow> <mi>t</mi> </msubsup> <mo>)</mo> </mrow> </mrow>

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

结点损伤因子是结点最大主应变ε的函数，其中r结点的损伤因子计算方法如下，由于只讨论

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

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

当大于0，小于0时，则单元处于复杂应力状态，分别取和代入以下情况3的ε值；

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

<mrow> <mi>&amp;sigma;</mi> <mo>=</mo> <msub> <mi>E</mi> <mn>0</mn> </msub> <mo>{</mo> <msub> <mi>&amp;epsiv;</mi> <mi>f</mi> </msub> <mrow> <mo>(</mo> <mn>1</mn> <mo>-</mo> <msub> <mi>A</mi> <mi>t</mi> </msub> <mo>)</mo> </mrow> <mo>+</mo> <mfrac> <mrow> <msub> <mi>A</mi> <mi>t</mi> </msub> <msub> <mi>&amp;epsiv;</mi> <mi>f</mi> </msub> </mrow> <mrow> <mi>exp</mi> <mo>&amp;lsqb;</mo> <msub> <mi>B</mi> <mi>t</mi> </msub> <mrow> <mo>(</mo> <mi>&amp;epsiv;</mi> <mo>-</mo> <msub> <mi>&amp;epsiv;</mi> <mi>f</mi> </msub> <mo>)</mo> </mrow> <mo>&amp;rsqb;</mo> </mrow> </mfrac> <mo>}</mo> </mrow>

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

<mrow> <msubsup> <mi>D</mi> <mrow> <mi>I</mi> <mo>+</mo> <mn>1</mn> </mrow> <mi>r</mi> </msubsup> <mo>=</mo> <mi>g</mi> <mrow> <mo>(</mo> <mi>&amp;epsiv;</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <mrow> <mi>&amp;epsiv;</mi> <mo>&amp;le;</mo> <msub> <mi>&amp;epsiv;</mi> <mi>f</mi> </msub> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mn>1</mn> <mo>-</mo> <mfrac> <mrow> <msub> <mi>&amp;epsiv;</mi> <mi>f</mi> </msub> <mrow> <mo>(</mo> <mn>1</mn> <mo>-</mo> <msub> <mi>A</mi> <mi>t</mi> </msub> <mo>)</mo> </mrow> </mrow> <mi>&amp;epsiv;</mi> </mfrac> <mo>-</mo> <mfrac> <msub> <mi>A</mi> <mi>t</mi> </msub> <mrow> <mi>exp</mi> <mo>&amp;lsqb;</mo> <msub> <mi>B</mi> <mi>t</mi> </msub> <mrow> <mo>(</mo> <mi>&amp;epsiv;</mi> <mo>-</mo> <msub> <mi>&amp;epsiv;</mi> <mi>f</mi> </msub> <mo>)</mo> </mrow> <mo>&amp;rsqb;</mo> </mrow> </mfrac> </mrow> </mtd> <mtd> <mrow> <mi>&amp;epsiv;</mi> <mo>&gt;</mo> <msub> <mi>&amp;epsiv;</mi> <mi>f</mi> </msub> </mrow> </mtd> </mtr> </mtable> </mfenced> </mrow>

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

<mrow> <mi>&amp;sigma;</mi> <mo>=</mo> <msub> <mi>E</mi> <mn>0</mn> </msub> <mo>{</mo> <mfrac> <mrow> <msub> <mi>&amp;epsiv;</mi> <mi>f</mi> </msub> <mrow> <mo>(</mo> <mn>1</mn> <mo>-</mo> <msub> <mi>A</mi> <mi>c</mi> </msub> <mo>)</mo> </mrow> </mrow> <mrow> <mo>-</mo> <msqrt> <mn>2</mn> </msqrt> <mi>v</mi> </mrow> </mfrac> <mo>+</mo> <mfrac> <mrow> <msub> <mi>A</mi> <mi>c</mi> </msub> <msub> <mi>&amp;epsiv;</mi> <mi>f</mi> </msub> </mrow> <mrow> <mi>exp</mi> <mo>&amp;lsqb;</mo> <msub> <mi>B</mi> <mi>c</mi> </msub> <mrow> <mo>(</mo> <mo>-</mo> <msqrt> <mn>2</mn> </msqrt> <mi>v</mi> <mi>&amp;epsiv;</mi> <mo>-</mo> <msub> <mi>&amp;epsiv;</mi> <mi>f</mi> </msub> <mo>)</mo> </mrow> <mo>&amp;rsqb;</mo> </mrow> </mfrac> <mo>}</mo> </mrow>

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

<mrow> <msubsup> <mi>D</mi> <mrow> <mi>I</mi> <mo>+</mo> <mn>1</mn> </mrow> <mi>r</mi> </msubsup> <mo>=</mo> <mi>g</mi> <mrow> <mo>(</mo> <mi>&amp;epsiv;</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <mrow> <mi>&amp;epsiv;</mi> <mo>&amp;le;</mo> <msub> <mi>&amp;epsiv;</mi> <mi>f</mi> </msub> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mn>1</mn> <mo>-</mo> <mfrac> <mrow> <msub> <mi>&amp;epsiv;</mi> <mi>f</mi> </msub> <mrow> <mo>(</mo> <mn>1</mn> <mo>-</mo> <msub> <mi>A</mi> <mi>c</mi> </msub> <mo>)</mo> </mrow> </mrow> <mi>&amp;epsiv;</mi> </mfrac> <mo>-</mo> <mfrac> <msub> <mi>A</mi> <mi>c</mi> </msub> <mrow> <mi>exp</mi> <mo>&amp;lsqb;</mo> <msub> <mi>B</mi> <mi>c</mi> </msub> <mrow> <mo>(</mo> <mi>&amp;epsiv;</mi> <mo>-</mo> <msub> <mi>&amp;epsiv;</mi> <mi>f</mi> </msub> <mo>)</mo> </mrow> <mo>&amp;rsqb;</mo> </mrow> </mfrac> </mrow> </mtd> <mtd> <mrow> <mi>&amp;epsiv;</mi> <mo>&gt;</mo> <msub> <mi>&amp;epsiv;</mi> <mi>f</mi> </msub> </mrow> </mtd> </mtr> </mtable> </mfenced> </mrow>

情况3：复杂应力状态情况下，拉、压损伤之间存在耦合效应，Mazars建议用下时描述单元总损伤

<mrow> <msubsup> <mi>D</mi> <mrow> <mi>I</mi> <mo>+</mo> <mn>1</mn> </mrow> <mi>r</mi> </msubsup> <mo>=</mo> <msub> <mi>&amp;alpha;</mi> <mi>T</mi> </msub> <msub> <mi>D</mi> <mi>T</mi> </msub> <mo>+</mo> <msub> <mi>&amp;alpha;</mi> <mi>C</mi> </msub> <msub> <mi>D</mi> <mi>C</mi> </msub> </mrow>

<mrow> <msub> <mi>D</mi> <mi>T</mi> </msub> <mo>=</mo> <mi>g</mi> <mrow> <mo>(</mo> <mi>&amp;epsiv;</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <mrow> <mi>&amp;epsiv;</mi> <mo>&amp;le;</mo> <msub> <mi>&amp;epsiv;</mi> <mi>f</mi> </msub> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mn>1</mn> <mo>-</mo> <mfrac> <mrow> <msub> <mi>&amp;epsiv;</mi> <mi>f</mi> </msub> <mrow> <mo>(</mo> <mn>1</mn> <mo>-</mo> <msub> <mi>A</mi> <mi>t</mi> </msub> <mo>)</mo> </mrow> </mrow> <mi>&amp;epsiv;</mi> </mfrac> <mo>-</mo> <mfrac> <msub> <mi>A</mi> <mi>t</mi> </msub> <mrow> <mi>exp</mi> <mo>&amp;lsqb;</mo> <msub> <mi>B</mi> <mi>t</mi> </msub> <mrow> <mo>(</mo> <mi>&amp;epsiv;</mi> <mo>-</mo> <msub> <mi>&amp;epsiv;</mi> <mi>f</mi> </msub> <mo>)</mo> </mrow> <mo>&amp;rsqb;</mo> </mrow> </mfrac> </mrow> </mtd> <mtd> <mrow> <mi>&amp;epsiv;</mi> <mo>&gt;</mo> <msub> <mi>&amp;epsiv;</mi> <mi>f</mi> </msub> </mrow> </mtd> </mtr> </mtable> </mfenced> </mrow>

<mrow> <msub> <mi>D</mi> <mi>C</mi> </msub> <mo>=</mo> <mi>g</mi> <mrow> <mo>(</mo> <mi>&amp;epsiv;</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <mrow> <mi>&amp;epsiv;</mi> <mo>&amp;le;</mo> <msub> <mi>&amp;epsiv;</mi> <mi>f</mi> </msub> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mn>1</mn> <mo>-</mo> <mfrac> <mrow> <msub> <mi>&amp;epsiv;</mi> <mi>f</mi> </msub> <mrow> <mo>(</mo> <mn>1</mn> <mo>-</mo> <msub> <mi>A</mi> <mi>c</mi> </msub> <mo>)</mo> </mrow> </mrow> <mi>&amp;epsiv;</mi> </mfrac> <mo>-</mo> <mfrac> <msub> <mi>A</mi> <mi>c</mi> </msub> <mrow> <mi>exp</mi> <mo>&amp;lsqb;</mo> <msub> <mi>B</mi> <mi>c</mi> </msub> <mrow> <mo>(</mo> <mi>&amp;epsiv;</mi> <mo>-</mo> <msub> <mi>&amp;epsiv;</mi> <mi>f</mi> </msub> <mo>)</mo> </mrow> <mo>&amp;rsqb;</mo> </mrow> </mfrac> </mrow> </mtd> <mtd> <mrow> <mi>&amp;epsiv;</mi> <mo>&gt;</mo> <msub> <mi>&amp;epsiv;</mi> <mi>f</mi> </msub> </mrow> </mtd> </mtr> </mtable> </mfenced> </mrow>

其中，α_T和α_C为耦合系数，纯压缩时取α_T＝0，纯拉伸时取α_C＝0，组合情况α_T+α_C＝1，D_T和D_C为耦合损伤因子；

采用相同方法计算s、t结点得损伤因子

(6-3-6)求解第I+1步的局部坐标和整体坐标系下混凝土单元刚度矩阵：

第I+1步的局部坐系下标混凝土单元刚度矩阵计算方法如下：

<mrow> <mo>&amp;lsqb;</mo> <msubsup> <mi>k</mi> <mrow> <mi>I</mi> <mo>+</mo> <mn>1</mn> </mrow> <mi>c</mi> </msubsup> <mo>&amp;rsqb;</mo> <mo>=</mo> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <mrow> <mo>&amp;lsqb;</mo> <msubsup> <mi>k</mi> <mrow> <mi>I</mi> <mo>+</mo> <mn>1</mn> </mrow> <mi>p</mi> </msubsup> <mo>&amp;rsqb;</mo> </mrow> </mtd> <mtd> <mrow></mrow> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> </mtr> <mtr> <mtd> <mrow></mrow> </mtd> <mtd> <mrow></mrow> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mrow></mrow> </mtd> <mtd> <mrow></mrow> </mtd> <mtd> <mrow></mrow> </mtd> <mtd> <mn>0</mn> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mrow></mrow> </mtd> <mtd> <mrow> <mo>&amp;lsqb;</mo> <msubsup> <mi>k</mi> <mrow> <mi>I</mi> <mo>+</mo> <mn>1</mn> </mrow> <mi>b</mi> </msubsup> <mo>&amp;rsqb;</mo> </mrow> </mtd> <mtd> <mrow></mrow> </mtd> <mtd> <mn>0</mn> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mrow></mrow> </mtd> <mtd> <mrow></mrow> </mtd> <mtd> <mrow></mrow> </mtd> <mtd> <mn>0</mn> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> </mtr> </mtable> </mfenced> </mrow>

其中，

第I+1步的整体坐标系的混凝土单元刚度矩阵形成方法同上述步骤4所述；

(6-4)由局部坐标系下的计算钢筋平面应力单元的应力向量根据计算前I荷载步的钢筋平面应力单元累计应力向量并由求解前I荷载步累计主应力向量和累计主应变向量由累计主应变向量判断钢筋是否屈服并求解第I+1荷载步钢筋单元刚度矩阵，具体如下：

(6-4-1)由局部坐标系下的计算钢筋单元平面应力单元的应力向量钢筋平面应力单元的应力向量为：

<mrow> <mo>{</mo> <msubsup> <mover> <mi>&amp;sigma;</mi> <mo>&amp;OverBar;</mo> </mover> <mi>I</mi> <mi>p</mi> </msubsup> <mo>}</mo> <mo>=</mo> <msup> <mrow> <mo>{</mo> <msup> <mrow> <mo>{</mo> <msubsup> <mover> <mi>&amp;sigma;</mi> <mo>&amp;OverBar;</mo> </mover> <mrow> <mi>I</mi> <mi>r</mi> </mrow> <mi>p</mi> </msubsup> <mo>}</mo> </mrow> <mi>T</mi> </msup> <mo>,</mo> <msup> <mrow> <mo>{</mo> <msubsup> <mover> <mi>&amp;sigma;</mi> <mo>&amp;OverBar;</mo> </mover> <mrow> <mi>I</mi> <mi>s</mi> </mrow> <mi>p</mi> </msubsup> <mo>}</mo> </mrow> <mi>T</mi> </msup> <mo>,</mo> <msup> <mrow> <mo>{</mo> <msubsup> <mover> <mi>&amp;sigma;</mi> <mo>&amp;OverBar;</mo> </mover> <mrow> <mi>I</mi> <mi>t</mi> </mrow> <mi>p</mi> </msubsup> <mo>}</mo> </mrow> <mi>T</mi> </msup> <mo>}</mo> </mrow> <mi>T</mi> </msup> </mrow>

其中， 表示钢筋平面应力单元的弹性矩阵，表示钢筋平面应力单元的应变矩阵；

(6-4-2)根据(6-4-1)得到的第I荷载步钢筋单元平面应力单元的应力向量计算前I荷载步钢筋单元的累计应力向量具体计算方法如下：

计算前I荷载步钢筋单元的累计应力向量

<mrow> <mo>{</mo> <msubsup> <mover> <mi>&amp;sigma;</mi> <mo>&amp;OverBar;</mo> </mover> <mi>I</mi> <mo>*</mo> </msubsup> <mo>}</mo> <mo>=</mo> <msup> <mrow> <mo>{</mo> <msup> <mrow> <mo>{</mo> <msubsup> <mover> <mi>&amp;sigma;</mi> <mo>&amp;OverBar;</mo> </mover> <mrow> <mi>I</mi> <mi>r</mi> </mrow> <mo>*</mo> </msubsup> <mo>}</mo> </mrow> <mi>T</mi> </msup> <mo>,</mo> <msup> <mrow> <mo>{</mo> <msubsup> <mover> <mi>&amp;sigma;</mi> <mo>&amp;OverBar;</mo> </mover> <mrow> <mi>I</mi> <mi>s</mi> </mrow> <mo>*</mo> </msubsup> <mo>}</mo> </mrow> <mi>T</mi> </msup> <mo>,</mo> <msup> <mrow> <mo>{</mo> <msubsup> <mover> <mi>&amp;sigma;</mi> <mo>&amp;OverBar;</mo> </mover> <mrow> <mi>I</mi> <mi>t</mi> </mrow> <mo>*</mo> </msubsup> <mo>}</mo> </mrow> <mi>T</mi> </msup> <mo>}</mo> </mrow> <mi>T</mi> </msup> </mrow>

其中，钢筋单元r结点的累计应力向量s结点的累计应力向量t结点的累计应力向量

(6-4-3)根据(6-4-2)得到的前I荷载步钢筋单元的累计应力向量求前I荷载步钢筋单元累计主应力向量

<mrow> <msup> <mrow> <mo>{</mo> <msubsup> <mover> <mi>&amp;sigma;</mi> <mo>&amp;OverBar;</mo> </mover> <mi>M</mi> <mo>*</mo> </msubsup> <mo>}</mo> </mrow> <mi>T</mi> </msup> <mo>=</mo> <mo>{</mo> <msup> <mrow> <mo>{</mo> <msubsup> <mover> <mi>&amp;sigma;</mi> <mo>&amp;OverBar;</mo> </mover> <mrow> <mi>M</mi> <mi>r</mi> </mrow> <mo>*</mo> </msubsup> <mo>}</mo> </mrow> <mi>T</mi> </msup> <mo>,</mo> <msup> <mrow> <mo>{</mo> <msubsup> <mover> <mi>&amp;sigma;</mi> <mo>&amp;OverBar;</mo> </mover> <mrow> <mi>M</mi> <mi>s</mi> </mrow> <mo>*</mo> </msubsup> <mo>}</mo> </mrow> <mi>T</mi> </msup> <mo>,</mo> <msup> <mrow> <mo>{</mo> <msubsup> <mover> <mi>&amp;sigma;</mi> <mo>&amp;OverBar;</mo> </mover> <mrow> <mi>M</mi> <mi>t</mi> </mrow> <mo>*</mo> </msubsup> <mo>}</mo> </mrow> <mi>T</mi> </msup> <mo>}</mo> </mrow>

前I荷载步钢筋单元r、s、t结点累计应力向量和r、s、t结点主应力向量如下：

r结点：

第一主应力：

第二主应力：

第三主应力：

其中，r结点的累计应力向量用u方向累计应力v方向累计应力和垂直u、v平面的w方向累计应力表达为

s结点：

第一主应力：

第二主应力：

第三主应力：

其中，s结点的累计应力向量用u方向累计应力v方向累计应力和垂直u、v平面的w方向累计应力表达为

t结点：

第一主应力：

第二主应力：

第三主应力：

其中，t结点的累计应力向量用u方向累计应力v方向累计应力和垂直u、v平面的w方向累计应力表达为

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

<mrow> <mo>{</mo> <msubsup> <mover> <mi>&amp;epsiv;</mi> <mo>&amp;OverBar;</mo> </mover> <mi>M</mi> <mo>*</mo> </msubsup> <mo>}</mo> <mo>=</mo> <msup> <mrow> <mo>{</mo> <msup> <mrow> <mo>{</mo> <msubsup> <mover> <mi>&amp;epsiv;</mi> <mo>&amp;OverBar;</mo> </mover> <mrow> <mi>M</mi> <mi>r</mi> </mrow> <mo>*</mo> </msubsup> <mo>}</mo> </mrow> <mi>T</mi> </msup> <mo>,</mo> <msup> <mrow> <mo>{</mo> <msubsup> <mover> <mi>&amp;epsiv;</mi> <mo>&amp;OverBar;</mo> </mover> <mrow> <mi>M</mi> <mi>s</mi> </mrow> <mo>*</mo> </msubsup> <mo>}</mo> </mrow> <mi>T</mi> </msup> <mo>,</mo> <msup> <mrow> <mo>{</mo> <msubsup> <mover> <mi>&amp;epsiv;</mi> <mo>&amp;OverBar;</mo> </mover> <mrow> <mi>M</mi> <mi>t</mi> </mrow> <mo>*</mo> </msubsup> <mo>}</mo> </mrow> <mi>T</mi> </msup> <mo>}</mo> </mrow> <mi>T</mi> </msup> </mrow>

其中，钢筋单元r结点的结点主应变记为向量 分别表示r结点第一、第二、第三主应变， s结点的结点主应变记为向量分别表示s结点第一、第二、第三主应变，t结点的结点主应变记为向量分别表示t结点第一、第二、第三主应变，E_s为钢筋的弹性模量；

(6-4-5)求解第I+1步钢筋单元的弹性模量和刚度矩阵：

取钢筋单元主累计应变中绝对值最大值为钢筋单元最大主应变为ε_s；

若ε_s≤f_y/E_s，则钢筋单元未屈服，第I+1步钢筋单元的弹性模量为E_s，钢筋单元的刚度矩阵为：

<mrow> <mo>&amp;lsqb;</mo> <msubsup> <mi>k</mi> <mrow> <mi>I</mi> <mo>+</mo> <mn>1</mn> </mrow> <mi>s</mi> </msubsup> <mo>&amp;rsqb;</mo> <mo>=</mo> <mo>&amp;Integral;</mo> <mo>&amp;Integral;</mo> <mo>&amp;lsqb;</mo> <msup> <mover> <mi>B</mi> <mo>&amp;OverBar;</mo> </mover> <mi>p</mi> </msup> <mo>&amp;rsqb;</mo> <mo>&amp;lsqb;</mo> <msup> <mover> <mi>D</mi> <mo>&amp;OverBar;</mo> </mover> <mi>p</mi> </msup> <mo>&amp;rsqb;</mo> <mo>&amp;lsqb;</mo> <msup> <mover> <mi>B</mi> <mo>&amp;OverBar;</mo> </mover> <mi>p</mi> </msup> <mo>&amp;rsqb;</mo> <mi>d</mi> <mi>A</mi> </mrow>

<mrow> <mo>&amp;lsqb;</mo> <msup> <mover> <mi>D</mi> <mo>&amp;OverBar;</mo> </mover> <mi>p</mi> </msup> <mo>&amp;rsqb;</mo> <mo>=</mo> <mfrac> <msub> <mi>E</mi> <mi>s</mi> </msub> <mrow> <mn>1</mn> <mo>-</mo> <msup> <mi>v</mi> <mn>2</mn> </msup> </mrow> </mfrac> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <mn>1</mn> </mtd> <mtd> <mi>v</mi> </mtd> <mtd> <mn>0</mn> </mtd> </mtr> <mtr> <mtd> <mi>v</mi> </mtd> <mtd> <mn>1</mn> </mtd> <mtd> <mn>0</mn> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mfrac> <mrow> <mn>1</mn> <mo>-</mo> <mi>v</mi> </mrow> <mn>2</mn> </mfrac> </mtd> </mtr> </mtable> </mfenced> </mrow>

其中，表示钢筋平面应力单元的弹性矩阵；

若ε_s＞f_y/E_s，则钢筋单元屈服，第I+1步钢筋单元的弹性模量为E_s2，钢筋单元的刚度矩阵为：

<mrow> <mo>&amp;lsqb;</mo> <msubsup> <mi>k</mi> <mrow> <mi>I</mi> <mo>+</mo> <mn>1</mn> </mrow> <mi>s</mi> </msubsup> <mo>&amp;rsqb;</mo> <mo>=</mo> <mo>&amp;Integral;</mo> <mo>&amp;Integral;</mo> <mo>&amp;lsqb;</mo> <msup> <mover> <mi>B</mi> <mo>&amp;OverBar;</mo> </mover> <mi>p</mi> </msup> <mo>&amp;rsqb;</mo> <mo>&amp;lsqb;</mo> <msup> <mover> <mi>D</mi> <mo>&amp;OverBar;</mo> </mover> <mi>p</mi> </msup> <mo>&amp;rsqb;</mo> <mo>&amp;lsqb;</mo> <msup> <mover> <mi>B</mi> <mo>&amp;OverBar;</mo> </mover> <mi>p</mi> </msup> <mo>&amp;rsqb;</mo> <mi>d</mi> <mi>A</mi> </mrow>

<mrow> <mo>&amp;lsqb;</mo> <msup> <mover> <mi>D</mi> <mo>&amp;OverBar;</mo> </mover> <mi>p</mi> </msup> <mo>&amp;rsqb;</mo> <mo>=</mo> <mfrac> <msub> <mi>E</mi> <mrow> <mi>s</mi> <mn>2</mn> </mrow> </msub> <mrow> <mn>1</mn> <mo>-</mo> <msup> <mi>v</mi> <mn>2</mn> </msup> </mrow> </mfrac> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <mn>1</mn> </mtd> <mtd> <mi>v</mi> </mtd> <mtd> <mn>0</mn> </mtd> </mtr> <mtr> <mtd> <mi>v</mi> </mtd> <mtd> <mn>1</mn> </mtd> <mtd> <mn>0</mn> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mfrac> <mrow> <mn>1</mn> <mo>-</mo> <mi>v</mi> </mrow> <mn>2</mn> </mfrac> </mtd> </mtr> </mtable> </mfenced> </mrow>

其中，表示钢筋平面应力单元的弹性矩阵；

(6-5)求解第I+1步钢筋混凝土单元刚度矩阵[k_I+1]：

<mrow> <mo>&amp;lsqb;</mo> <msub> <mi>k</mi> <mrow> <mi>I</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>&amp;rsqb;</mo> <mo>=</mo> <mo>&amp;lsqb;</mo> <msubsup> <mi>k</mi> <mrow> <mi>I</mi> <mo>+</mo> <mn>1</mn> </mrow> <mi>c</mi> </msubsup> <mo>&amp;rsqb;</mo> <mo>+</mo> <mo>&amp;lsqb;</mo> <msubsup> <mi>k</mi> <mrow> <mi>I</mi> <mo>+</mo> <mn>1</mn> </mrow> <mi>s</mi> </msubsup> <mo>&amp;rsqb;</mo> </mrow>

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

在整体坐标下，记第J迭代步的位移误差向量为且J∈[1，M]；M为迭代次数，在程序运行时设定；

对于第I荷载步第1迭代步的公式如下：

<mrow> <mo>&amp;lsqb;</mo> <msub> <mi>k</mi> <mrow> <mi>I</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>&amp;rsqb;</mo> <mo>{</mo> <msubsup> <mover> <mi>&amp;delta;</mi> <mo>&amp;OverBar;</mo> </mover> <mn>1</mn> <mo>&amp;prime;</mo> </msubsup> <mo>}</mo> <mo>=</mo> <mfrac> <mn>1</mn> <mi>N</mi> </mfrac> <mo>{</mo> <mi>p</mi> <mo>}</mo> <mo>-</mo> <mo>&amp;lsqb;</mo> <msub> <mi>k</mi> <mrow> <mi>I</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>&amp;rsqb;</mo> <mo>{</mo> <msubsup> <mi>&amp;delta;</mi> <mi>I</mi> <mo>&amp;prime;</mo> </msubsup> <mo>}</mo> </mrow>

可以求解{δ′₁}，对于第J迭代步公式如下：

<mrow> <mo>&amp;lsqb;</mo> <msub> <mi>k</mi> <mrow> <mi>I</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>&amp;rsqb;</mo> <mo>{</mo> <msubsup> <mover> <mi>&amp;delta;</mi> <mo>&amp;OverBar;</mo> </mover> <mi>j</mi> <mo>&amp;prime;</mo> </msubsup> <mo>}</mo> <mo>=</mo> <mfrac> <mn>1</mn> <mi>N</mi> </mfrac> <mo>{</mo> <mi>p</mi> <mo>}</mo> <mo>-</mo> <mo>&amp;lsqb;</mo> <msub> <mi>k</mi> <mrow> <mi>I</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>&amp;rsqb;</mo> <mrow> <mo>(</mo> <mo>{</mo> <msubsup> <mi>&amp;delta;</mi> <mi>I</mi> <mo>&amp;prime;</mo> </msubsup> <mo>}</mo> <mo>+</mo> <mo>{</mo> <msubsup> <mover> <mi>&amp;delta;</mi> <mo>&amp;OverBar;</mo> </mover> <mn>1</mn> <mo>&amp;prime;</mo> </msubsup> <mo>}</mo> <mo>+</mo> <mo>{</mo> <msubsup> <mover> <mi>&amp;delta;</mi> <mo>&amp;OverBar;</mo> </mover> <mn>2</mn> <mo>&amp;prime;</mo> </msubsup> <mo>}</mo> <mo>+</mo> <mn>...</mn> <mo>+</mo> <mo>{</mo> <msubsup> <mover> <mi>&amp;delta;</mi> <mo>&amp;OverBar;</mo> </mover> <mrow> <mi>j</mi> <mo>-</mo> <mn>1</mn> </mrow> <mo>&amp;prime;</mo> </msubsup> <mo>}</mo> <mo>)</mo> </mrow> </mrow>

对于第M迭代步的迭代公式如下：

<mrow> <mo>&amp;lsqb;</mo> <msub> <mi>k</mi> <mrow> <mi>I</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>&amp;rsqb;</mo> <mo>{</mo> <msubsup> <mover> <mi>&amp;delta;</mi> <mo>&amp;OverBar;</mo> </mover> <mi>M</mi> <mo>&amp;prime;</mo> </msubsup> <mo>}</mo> <mo>=</mo> <mfrac> <mn>1</mn> <mi>N</mi> </mfrac> <mo>{</mo> <mi>p</mi> <mo>}</mo> <mo>-</mo> <mo>&amp;lsqb;</mo> <msub> <mi>k</mi> <mrow> <mi>I</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>&amp;rsqb;</mo> <mrow> <mo>(</mo> <mo>{</mo> <msubsup> <mi>&amp;delta;</mi> <mi>I</mi> <mo>&amp;prime;</mo> </msubsup> <mo>}</mo> <mo>+</mo> <mo>{</mo> <msubsup> <mover> <mi>&amp;delta;</mi> <mo>&amp;OverBar;</mo> </mover> <mn>1</mn> <mo>&amp;prime;</mo> </msubsup> <mo>}</mo> <mo>+</mo> <mo>{</mo> <msubsup> <mover> <mi>&amp;delta;</mi> <mo>&amp;OverBar;</mo> </mover> <mn>2</mn> <mo>&amp;prime;</mo> </msubsup> <mo>}</mo> <mo>+</mo> <mn>...</mn> <mo>+</mo> <mo>{</mo> <msubsup> <mover> <mi>&amp;delta;</mi> <mo>&amp;OverBar;</mo> </mover> <mrow> <mi>M</mi> <mo>-</mo> <mn>1</mn> </mrow> <mo>&amp;prime;</mo> </msubsup> <mo>}</mo> <mo>)</mo> </mrow> </mrow>

(6-7)由上述步骤可以得到同时考虑Mazars损伤模型和钢筋作用的{p}荷载作用下的第I荷载步的结点位移，记为{W_I}：

<mrow> <mo>{</mo> <msub> <mi>W</mi> <mi>I</mi> </msub> <mo>}</mo> <mo>=</mo> <mo>{</mo> <msubsup> <mi>&amp;delta;</mi> <mi>I</mi> <mo>&amp;prime;</mo> </msubsup> <mo>}</mo> <mo>+</mo> <mo>{</mo> <msubsup> <mover> <mi>&amp;delta;</mi> <mo>&amp;OverBar;</mo> </mover> <mn>1</mn> <mo>&amp;prime;</mo> </msubsup> <mo>}</mo> <mo>+</mo> <mo>{</mo> <msubsup> <mover> <mi>&amp;delta;</mi> <mo>&amp;OverBar;</mo> </mover> <mn>2</mn> <mo>&amp;prime;</mo> </msubsup> <mo>}</mo> <mo>+</mo> <mo>...</mo> <mo>+</mo> <mo>{</mo> <msubsup> <mover> <mi>&amp;delta;</mi> <mo>&amp;OverBar;</mo> </mover> <mi>M</mi> <mo>&amp;prime;</mo> </msubsup> <mo>}</mo> </mrow>

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

最终将N个荷载步的所有结点位移累加，就可以的得到同时考虑Mazars损伤模型和钢筋作用的{p}荷载作用下的结点位移

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

7.根据权利要求1所述的一种考虑材料损伤的自由曲面形态创建方法，其特征在于，步骤8所述的对考虑Mazars损伤模型和钢筋作用的自由曲面控制点高度调整，具体为：

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

<mrow> <msup> <mi>&amp;Delta;z</mi> <mrow> <mo>(</mo> <mi>L</mi> <mo>)</mo> </mrow> </msup> <mo>=</mo> <msup> <mi>&amp;lambda;</mi> <mrow> <mo>(</mo> <mi>L</mi> <mo>)</mo> </mrow> </msup> <mo>&amp;CenterDot;</mo> <munderover> <mi>&amp;Sigma;</mi> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mrow> <mi>m</mi> <mo>&amp;times;</mo> <mi>n</mi> </mrow> </munderover> <mfrac> <mrow> <mo>&amp;part;</mo> <mi>C</mi> <mrow> <mo>(</mo> <msup> <mi>z</mi> <mrow> <mo>(</mo> <mi>L</mi> <mo>)</mo> </mrow> </msup> <mo>)</mo> </mrow> </mrow> <mrow> <mo>&amp;part;</mo> <msub> <mi>z</mi> <mi>i</mi> </msub> </mrow> </mfrac> <msub> <mi>&amp;Delta;z</mi> <mi>i</mi> </msub> </mrow>

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

<mrow> <msubsup> <mi>z</mi> <mi>i</mi> <mrow> <mo>(</mo> <mi>L</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> </msubsup> <mo>=</mo> <msubsup> <mi>z</mi> <mi>i</mi> <mrow> <mo>(</mo> <mi>L</mi> <mo>)</mo> </mrow> </msubsup> <mo>+</mo> <msup> <mi>&amp;Delta;z</mi> <mrow> <mo>(</mo> <mi>L</mi> <mo>)</mo> </mrow> </msup> <mo>=</mo> <msubsup> <mi>z</mi> <mi>i</mi> <mrow> <mo>(</mo> <mi>L</mi> <mo>)</mo> </mrow> </msubsup> <mo>+</mo> <msup> <mi>&amp;lambda;</mi> <mrow> <mo>(</mo> <mi>L</mi> <mo>)</mo> </mrow> </msup> <munderover> <mi>&amp;Sigma;</mi> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mrow> <mi>m</mi> <mo>&amp;times;</mo> <mi>n</mi> </mrow> </munderover> <mfrac> <mrow> <mo>&amp;part;</mo> <mi>C</mi> <mrow> <mo>(</mo> <msup> <mi>z</mi> <mrow> <mo>(</mo> <mi>L</mi> <mo>)</mo> </mrow> </msup> <mo>)</mo> </mrow> </mrow> <mrow> <mo>&amp;part;</mo> <msub> <mi>z</mi> <mi>i</mi> </msub> </mrow> </mfrac> <msub> <mi>&amp;Delta;z</mi> <mi>i</mi> </msub> </mrow>

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

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

结构应变能C是z^(L)的函数，假定应变能C是关于z^(L)的连续可微函数，将C(z^(L))进行Taylor展开：

<mrow> <mi>C</mi> <mrow> <mo>(</mo> <msup> <mi>z</mi> <mrow> <mo>(</mo> <mi>L</mi> <mo>)</mo> </mrow> </msup> <mo>+</mo> <msub> <mi>&amp;Delta;z</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mi>C</mi> <mrow> <mo>(</mo> <msup> <mi>z</mi> <mrow> <mo>(</mo> <mi>L</mi> <mo>)</mo> </mrow> </msup> <mo>)</mo> </mrow> <mo>+</mo> <munderover> <mi>&amp;Sigma;</mi> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mrow> <mi>m</mi> <mo>&amp;times;</mo> <mi>n</mi> </mrow> </munderover> <mfrac> <mrow> <mo>&amp;part;</mo> <mi>C</mi> <mrow> <mo>(</mo> <msup> <mi>z</mi> <mrow> <mo>(</mo> <mi>L</mi> <mo>)</mo> </mrow> </msup> <mo>)</mo> </mrow> </mrow> <mrow> <mo>&amp;part;</mo> <msub> <mi>z</mi> <mi>i</mi> </msub> </mrow> </mfrac> <msub> <mi>&amp;Delta;z</mi> <mi>i</mi> </msub> <mo>+</mo> <mi>o</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;Delta;z</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> </mrow>

<mrow> <mi>C</mi> <mrow> <mo>(</mo> <msup> <mi>z</mi> <mrow> <mo>(</mo> <mi>L</mi> <mo>)</mo> </mrow> </msup> <mo>-</mo> <msub> <mi>&amp;Delta;z</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mi>C</mi> <mrow> <mo>(</mo> <msup> <mi>z</mi> <mrow> <mo>(</mo> <mi>L</mi> <mo>)</mo> </mrow> </msup> <mo>)</mo> </mrow> <mo>-</mo> <munderover> <mi>&amp;Sigma;</mi> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mrow> <mi>m</mi> <mo>&amp;times;</mo> <mi>n</mi> </mrow> </munderover> <mfrac> <mrow> <mo>&amp;part;</mo> <mi>C</mi> <mrow> <mo>(</mo> <msup> <mi>z</mi> <mrow> <mo>(</mo> <mi>L</mi> <mo>)</mo> </mrow> </msup> <mo>)</mo> </mrow> </mrow> <mrow> <mo>&amp;part;</mo> <msub> <mi>z</mi> <mi>i</mi> </msub> </mrow> </mfrac> <msub> <mi>&amp;Delta;z</mi> <mi>i</mi> </msub> <mo>+</mo> <mi>o</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;Delta;z</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> </mrow>

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

<mrow> <munderover> <mi>&amp;Sigma;</mi> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mrow> <mi>m</mi> <mo>&amp;times;</mo> <mi>n</mi> </mrow> </munderover> <mfrac> <mrow> <mo>&amp;part;</mo> <mi>C</mi> <mrow> <mo>(</mo> <msup> <mi>z</mi> <mrow> <mo>(</mo> <mi>L</mi> <mo>)</mo> </mrow> </msup> <mo>)</mo> </mrow> </mrow> <mrow> <mo>&amp;part;</mo> <msub> <mi>z</mi> <mi>i</mi> </msub> </mrow> </mfrac> <msub> <mi>&amp;Delta;z</mi> <mi>i</mi> </msub> <mo>=</mo> <mfrac> <mrow> <mi>C</mi> <mrow> <mo>(</mo> <msup> <mi>z</mi> <mrow> <mo>(</mo> <mi>L</mi> <mo>)</mo> </mrow> </msup> <mo>+</mo> <msub> <mi>&amp;Delta;z</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> <mo>-</mo> <mi>C</mi> <mrow> <mo>(</mo> <msup> <mi>z</mi> <mrow> <mo>(</mo> <mi>L</mi> <mo>)</mo> </mrow> </msup> <mo>-</mo> <msub> <mi>&amp;Delta;z</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> </mrow> <mrow> <mn>2</mn> <msub> <mi>&amp;Delta;z</mi> <mi>i</mi> </msub> </mrow> </mfrac> <mo>-</mo> <mfrac> <msup> <mrow> <mo>(</mo> <msub> <mi>&amp;Delta;z</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mn>6</mn> </mfrac> <msup> <mi>C</mi> <mrow> <mo>(</mo> <mn>3</mn> <mo>)</mo> </mrow> </msup> <mrow> <mo>(</mo> <mi>&amp;xi;</mi> <mo>)</mo> </mrow> </mrow>

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

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

其计算步骤如下，其中F表示黄金分割法的第F迭代步，且F属于正整数集，λ_F为第F次迭代求得的步长值，当程序循环至满足step2的条件时将λ_F的值赋给λ^(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。