CN112966421B - 使用p型有限元法计算薄板结构屈曲载荷因子和相应屈曲形状的计算方法 - Google Patents
使用p型有限元法计算薄板结构屈曲载荷因子和相应屈曲形状的计算方法 Download PDFInfo
- Publication number
- CN112966421B CN112966421B CN202110283019.1A CN202110283019A CN112966421B CN 112966421 B CN112966421 B CN 112966421B CN 202110283019 A CN202110283019 A CN 202110283019A CN 112966421 B CN112966421 B CN 112966421B
- Authority
- CN
- China
- Prior art keywords
- buckling
- finite element
- calculation
- calculating
- load factor
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 46
- 238000004364 calculation method Methods 0.000 title claims abstract description 27
- 239000011159 matrix material Substances 0.000 claims abstract description 29
- 239000000463 material Substances 0.000 claims description 12
- 239000013598 vector Substances 0.000 claims description 4
- 238000002474 experimental method Methods 0.000 claims description 2
- 238000006073 displacement reaction Methods 0.000 description 8
- 229910000831 Steel Inorganic materials 0.000 description 4
- 238000001228 spectrum Methods 0.000 description 4
- 239000010959 steel Substances 0.000 description 4
- 238000005452 bending Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000006243 chemical reaction Methods 0.000 description 1
- 230000006835 compression Effects 0.000 description 1
- 238000007906 compression Methods 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 238000010008 shearing Methods 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
- G06T17/20—Finite element generation, e.g. wire-frame surface description, tesselation
-
- G—PHYSICS
- G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
- G16C—COMPUTATIONAL CHEMISTRY; CHEMOINFORMATICS; COMPUTATIONAL MATERIALS SCIENCE
- G16C10/00—Computational theoretical chemistry, i.e. ICT specially adapted for theoretical aspects of quantum chemistry, molecular mechanics, molecular dynamics or the like
-
- G—PHYSICS
- G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
- G16C—COMPUTATIONAL CHEMISTRY; CHEMOINFORMATICS; COMPUTATIONAL MATERIALS SCIENCE
- G16C60/00—Computational materials science, i.e. ICT specially adapted for investigating the physical or chemical properties of materials or phenomena associated with their design, synthesis, processing, characterisation or utilisation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2111/00—Details relating to CAD techniques
- G06F2111/04—Constraint-based CAD
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2113/00—Details relating to the application field
- G06F2113/24—Sheet material
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
Landscapes
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Physics & Mathematics (AREA)
- Computing Systems (AREA)
- Geometry (AREA)
- Life Sciences & Earth Sciences (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Bioinformatics & Computational Biology (AREA)
- General Physics & Mathematics (AREA)
- Evolutionary Computation (AREA)
- General Engineering & Computer Science (AREA)
- Computer Hardware Design (AREA)
- Computer Graphics (AREA)
- Software Systems (AREA)
- Health & Medical Sciences (AREA)
- General Health & Medical Sciences (AREA)
- Spectroscopy & Molecular Physics (AREA)
- Investigating Strength Of Materials By Application Of Mechanical Stress (AREA)
Abstract
本发明涉及一种使用p型有限元法计算薄板结构屈曲载荷因子和相应屈曲形状的计算方法,属于结构稳定性校核技术领域。本发明将p型有限元法应用于结构稳定性分析领域,求解在相应几何和约束条件下的应力场,得到相应的几何刚度矩阵并应用Lanczos迭代方法计算屈曲特征值,进而得到屈曲载荷因子和屈曲形状。该方法能减少计算成本,提高计算的收敛率和计算精度。
Description
技术领域
本发明涉及一种使用p型有限元法计算薄板结构屈曲载荷因子和相应屈曲形状的计算方法,属于结构稳定性校核技术领域。
背景技术
屈曲是结构常见的破坏形式之一,尤其是对于一些长、薄结构,具备这样性质的板壳结构广泛应用于航空、航天、航海、军工等领域中。飞机机翼、尾翼上的翼面壁板、梁腹板和机身上的蒙皮、隔框等常采用薄壁结构,这些结构都具有一个共同的特征,即一个方向的尺寸大大小于另外两个方向的尺寸,可以简化为板壳结构,当这种结构受到压缩、弯曲和剪切等外荷载单独或联合作用时,屈曲失稳是最常见的失效模式之一。薄板的屈曲分析在工程设计中至关重要,对于大多数薄板结构,屈曲失稳先于强度破坏。
对于屈曲分析,主要应用伽辽金法、摄动法和有限元法。伽辽金法在求解过程中不涉及变分问题,简单易行,应用较广,但是在微分方程算子不是线性算子时很难求解;摄动法在求解过程中需要解渐近展开式,由于渐近级数一般是发散的,在求解时可能在级数展开项的前几项表现为逐渐收敛形式,而在级数项增加到一定程度后又表现为发散形式,要解决这一问题就需要较高的数学专业技能,增加了求解难度;有限元法采用将复杂结构离散化为一定数量的单元求数值解,为了提高数值解的计算精度,必须划分更小的单元。传统有限元法划分单元较多,计算效率偏低。
本发明专利运用p型有限元法进行薄板结构的特征值屈曲分析,利用p型有限元法单元数量少、前处理少及收敛速度快的优点,提出了一种基于p型有限元法计算薄板结构屈曲载荷系数和相应屈曲形状的方法。
发明内容
本发明将p型有限元法应用于结构稳定性分析领域,求解在相应几何和约束条件下的应力场,得到相应的几何刚度矩阵并应用Lanczos迭代方法计算屈曲特征值,进而得到屈曲载荷因子和屈曲形状。该方法能减少计算成本,提高计算的收敛率和计算精度。本发明通过以下技术方案实现。
使用p型有限元法计算薄板结构屈曲载荷因子和相应屈曲形状的计算方法,其步骤如下:
步骤1、确定薄板结构的几何尺寸和材料参数,建立三维有限元模型,该步骤包括:
步骤1.1、根据结构的实际受力情况,确定结构的受力特征;
步骤1.2、根据结构所使用的材料确定材料属性,如有必要应通过实验测定材料的弹性模量、泊松比。
步骤1.3、建立三维有限元模型。
步骤1.4、施加载荷和边界条件:根据确定的载荷和边界条件,在相应的边界上施加相对应的载荷和边界条件;
步骤2使用p型有限元法计算在对应荷载和约束条件下的线性解。具体步骤如下:
步骤2.1、结合所指定的插值多项式阶数,依据方程
Ka=F (1)
求解结构位移列阵a,其中,K为整体刚度矩阵,K=∑GTKeG,Ke为单元刚度矩阵,a为结构位移列阵,F为结构节点载荷列阵,Ke,G,表示结构节点自由度与单元节点自由度的转换矩阵。
步骤2.2、依据单元刚度矩阵公式
单元等效节点载荷列阵
其中:单元内部节点力
外部节点力
上式(2)至上式(5)中,Ω表示为单元内部;B=LNI,B表示应变矩阵,L表示微分算子,NI表示高阶形函数矩阵,I表示单元的标号,Γt表示单元的外部边界;D表示应力矩阵,b表示体力,表示荷载边界条件,这些量由步骤1中确定的材料属性、载荷及边界条件计算得到;把上式(2)至上式(5)代入到上式(1)中,再对线性方程组(1)进行求解,则得到了结构位移列阵a;
步骤2.3、依据公式aI=Ga进一步得到了单元位移列阵aI;
步骤2.4、依据公式u=NIaI,ε=Lu,σ=Dε=DBaI,求解得到应力场σ、应变场ε、位移场u。
本发明采用p型有限元法求解位移场和应力应变场。p型有限元法可在固定网格下,通过采用阶谱类型的高阶形状函数来提高计算精度,其主要特征在于在三维标准六面体单元上采用了一种阶谱类型的高阶形状函数。本发明采用六面体网格单元对薄板屈曲载荷因子和屈曲形状进行数值模拟分析。
阶谱类型的三维标准六面体单元高阶形状函数的构造形式如下:
三维标准六面体单元如图1所示,其形函数为:
(1)节点模式形状函数:
节点模式形状函数共有8个,其定义如下:
(2)边模式形状函数:
边模式形状函数共有12(p-1)个,其定义如下:
与连接节点1和节点2的边相关联的边模式形状函数为:
其中,
类似地,与连接节点2和节点3的边相关联的边模式形状函数为:
其余边模式形状函数以此类推可得。
(3)面模式形状函数:
面模式形状函数共有3(p-2)(p-3),p≥4个,其定义如下:
与连接节点1,2,5,6的面相关联的面模式形状函数为:
其中,下标m=m(i,j)由编号约定决定。
其余面模式形状函数以此类推可得。
(4)内部模式形状函数:
内部模式形状函数共有p≥6个,其定义如下:
其中,下标m=m(i,j,k),i,j,k=2,3,…,p-4,i+j+k=6,7,…,p。
插值多项式的阶次可依次从p=1逐渐提高,提升插值多项式低阶的刚度矩阵可继续沿用,只需计算高阶部分,避免了刚度矩阵低阶部分的重复计算,具有良好的承袭性,节约了前处理的成本。
步骤3:通过所求得的线性解中的应力场,计算所对应的几何刚度矩阵;再使用几何刚度矩阵进行屈曲载荷因子特征值求解:具体为
步骤3.1、通过计算得到的应力场求得几何刚度矩阵;
步骤3.2、将所求得的几何刚度矩阵用于特征值计算,将广义的特征值问题Kφ=λMφ转化为标准特征值问题Kφ=λφ,即将M化为单位阵;
其中K为刚度矩阵,M质量矩阵,φ为特征矢量,λ即为所求得的屈曲载荷因子特征值;
采用Lanczos迭代方法进行特征值计算;标准特征值问题:Kφ=λφ
任何初始向量U,||U||k=1,U0=0,其中||||k为k范数
采用三项递推公式:{Uk+1}=(K{Uk}-αk{Uk}-βk{Uk-1})/βk+1
其中,β1=0
αk={Uk}TK{Uk}
βk+1=||K-{Uk}-αk{Uk}-βk{Uk-1}||2
这里k=1,2,3,...,m-1≤n,||||2代表2范数,m为最高阶特征值个数,n为方阵阶数,完成迭代计算后,最终得到两组解:用来建立预屈曲应力状态的线性解和特征值屈曲载荷因子解;
步骤4:根据所求得的屈曲载荷因子是否满足精度要求,并判断能量范数误差是否满足一定的要求,如屈曲载荷因子未满足精度要求或者能量范数误差过大,提高插值多项式的阶数并增加网格剖分数量,返回步骤2;屈曲载荷因子满足精度要求,根据屈曲载荷因子获得相应屈曲形状。
所述步骤2中p型有限元法在单元数量确定的情况下,通过提升插值多项式的阶数来提高计算结果的精度。当单元厚度较小时,使用低阶插值函数会使剪切刚度占主导地位,发生剪切锁定,因此插值多项式的阶数应大于等于3。
上述p型有限元法的插值多项式是阶谱类型形状函数,提升插值多项式阶数时,低阶时的刚度矩阵可继续沿用,只需计算新增高阶项,装配形成新的总体刚度矩阵,具有较好的承袭性
上述相应屈曲形状是通过以下方式得到:屈曲载荷因子乘以步骤1中所施加荷载得到屈曲临界荷载,再以屈曲临界荷载作为外荷载以步骤2.1—2.4所述方式计算得到位移场,据此绘出薄板位移图即为相应的屈曲形状。
本发明的有益效果是:该方法能减少计算成本,提高计算的收敛率和计算精度。
附图说明
图1是三维有限元网格单元示意图;
图2是薄钢板结构示意图;
图3是薄钢板屈曲分析网格示意图;
图4是实施例1和实施例2屈曲形状示意图。
具体实施方式
下面结合附图和具体实施方式,对本发明作进一步说明。
实施例1
根据图2所示薄钢板构件,其中薄板长为a=3000mm,宽b=a/m,厚h=10mm,薄板两端受压应力σ=1.0MPa作用,边界条件为四边简支约束。薄钢板的材料参数:杨氏模量E=205800MPa,泊松比ν=0.3,m为变化参数,本实施例实例选择从1到2。
模型的网格划分如图3所示。当薄板的几何参数变化时,网格采用均匀剖分,并且单元数量均为9。
由采用平衡法求解的理论公式可得此时对应的屈曲载荷因子的理论解。
上式中:k为屈曲系数,与薄板的长宽比m=a/b有关,当m为整数时,k=kmin=4,E为弹性模量,ν为泊松比,t为薄板厚度,b为薄板的宽度,下表为m=1,…,5时临界屈曲载荷的理论值。
表1
基于p型有限元求解该薄板构件的线性解,并通过线性解得到几何刚度矩阵,再得到临界屈曲载荷因子,将屈曲载荷因子与所施加的外力相乘最终得到屈曲荷载。下表给出了p=3,4,…,8且m=1,2时有限元解的能量范数误差和对应的约束条件下得到的临界屈曲荷载,及与理论解相比的相对误差。
表2
基于p型有限元法计算临界屈曲载荷,随着插值多项式阶数p的增加,有限元解及有限元解能量范数误差逐渐收敛。基于上表可看出,本发明计算临界屈曲载荷时,具有前处理少,收敛速度快,采用较低的计算成本,便可获得较高精度的数值解。
实施例2
实例2模型与实例1类似,在相同网格和材料属性条件下,改变薄板的长宽比,继续计算m=4,5时的临界屈曲载荷,具体数据如下表:
表3
由上表可看出改变薄板的长宽比,网格条件相同条件下,对比理论解,此时临界屈曲载荷的精度较高,且具有良好的收敛性。
实施例1与实施例2屈曲形状如图4所示。
以上对本发明的具体实施方式作了详细说明,但是本发明并不限于上述实施方式,在本领域普通技术人员所具备的知识范围内,还可以在不脱离本发明宗旨的前提下作出各种变化。
Claims (2)
1.一种使用p型有限元法计算薄板结构屈曲载荷因子和相应屈曲形状的计算方法,其特征在于步骤如下:
步骤1:构建薄板构件三维有限元模型,该步首先创建薄板构件的几何模型,接着对所建立几何模型进行网格剖分;其次对要进行屈曲分析的薄板构件进行材料参数的设定;然后施加对应的荷载与边界条件,其步骤如下:
步骤1、确定薄板结构的几何尺寸和材料参数,建立三维有限元模型,该步骤包括:
步骤1.1、根据结构的实际受力情况,确定结构的受力特征;
步骤1.2、根据结构所使用的材料确定材料属性,通过实验测定材料的弹性模量、泊松比;
步骤1.3、建立三维有限元模型;
步骤2:使用p型有限元法计算在对应荷载和约束条件下的线性解中的应力场;
步骤3:通过所求得的线性解中的应力场,计算所对应的几何刚度矩阵;再使用几何刚度矩阵进行屈曲载荷因子特征值求解:具体为
步骤3.1、通过计算得到的应力场求得几何刚度矩阵;
步骤3.2、将所求得的几何刚度矩阵用于特征值计算,将广义的特征值问题Kφ=λMφ转化为标准特征值问题Kφ=λφ,即将M化为单位阵;
其中K为刚度矩阵,M质量矩阵,φ为特征矢量,λ即为所求得的屈曲载荷因子特征值;
采用Lanczos迭代方法进行特征值计算;标准特征值问题:Kφ=λφ
任何初始向量U,||U||k=1,U0=0,其中|| ||k为k范数
采用三项递推公式:{Uk+1}=(K{Uk}-αk{Uk}-βk{Uk-1})/βk+1
其中,β1=0
αk={Uk}TK{Uk}
βk+1=||K-{Uk}-αk{Uk}-βk{Uk-1}||2
这里k=1,2,3,...,m-1≤n,|| ||2代表2范数,m为最高阶特征值个数,n为方阵阶数,完成迭代计算后,最终得到两组解:用来建立预屈曲应力状态的线性解和特征值屈曲载荷因子解;
步骤4:根据所求得的屈曲载荷因子是否满足精度要求,并判断能量范数误差是否满足一定的要求,如屈曲载荷因子未满足精度要求或者能量范数误差过大,提高插值多项式的阶数并增加网格剖分数量,返回步骤2;屈曲载荷因子满足精度要求,根据屈曲载荷因子获得相应屈曲形状。
2.根据权利要求1所述的使用p型有限元法计算薄板结构屈曲载荷因子和相应屈曲形状的计算方法,其特征在于:所述步骤2中p型有限元法在单元数量确定的情况下,通过提升插值多项式的阶数来提高计算结果的精度,插值多项式的阶数应大于等于3。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202110283019.1A CN112966421B (zh) | 2021-03-16 | 2021-03-16 | 使用p型有限元法计算薄板结构屈曲载荷因子和相应屈曲形状的计算方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202110283019.1A CN112966421B (zh) | 2021-03-16 | 2021-03-16 | 使用p型有限元法计算薄板结构屈曲载荷因子和相应屈曲形状的计算方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN112966421A CN112966421A (zh) | 2021-06-15 |
CN112966421B true CN112966421B (zh) | 2023-12-26 |
Family
ID=76277922
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202110283019.1A Active CN112966421B (zh) | 2021-03-16 | 2021-03-16 | 使用p型有限元法计算薄板结构屈曲载荷因子和相应屈曲形状的计算方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN112966421B (zh) |
Families Citing this family (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN113705040B (zh) * | 2021-08-03 | 2024-03-22 | 大连理工大学 | 结构损伤分析的近场有限元法及在商用软件中的实现方法 |
CN117390778A (zh) * | 2022-12-22 | 2024-01-12 | 北京强度环境研究所 | 一种考虑薄板应变强化效应的塑性安定上限载荷计算方法 |
Citations (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN108108578A (zh) * | 2018-01-30 | 2018-06-01 | 南京理工大学 | 基于无网格法的fg-grc板屈曲载荷因子的数值算法 |
CN109214041A (zh) * | 2018-07-19 | 2019-01-15 | 东南大学 | 一种考虑力载荷的板结构屈曲温度分析方法 |
CN109918712A (zh) * | 2019-01-23 | 2019-06-21 | 昆明理工大学 | 一种基于p型有限元法和围线积分法求解复合型应力强度因子的计算方法 |
CN112014018A (zh) * | 2020-09-01 | 2020-12-01 | 西南交通大学 | 一种基于超声层析成像的应力场测量方法 |
-
2021
- 2021-03-16 CN CN202110283019.1A patent/CN112966421B/zh active Active
Patent Citations (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN108108578A (zh) * | 2018-01-30 | 2018-06-01 | 南京理工大学 | 基于无网格法的fg-grc板屈曲载荷因子的数值算法 |
CN109214041A (zh) * | 2018-07-19 | 2019-01-15 | 东南大学 | 一种考虑力载荷的板结构屈曲温度分析方法 |
CN109918712A (zh) * | 2019-01-23 | 2019-06-21 | 昆明理工大学 | 一种基于p型有限元法和围线积分法求解复合型应力强度因子的计算方法 |
CN112014018A (zh) * | 2020-09-01 | 2020-12-01 | 西南交通大学 | 一种基于超声层析成像的应力场测量方法 |
Also Published As
Publication number | Publication date |
---|---|
CN112966421A (zh) | 2021-06-15 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN112966421B (zh) | 使用p型有限元法计算薄板结构屈曲载荷因子和相应屈曲形状的计算方法 | |
Li et al. | Nonlinear dynamic response of sandwich beams with functionally graded negative Poisson’s ratio honeycomb core | |
Garcia | Numerical investigation of nonlinear aeroelastic effects on flexible high-aspect-ratio wings | |
CN105183996A (zh) | 面元修正与网格预先自适应计算方法 | |
CN111597632B (zh) | 一种基于刚性多连杆机构驱动的变形翼结构设计方法 | |
CN113886967B (zh) | 多巡航工况的大型飞机机翼气动弹性优化方法 | |
CN102514709B (zh) | 一种采用格栅结构的飞行器机翼盒段及设计方法 | |
CN113191040A (zh) | 一种考虑结构稳定性的单材料结构拓扑优化方法和系统 | |
CN113723027A (zh) | 一种针对弹性飞机的静气动弹性计算方法 | |
CN102722606A (zh) | 一种降低直升机旋翼桨毂振动载荷的方法 | |
Nampally et al. | Nonlinear finite element analysis of lattice core sandwich beams | |
CN111159636A (zh) | 基于绝对节点坐标描述的柔性多体系统动力学半解析灵敏度分析方法 | |
CN110704953B (zh) | 一种大展弦比机翼静气弹性能设计敏度的分析方法 | |
Palacios et al. | Static nonlinear aeroelasticity of flexible slender wings in compressible flow | |
Kaplunov et al. | A low-frequency model for dynamic motion in pre-stressed incompressible elastic structures | |
CN103942381B (zh) | 用于飞机铝合金结构性能预测的状态近场动力学方法 | |
Zhang et al. | A morphing wing with cellular structure of non-uniform density | |
Gandhi et al. | Design optimization of a controllable camber rotor airfoil | |
Liu et al. | Solving topology optimization problems by the guide-weight method | |
CN117094114A (zh) | 一种理想弹塑性薄板塑性安定上下限载荷的数值计算方法 | |
Cramer et al. | Modeling of tunable elastic ultralight aircraft | |
CN113486512B (zh) | 一种功能梯度变厚度叶片模型的颤振分析方法 | |
Boston et al. | Multistable honeycomb architecture for spanwise wing morphing | |
Tsushima et al. | Structural and aerodynamic models for aeroelastic analysis of corrugated morphing wings | |
Li et al. | Effects of unbalanced lamination parameters on the static aeroelasticity of a high aspect ratio wing |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |