一种曲线加筋板有限元分析方法与流程

本发明涉及一种曲线加筋板有限元分析方法，属于加筋板结构建模分析领域。

背景技术：

加筋板是一种由基体结构和加强筋通过一定的铺设方式连接组合而成的结构组件。因其具有加工容易、比强度大、承载力强等优点，且可在同等重量的前提下提高结构屈曲破坏临界压力，加筋板成为工业耐压结构的主要结构形式之一。传统加筋板通常采用直线加筋形式，通过横向、纵向或按照某一角度铺设。直线形式无法提高材料利用率，不便于优化设计。

商业有限元软件在对加筋板结构进行有限元分析时，必须保证筋板节点重合建模，当筋条在板内铺设方案改变时，需重新对板单元划分网格，因此极大降低了建模效率。在对加筋板单元建模中，常用理论为经典薄板理论与一阶剪切变形理论，均无法构造一种面对加筋厚薄板通用的建模分析方法。

技术实现要素：

发明目的：针对传统有限元软件要求筋板节点重合导致建模效率低、且现有建模方法无法实现厚薄板通用的缺陷，本发明提供一种曲线加筋板有限元分析方法，该方法可以实现对曲线加筋厚薄板通用分析，并能实现筋条在板内的任意布置，可通过调整筋条曲率及铺设位置来提高结构整体刚度。

技术方案：本发明所述的一种曲线加筋板有限元分析方法，包括如下步骤：

步骤1，根据DST-BK单元理论得到板单元刚度矩阵Kp及板单元质量矩阵Mp；

步骤2，在筋板交界面处根据筋板位移耦合条件利用板节点位移表达筋节点位移，得到筋板位移协调后的筋单元刚度矩阵Ks和筋单元质量矩阵Ms；

步骤3，将筋单元刚度矩阵Ks、筋单元质量矩阵Ms与板单元刚度矩阵Kp、板单元质量矩阵Mp相加得到加筋板整体结构刚度矩阵和质量矩阵，进行特征值分析，得到曲线加筋板动态特性。

上述步骤2中，利用板节点位移表达筋节点位移包括下述步骤：

步骤21，采用Distmes划分板单元网格，获得板单元整体坐标及自然坐标；采用3阶B样条曲线划分筋单元网格，获得局部坐标系下筋节点的自然坐标并转换为整体坐标；

步骤22，由筋节点的整体坐标和板节点的整体坐标，找出每个筋节点所对应的板单元，并由所对应的板单元的形函数计算出该筋节点对应在板单元内的自然坐标；

步骤23，将筋节点对应在板单元内的自然坐标，结合板单元的形函数，插值得到该处板单元位移，根据筋板位移耦合条件，利用筋单元形函数由板单元节点位移插值得到落在板单元的筋节点位移，并得到筋板位移耦合形函数。

用板节点位移表达的筋节点位移如下式(4)：

式(4)中，us、υs和ws分别为曲线筋对应于x,y,z三个方向的线位移，βsx和βsy为曲线筋的转角；Ns,i为筋单元形函数，Ns,1＝1-ξs-ηs,Ns,2＝ξs,Ns,3＝ηs，ξs和ηs为筋节点的自然坐标；Np,j为板单元形函数，Np,1＝1-ξp-ηp,Np,2＝ξp,Np,3＝ηp，ξp和ηp为板节点的自然坐标；up,j、υp,j和wp,j分别为板单元角节点对应于x,y,z三个方向的线位移，βpx,j和βpy,j为板单元角节点的转角；αk为边中点的转角，Pk表示一组高阶函数，P4＝4ξp(1-ξp-ηp)+4ηp/3-1/2，P5＝4ξpηp+4(1-ξp-ηp)/3-1/2，P6＝4(1-ξp-ηp)ηp+4ξp/3-1/2，Ck和Sk分别为三角形边与x轴所成角度的余弦和正弦值。

将式(4)简写为式(5)，得到筋板位移耦合形函数Nsp，筋板之间的位移关系和坐标关系可表示为：

dsg＝Nspdp,rs＝Nsprp (5)，

