结构极限强度分析的三种有限元解法研究

彭大炜, 张世联

(上海交通大学,上海200030)

结构极限强度分析的三种有限元解法研究

彭大炜, 张世联

(上海交通大学,上海200030)

介绍了弧长法、阻尼因子法和准静态法等3种计算结构极限强度的非线性有限元解法及其求解思路。通过Reckling No.23模型数值计算,对上述3种解法各自的计算过程、求解特点和要点进行了比较和归纳。为更准确、高效地运用非线性有限元程序计算分析结构极限强度提出了合理的建议。

弧长法;阻尼因子法;准静态法;极限承载能力

0 前言

近年来,结构设计的思路已经逐步从传统的许用应力设计向极限状态设计发展[1]。研究结构的极限强度水平,确定其与设计载荷水平之间的确切裕度,已成为结构理性设计的核心和基础。

目前,船体结构极限承载能力研究的方法主要有[2]:实船事故调查和模型试验法、直接法(解析法)、以Smith法为代表的逐步崩溃法和非线性有限元法。随着计算机技术的不断发展,非线性有限元法现已成为计算和评估结构极限承载能力最理想的方法,其分析问题的主要流程如图1所示。

在非线性有限元法中,有3种不同的分析结构极限强度的解法。它们分别是弧长法、阻尼因子法和准静态法。本文基于ABAQUS大型通用非线性有限元程序,以Reckling-No.23试验模型为对象,分别采用上述3种解法进行极限强度的数值分析,并研究不同解法的求解过程,对不同解法的应用特点进行了归纳和总结,且还针对如何更准确、高效地运用非线性有限元法计算分析结构极限强度的问题提出合理的建议。

1 理论介绍

从求解本质而言,弧长法和阻尼因子法是以结构非线性静态平衡方程式(1)的求解(牛顿拉-普森迭代)为基础。而准静态方法是以结构非线性运动方程式(2)的显式求解(中心差分法)为基础。

式中:{P}为载荷列阵;{I}为内力列阵;[M]为质量矩阵;{¨u}为加速度列阵。

以下分别简要介绍这3种解法的基本思路。

1.1 弧长法

弧长法的基本思路是通过设置一个参数(弧长l)来控制平衡方程的增量迭代和收敛。将(1)式写成增量格式:

式中:[KT]为切线刚度矩阵;{Δu}为位移增量;{ΔP}为载荷增量;{R}为残差力。

在弧长法中,第i步迭代的载荷增量{ΔP}i,由载荷增量因子Δ λi和参考载荷{Pref}来控制,即

将(4)式代入(3)式即得到弧长法的第i步迭代的增量形格式:

图2 弧长法示意图

弧长法求解时,是以前一步增量计算得到的平衡点为圆心,以弧长增量Δli为半径,通过牛顿拉-普森迭代搜索下一个平衡点,如图2所示。每一步的弧长增量Δli、载荷增量因子Δ λi和位移增量{Δu}i通过下面的约束方程来控制:

通过迭代,直至残差力在容差{R}i范围内。当第i步迭代完成时,有

由于在弧长增量Δli中同时包含了载荷增量{ΔPi}和位移增量{Δui}的信息,运用弧长法能够全程跟踪结构在“加载”(Δ λi＞0)和“卸载”(Δ λi＜0)过程中的载荷-位移平衡路径。但在实际运用时,存在跟踪失败[3]和收敛性问题。目前弧长法还在不断地研究改进。

1.2 阻尼因子法

阻尼因子法仍采用传统的载荷步长增量来进行非线性迭代求解。其思路是通过在平衡方程式(1)中引入阻尼力项来求解结构不稳定崩溃的问题,求解的平衡方程为

其中{Fv}为阻尼力列阵,取决于求解时结构变形的速度,由广义节点速度{v}来表征:

式中:c为阻尼系数[4];M＊为人工质量矩阵;广义节点速度{v}={Δu}/Δt。

加载的初始阶段,结构尚处于稳定状态,此时广义节点速度{v}很小,故阻尼力项{Fv}对平衡方程式(10)几乎没有影响,可以忽略。随着载荷的不断增加,结构趋向于不稳定,当外载{P}已经不能完全由结构内力{I}来平衡时,结构达到极限状态,相当部分的应变能将释放转化为动能,广义节点速率{v}迅速增大。此时,阻尼力项{Fv}通过做功消耗释放的应变能,在平衡方程式(10)中起到维持求解系统的“平衡”作用。

由于采用传统的载荷步长来进行非线性迭代,阻尼因子法无法继续有效跟踪结构在“卸载”过程中的载荷-位移路径,取而代之的是一条几乎水平(斜率为0)的直线。通过考察阻尼力项{Fv}为维持系统“平衡”所做的功所占结构应变能百分比的历史变化曲线,就能够确定结构的极限承载能力。相比弧长法、阻尼因子法的数值收敛性要好一些,能够解决的问题也更广泛。

1.3 准静态法

与弧长法和阻尼因子法的静态求解方法不同,准静态法从本质上讲是一个结构动态求解的过程。在求解时,对式(2)采用中心差分法进行显示的时间积分,由一个增量步的动力学条件计算下一个增量步的动力学条件,直至求解时间结束。准静态法的基本思路就是用慢速加载的动态分析来模拟静态问题,所以,求解关键是要设置合适加载速率。加载速率过快会导致求解结果的局部性(剧烈的结构局部变形),使计算结果偏离“准静态”的要求;而过慢的加载速率意味着较长的加载时间,从而使计算时间大幅度增加。故分析时,通常会取从快到慢多个加载速率进行分析比较,以选定一个合适的加载速率。判断加载速率是否合适的一个重要标准就是分析过程中结构模型的动能与其应变能之比,一般准静态的要求是小于5%。

由于中心差分法是条件稳定的算法[5],在分析时,时间步长Δt必须小于稳定性限制Δtstable才能保证求解的稳定性:

式中:ωmax为结构最高阶固有频率;Δte为结构模型中最小尺寸的单元的稳定时间步长。Δte与单元的特征尺度Le、弹性模量E和材料密度ρ有关,是稳定性限制Δtstable的上界[6]:

在计算中,由于结构最高阶频率ωmax不易求得,故时间步长Δt就取为Δte。

相比弧长法和阻尼因子法,准静态方法最大的优势在于采用中心差分法进行显式时间积分不存在收敛性的问题。因此,准静态法能够很好地求解更复杂结构崩溃问题,如结构的自接触和材料的失效问题。在求解极限状态问题时,时间步长Δt往往较小,准静态“缓慢”加载可能导致求解的机时很长,这时可采用质量放大等方法进行调整。

2 数值计算分析

2.1 Reckling No.23模型试验

Reckling(1979)[7]采用箱型剖面模拟实船剖面进行了系列总纵极限承载能力试验,以研究船体梁在极限状态下各种崩溃模式和剖面应力分布。本文选取加强筋较多的Reckling No.23号模型作为分析对象,模型跨长l=500 mm。表1列出了模型截面尺寸和材料属性。

图3 Reckling有限元模型一阶屈曲模态

表1 Reckling No.23模型截面尺寸和材料属性

在进行极限强度分析时,可将如图3所示模型的第一阶线性屈曲模态的变形形式作为结构的初始缺陷引入模型网格。

2.2 弧长法解

采用弧长法时,在模型两端设置大小相同、方向相反的参考载荷(中垂弯矩),对应各弧长增量步的载荷与位移的表达式如式(7)～式(9)所示。得到的端面弯矩-转角曲线如图4所示:考察曲线峰值点对应的弯矩,即为采用弧长法得到的极限载荷。

2.3 阻尼因子法解

图4 载荷-位移平衡路径曲线(弧长法)

采用阻尼因子法时,在模型端面设置的载荷应超过结构的极限载荷,一般可取为端面的塑性弯矩。在结构达到极限状态前,采用阻尼因子法和采用弧长法得到的端面弯矩-转角曲线几乎重合;当结构达到极限状态时,采用阻尼因子法得到的端面弯矩-转角曲线变成一条水平线,如图5(a)所示。图5(b)、图5(c)表达了中垂加载过程中,阻尼力项所做的功和结构应变能随端部施加的中垂弯矩和端部转角的变化曲线。通过综合考察图5所示的端面弯矩-转角曲线、弯矩-能量曲线和转角-能量曲线来确定结构的极限承载能力。

图5 阻尼因子法

2.4 准静态法解

采用准静态法时,在模型两端施加随时间光滑变化的强迫位移(转角)。当总的强迫位移值一定时,加载速率取决于加载的总时间。当加载时间为1 s时的端部强迫位移加载曲线如图6(a)所示。这里取加载的总时间为1/100 s、1/20 s和1/1 s由快到慢三个加载速率进行计算,输出端面的弯矩反力,得到如图6(b)所示的弯矩-转角曲线。同时考察如图6(c)所示的结构动能和应变能历史曲线,确定合适的加载速率。

图6 准静态法

2.5 讨论

表2列出了三种数值解法的计算结果及其相对试验值的误差。计算表明,虽然三种有限元解法具体的计算思路和过程有所不同,但最终得出的结果基本一致,且与试验值相吻合。下面对不同解法的应用特点与分析的过程进行讨论:

表2 数值计算结果与试验值比较

(1)采用弧长法应对结构模型的规模加以控制,以避免出收敛性问题。

(2)采用阻尼因子法应同时考察载荷-位移曲线和阻尼力项所做的功占应变能百分比的变化曲线。当结构到达极限状态时,阻尼力项所做的功将迅速增加,甚至超过结构的应变能。这表明阻尼力项在平衡方程中占据了主导作用,此时对应的载荷即为结构的极限载荷。

(3)采用准静态方法应取由快到慢若干个加载速率进行分析。加载曲线要求其一阶和二阶导数都是光滑的,从而避免由于加载不连续引起的求解波动。通过考察加载过程中结构动能与应变能之比来判断加载速率是否合适。计算表明,采用合适的加载速率得到的准静态结果与静态计算的结果接近。

3 总结

本文对求解结构极限承载能力的非线性有限元法中的3种解法:弧长法、阻尼因子法和准静态法进行了理论介绍,并通过对Reckling No.23模型试验的数值计算,对3种解法的求解思路、计算特点和关键点进行了归纳。计算分析表明,3种解法都是有效的数值方法,可以相互作为辅助和补充。当求解模型规模较大时,可以考虑采用阻尼因子法和准静态方法进行求解;当求解涉及复杂的结构接触、材料失效等不容易收敛因素时,采用准静态方法可能更为有效。表3总结了3种非线性有限元解法各自的特点。

表3 3种非线性有限元解法的特点

[1] Jeom Kee Paik,Anil Kumar Thayamballi.Ultimate Limit State Design of Steel-Plated Structures[M].Chichester,U K:Wiley,2003.

[2] ISSC.Report of ISSC Special Task Committee VI.2-Ultimate Hull Girder Strength[R].Proceedings of ISSC 2000[C].Nagasaki,Japan,2000.

[3] 李元齐,沈祖炎.弧长控制类方法使用中若干问题的探讨与改进[J].计算力学学报,1998,15(4):414-422.

[4] ABAQUS Analysis User’s Manual.7.1.1 SOLVING NONLINEAR PROBL EM[S].Using the damping factor.

[5] 王勖成.有限单元法[M].北京:清华大学出版社,2003.

[6] 庄茁,张帆,岑松.ABAQUS非线性有限元分析与实例(ABAQUS数码工程师系列丛书)[M].北京:科学出版社,2005.

[7] Reckling K A.Behavior of box girder under bending and shear[R].Proceeding of the International Ship&Offshore Structures Congress(ISSC),Pairs,1979,2:46-49.

Study on Three FEA Methods of Structure Ultimate Strength Analysis

PENG Da-wei, ZHANG Shi-lian

(Shanghai Jiao Tong University,Shanghai 200030,China)

Three algorithms,Arc-Length method,Damping Factor method and Quasi-Static method in nonlinear FEA method of ultimate strength analysis are presented in this paper.The three algorithms have different analysis processes,features and key points,which are discussed and compared with one another through a numeral calculation by Reckling No23 test model.The result of the comparisons gives a reasonable suggestion for better use of nonlinear FEA method of ultimate strength analysis.

Arc-Length method,Damping Factor method,Quasi-Static method,ultimate strength

U661.41

A

图1 非线性有限元分析流程

2009-11-03;修改稿收到日期:2009-12-25

彭大炜(1984-),男,硕士研究生,主要从事船体结构极限强度分析研究。

1001-4500(2010)02-0001-05