本发明涉及一种考虑材料损伤和钢筋作用的自由曲面结构形态创建方法。
背景技术:
国内外目前自由曲面结构的形态创建方法可以概括为以下两种。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];pi,j是u方向编号为i且v方向编号为j的控制点坐标;wi,j是u方向编号为i且v方向编号为j的控制点的权因子;ni,p(u)和nj,q(v)分别是定义在矢量u和v上的样条基函数,其表达式如下所示,ui为u方向的结点矢量且ui∈u,
所述步骤3计算步骤2得到的钢筋混凝土三角形壳单元的局部坐标系下的刚度矩阵具体步骤为:
钢筋混凝土三角形单元分为混凝土单元和钢筋单元,混凝土单元和钢筋单元共结点,两者刚度矩阵叠加得到钢筋混凝土三角形单元的刚度矩阵。混凝土单元采用三角形薄壳单元模拟,单元含3个结点,每个结点含6个自由度。钢筋单元采用三角形平面应力单元模拟,单元含3个结点,每个结点含2个自由度。
(3-1)混凝土单元采用三角形薄壳单元模拟,三角形薄壳单元的刚度矩阵由平面应力状态和弯曲应力状态的刚度矩阵得到,三角形薄壳单元的三个结点分别记为r、s、t。
在混凝土三角形平面应力单元中,结点力和结点位移的关系如下:
其中,
在混凝土三角形弯曲应力单元中,结点力和结点位移的关系如下:
其中,
混凝土三角形薄壳单元的局部坐标系下的刚度矩阵[kc]为:
(3-2)钢筋单元采用三角形平面应力单元模拟,在钢筋三角形平面应力单元中,结点力和结点位移的关系如下:
其中,
(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属于正整数集;根据[ki]{δ′i}={p}/n,求第i步的结点位移{δ′i};当i=1时,混凝土初始损伤d1为0,钢筋不发生屈服,[k1]=[k],[k]为步骤5中得到的自由曲面结构的整体刚度矩阵。
(6-2)通过坐标转换,将{δ′i}转换为局部坐标系下的结点位移向量{δi},{δi}=[l]{δ′i};并将{δi}分解为平面应力单元对应的位移和弯曲应力单元对应的位移
(6-3)由局部坐标系下的
(6-3-1)由局部坐标系下的
混凝土平面应力单元的应力向量为:
其中,
混凝土弯曲应力单元的应力向量为:
其中,
第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结点的累计应力向量
(6-3-3)根据(6-3-2)得到的前i荷载步混凝土单元的累计应力向量
前i荷载步混凝土单元r、s、t结点累计应力向量
r结点:
第一主应力:
第二主应力:
第三主应力:
其中,r结点的累计应力向量
s结点:
第一主应力:
第二主应力:
第三主应力:
其中,s结点的累计应力向量
t结点:
第一主应力:
第二主应力:
第三主应力:
其中,t结点的累计应力向量
(6-3-4)根据(6-3-3)得到的混凝土单元累计单元主应力
其中,混凝土单元r结点的结点主应变记为
(6-3-5)求解第i+1步的损伤因子di+1:
其中,
结点损伤因子是结点最大主应变ε的函数,其中r结点的损伤因子
当
当
当
情况1:拉伸情况下,峰值应力之后其应力-应变曲线可以用下式来描述,式中at和bt为拉伸时的材料参数,εf为峰值应力的应变,对于混凝土材料取0.5<at<1.0,104<bt<105。
由上式可以得到拉伸状态下损伤演化方程,则
情况2:压缩情况下,峰值应力之后其应力-应变曲线可以用下式来描述,式中ac和bc为压缩时的材料参数,εf为峰值应力的应变,对于混凝土材料取1.0<ac<1.5,1000<bc<2000。
由上式可以得到单轴压缩损伤演化方程,则
情况3:复杂应力状态情况下,拉、压损伤之间存在耦合效应,mazars建议用下时描述单元总损伤
其中,αt和αc为耦合系数,纯压缩时取αt=0,纯拉伸时取αc=0,组合情况αt+αc=1,dt和dc为耦合损伤因子;
采用相同方法计算s、t结点得损伤因子
(6-3-6)求解第i+1步的局部坐标和整体坐标系下混凝土单元刚度矩阵:
第i+1步的局部坐系下标混凝土单元刚度矩阵计算方法如下:
其中,
第i+1步的整体坐标系的混凝土单元刚度矩阵形成方法同上述步骤4所述。
(6-4)由局部坐标系下的
(6-4-1)由局部坐标系下的
其中,
(6-4-2)根据(6-4-1)得到的第i荷载步钢筋单元平面应力单元的应力向量
计算前i荷载步钢筋单元的累计应力向量
其中,钢筋单元r结点的累计应力向量
(6-4-3)根据(6-4-2)得到的前i荷载步钢筋单元的累计应力向量
前i荷载步钢筋单元r、s、t结点累计应力向量
r结点:
第一主应力:
第二主应力:
第三主应力:
其中,r结点的累计应力向量
s结点:
第一主应力:
第二主应力:
第三主应力:
其中,s结点的累计应力向量
t结点:
第一主应力:
第二主应力:
第三主应力:
其中,t结点的累计应力向量
(6-4-4)由钢筋单元累计单元主应力
其中,钢筋单元r结点的结点主应变记为向量
(6-4-5)求解第i+1步钢筋单元的弹性模量和刚度矩阵:
取钢筋单元主累计应变
若εs≤fy/es,则钢筋单元未屈服,第i+1步钢筋单元的弹性模量为es,钢筋单元的刚度矩阵为:
其中,
若εs>fy/es,则钢筋单元屈服,第i+1步钢筋单元的弹性模量为es2,钢筋单元的刚度矩阵为:
其中,
(6-5)求解第i+1步钢筋混凝土单元刚度矩阵[ki+1]:
(6-6)迭代步消除不平衡力:
在整体坐标下,记第j迭代步的位移误差向量为
对于第i荷载步第1迭代步的公式如下:
可以求解{δ′i},对于第j迭代步公式如下:
对于第m迭代步的迭代公式如下:
(6-7)由上述步骤可以得到同时考虑mazars损伤模型和钢筋作用的{p}荷载作用下的第i荷载步的结点位移,记为{wi}:
求得第i荷载步的结点位移{wi}后,回到(6-1)求解第i+1荷载步的结点位移{wi+1},且i∈[1,n];
最终将n个荷载步的所有结点位移累加,就可以的得到同时考虑mazars损伤模型和钢筋作用的{p}荷载作用下的结点位移{w*}={w1}+{w2}+…+{wn}。
步骤8所述的对考虑mazars损伤模型和钢筋作用的自由曲面的控制点进行高度调整,具体步骤为:
记l为自由曲面结构的第l个优化步,且l属于正整数集,δz(l)表示第l个优化步和第l+1个优化步的控制点高度差,其求解方法如下所示:
由上式得到控制点高度差δz(l)后,可以求得第l+1优化步的控制点高度
其中,m和n分别为u和v方向的控制点个数;z(l)表示第l优化步控制点高度;
(8-1)应变能差分算法求解梯度:
结构应变能c是z(l)的函数,假定应变能c是关于z(l)的连续可微函数,将c(z(l))进行taylor展开:
可以得到梯度
式中:δzi自由曲面控制点上纵坐标增量;o(δzi)是高阶无穷小量;且ξ∈(0,δzi),c(3)(ξ)表示应变能关于ξ的三阶导数;
(8-2)采用黄金分割法求解步长:
其计算步骤如下,其中f表示黄金分割法的第f迭代步,且f属于正整数集,λf为第f次迭代求得的步长值,当程序循环至满足step2的条件时将λf的值赋给λ(l):
step1:置初始区间[a1,b1]和精度要求g>0,计算试探点λ1和μ1,计算c(λ1)和c(μ1),令f=1;
λ1=a1+0.382(b1-a1)
μ1=a1+0.618(b1-a1)
step2:若bk-ak<g,则停止计算;否则,当c(λk)>c(μk)时,转step3;当c(λf)<c(μf)时,转step4;
step3:置af+1=λf,bf+1=bf,λf+1=μf,μf+1=af+1+0.618(bf+1-af+1);计算应变能c(μf+1),转step5;
step4:置af+1=af,bf+1=μf,μf+1=λf,λf+1=af+1+0.382(bf+1-af+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×1010n·m-2,泊松比μ=0.2。混凝土壳厚度为t=0.1m,壳体受竖向均布荷载作用,两纵边简支约束。结构初始曲面形状如附图2所示,图(a)和(b)分别为该自由结构的正视图和俯视图。
自由曲面创建的控制点位置如附图2(b)空心圆点所示。优化目标为结构应变能最小。自由曲面结构取均布荷载q=120kn/m2,对该自由曲面结构进行创建。考虑和不考虑材料损伤和钢筋作用的曲面创建结果分别如附图3(a)和附图3(b)所示。
表1显示了考虑和不考虑损伤的自由曲面创建结果的结构最大损伤、位移和mises应力。
表1自由曲面结构的最大损伤、位移、mises应力的比较
附图4-图6分别为考虑不考虑材料损伤和钢筋作用的自由曲面创建结果的损伤因子、结点位移和mises应力对比,可以发现采用本方法创建的自由曲面在结构损伤因子、整体刚度和力学响应方面表现更优。
本发明方案所公开的技术手段不仅限于上述技术手段所公开的技术手段,还包括由以上技术特征任意组合所组成的技术方案。
以上述依据本发明的理想实施例为启示,通过上述的说明内容,相关工作人员完全可以在不偏离本项发明技术思想的范围内,进行多样的变更以及修改。本项发明的技术性范围并不局限于说明书上的内容,必须要根据权利要求范围来确定其技术性范围。