式(5)中，dsg表示整体坐标系下筋单元的位移，dp表示整体坐标系下板单元的位移，rs表示整体坐标系下筋单元的坐标，rp表示整体坐标系下板单元的坐标。

由此，筋板位移协调后的筋单元刚度矩阵Ks为：

式(6)中，Ts为筋单元坐标转换矩阵，Bs为筋单元应变矩阵，Ds为筋单元弹性矩阵，Js为筋单元的雅克比矩阵。

筋板位移协调后的筋单元质量矩阵Ms为：

式(7)中，Ts为筋单元坐标转换矩阵，ms为局部坐标系中筋条的质量矩阵，Js为筋单元的雅克比矩阵。

上述步骤3中，根据下式(8)进行特征值分析，得到曲线加筋板的动态特性：

[(Kp+Ks)-ω²(Mp+Ms)]d＝0 (8)

式(8)中，ω为固有频率，d为屈曲模态振型。

有益效果：与现有技术相比，本发明的优点在于：(1)本发明的曲线加筋板有限元分析方法采用板单元节点坐标与位移插值得到筋单元刚度矩阵，筋板单元节点无需重合，大大提高了曲线加筋板的建模效率；(2)利用厚薄板通用的分析方法对曲线加筋板进行建模分析，避免了一阶剪切变形理论在分析曲线加筋板结构时会遇到的剪切锁死问题。

附图说明

图1为本发明的一种曲线加筋板有限元分析方法的流程框图；

图2为实施例中四筋对称曲线加筋板的几何模型；

图3为实施例中曲线加筋板的有限元模型；

图4为实施例中获得的曲线加筋板第1阶振型云图；

图5为实施例中获得的曲线加筋板第2阶振型云图；

图6为实施例中获得的曲线加筋板第3阶振型云图。

具体实施方式

下面结合附图对本发明的技术方案作进一步说明。

本发明的一种曲线加筋板有限元分析方法，如图1，该方法主要包括下述步骤：

步骤1，根据DST-BK单元理论得到板单元刚度矩阵及板单元质量矩阵；

步骤2，在筋板交界面处根据筋板位移耦合条件利用板节点位移表达筋节点位移，得到筋板位移协调后的筋单元刚度矩阵Ks和筋单元质量矩阵Ms；

步骤3，将筋单元刚度矩阵Ks、筋单元质量矩阵Ms与板单元刚度矩阵Kp、板单元质量矩阵Mp相加得到加筋板整体结构刚度矩阵和质量矩阵，进行特征值分析，得到曲线加筋板动态特性。进一步进行收敛性分析，可得到准确模型下的曲线加筋板动态特性收敛结果。

该方法可以实现对曲线加筋厚薄板通用分析，并能实现筋条在板内的任意布置，可通过调整筋条曲率及铺设位置来提高结构整体刚度。

以一对称加筋的曲线加筋板为例，其几何模型如图2，利用本发明的曲线加筋板有限元分析方法分析该结构的动态特性，具体包括以下步骤：

1、根据DST-BK单元理论得到板单元刚度矩阵Kp及板单元质量矩阵Mp；

2、采用Distmesh进行板单元网格划分，采用3阶B样条曲线划分筋单元网格，得到曲线加筋板的有限元模型，如图3，同时得到筋节点自然坐标，计算筋节点对应在板单元内的自然坐标，计算筋板位移耦合形函数，用板节点位移表示筋位移场；

具体而言，先采用网格生成工具Distmesh生成3节点三角形板单元网格，获得板单元整体坐标及自然坐标，并采用3阶B样条曲线生成3节点筋单元网格，获得局部坐标系下筋节点的自然坐标并转换为整体坐标；然后由筋节点的整体坐标和板节点的整体坐标，找出每个筋节点所对应的板单元，并由所对应的板单元内3点的形函数计算出该筋节点对应在板单元内的自然坐标；由上一步所得筋节点所对应在板单元内的自然坐标，并结合3节点板单元的形函数，插值得到该处板单元位移，根据筋板交界处位移协调条件，认为该处筋板位移相同，因此采用三点处的筋板结合处位移，利用筋单元形函数插值得到筋条内任意点的位移场。

对于筋条来讲，采用Timoshenko梁单元理论，筋条内的任意一点位移场和坐标可由单元内的3点位移和坐标插值得到：

式(1)中，usg,vsg,wsg分别为筋单元任一点在整体坐标系下沿着x,y,z三个方向的线位移，βxsg,βysg为筋单元任一点在整体坐标系下的转角位移；usg,i,vsg,i,wsg,i为筋单元三个节点在整体坐标系下沿着x,y,z三个方向的线位移，βxsg,i和βysg,i为筋单元三个节点在整体坐标系下的转角位移；Ns,i表示筋单元形函数，其中，Ns,1＝1-ξs-ηs,Ns,2＝ξs,Ns,3＝ηs，ξs和ηs为筋节点的自然坐标。

在3节点9自由度DST-BK单元位移表示中引入线性插值的面内位移up,υp，单元内任意一点位移可由节点位移插值得到：

式(2)中，up,vp,wp分别为整体坐标系下板单元任一点沿着x,y,z方向的线位移，βpx,βpy为整体坐标系下板单元任一点的转角位移；Np,j表示板单元形函数，其中Np,1＝1-ξp-ηp,Np,2＝ξp,Np,3＝ηp，ξp和ηp为板节点的自然坐标；

up,j，υp,j和wp,j分别为板单元角节点对应于x,y,z三个方向的线位移；βpx,j和βpy,j为板单元角节点的转角，αk是边中点的转角，Pk表示一组高阶函数，P4＝4ξp(1-ξp-ηp)+4ηp/3-1/2，P5＝4ξpηp+4(1-ξp-ηp)/3-1/2，P6＝4(1-ξp-ηp)ηp+4ξp/3-1/2，Ck和Sk分别为三角形边与x轴所成角度的余弦和正弦值。

在筋板交界处，根据筋板位移协调条件，认为筋板位移相同，则落在板单元的筋节点位移可以由对应的板单元节点位移插值得到：

式(3)中，us,i,υs,i,ws,i表示落在某一板单元内的筋节点线位移，βsx,i,βsy,i表示落在某一板单元内的筋节点转角。

将(2)式、(3)式带入(1)式，可以得到：

式(4)中，us、υs和ws分别为曲线筋对应于x,y,z三个方向的线位移，βsx和βsy为曲线筋的转角。

将式(4)简写为式(5)，得到筋板位移耦合形函数Nsp，筋板之间的位移关系和坐标关系可表示为：

dsg＝Nspdp,rs＝Nsprp (5)，

式(5)中，dsg表示整体坐标系下筋单元的位移，dp表示整体坐标系下板单元的位移，rs表示整体坐标系下筋单元的坐标，rp表示整体坐标系下板单元的坐标。

3、采用上述筋板位移耦合形函数Nsp，结合筋单元坐标转换矩阵Ts、筋单元应变矩阵Bs和筋单元弹性矩阵Ds，通过高斯积分得到筋板位移协调后的筋单元刚度矩阵Ks：

式(6)中，Js为筋单元的雅克比矩阵。

4、采用上述筋板位移耦合形函数Nsp，结合筋单元坐标转换矩阵Ts，通过高斯积分得到筋单元质量矩阵Ms：

式(7)中，ms为局部坐标系中筋条的质量矩阵。

5、由上述得到的筋单元刚度矩阵Ks、筋单元质量矩阵Ms与板单元刚度矩阵Kp、板单元质量矩阵Mp相加得到曲线加筋板整体结构刚度矩阵和质量矩阵。

6、根据下式进行特征值分析，得到曲线加筋板动态特性。

[(Kp+Ks)-ω²(Mp+Ms)]d＝0 (8)，

式(8)中，ω为固有频率，d为屈曲模态振型。

7、进行收敛性分析，得到准确模型下的曲线加筋板动态特性收敛结果，图4～6为实施例中四筋对称曲线加筋板前三阶振型。