基于隔离非线性的实体单元模型与计算效率分析

李佳龙1,李 钢1,李宏男1,2

(1.大连理工大学海岸和近海工程国家重点实验室,辽宁,大连 116024;2.沈阳建筑大学土木工程学院,辽宁,沈阳 110168)

摘 要:实体有限元模型计算中往往需要较多的计算单元与结点数量,且这些单元状态判定以及大规模的刚度矩阵分解将消耗大量的计算资源,计算效率低。该文基于隔离非线性法理论建立了线性四面体与六面体等参单元分析模型,采用直接积分格式的6积分点替代六面体等参单元的8高斯点作为非线性应变插值点,能够在保证计算精度的同时提高单元状态判定效率。控制方程采用Woodbury公式与组合近似法联合求解,使得整个求解过程只有矩阵回代以及矩阵与向量的乘积,进一步提高了求解效率。基于时间复杂度的计算效率分析表明:随着结点自由度数目的增加,该文方法的计算效率相对传统变刚度法显著提高,数值算例验证了实体单元模型的正确性以及算法的高效性。

关键词:隔离非线性法;实体单元;直接积分;Woodbury公式;组合近似法;效率分析

实体单元作为一种常用的精细化有限元分析模型,理论上能够适用任意几何形状的结构且能全方面反映结构的微观状态,广泛应用于结构分析领域[1]。目前,结构非线性分析主要采用基于位移有限元的变刚度求解方法,非线性迭代求解时需要对单元状态进行判定并重新合成与分解新的刚度矩阵,实体单元模型较多的单元与结点使得上述过程消耗大量的计算资源,计算效率低[2]。在材料非线性分析中,变刚度法要求一旦有单元进入非线性状态即对整个刚度矩阵更新,通常结构仅有局部单元进入非线性状态而大部分单元仍保持弹性。对于局部非线性问题,Akgün等[3]提出快速修改结构切线刚度的Woodbury公式直接法,将结构的切线刚度表示成弹性刚度的低秩修正,避免了对整体刚度的重新合成与分解,但随着结构局部修改部分比例增加,其计算效率将大幅降低甚至低于完全分析。Triantafyllou等[4]采用初始刚度矩阵与滞回矩阵进行结构非线性计算,将非线性效应作为附加荷载作用于原始模型,高效且准确。相较于直接法,基于预处理共轭梯度法[5]发展而来的组合近似法[6-7],通过对初始点精确分析,将二项式级数展开的低阶项作为基向量,通过对刚度矩阵进行降阶处理,缩减了问题的计算规模,具有求解精度高、计算速度快等优点。

隔离非线性法[8](inelasticity-separated finite element method,IS FEM)作为一种结构非线性分析新方法,将结构局部非线性隔离开来,采用Woodbury公式直接求解控制方程,实现了局部非线性的问题的高效求解。其主要思想是,将单元应变分解为线弹性与非线性两部分,通过在单元中设置非线性应变插值点,以插值的形式建立单元的非线性应变场,以结点位移与引入的非线性应变同时作为控制方程的未知变量。控制方程的求解难点在于Woodbury公式中Schur补矩阵与向量之间的运算,规模较大时直接分解求解的计算效率低[9],若对Schur补矩阵进行稀疏化处理将产生较大误差,必须对计算结果进行修正方可接受。李钢等[8-10]使用该方法对纤维梁有限元进行了非线性分析,而对于大多数不能简化为梁柱单元的实体结构,有必要建立相应的实体单元分析模型,对实际工程结构的非线性分析具有重大意义。

隔离非线性法在分析实体单元模型时,Schur补矩阵的规模一般较大且能够稀疏化的程度较低,其常矩阵的预处理计算量大并将占据较大的存储空间。隔离非线性法以非线性应变插值点为基准点进行单元状态判定,实体单元模型较多的插值点使得单元判定效率较低。减缩积分单元[11]相对全积分单元虽可提高单元状态判定效率但对边界条件要求较高,且可能导致的零能模式使得计算结果失真。Ueda等[12]通过计算单元结点力来判断结点的弹塑性状态,避免了单元积分点状态的判定。牛辉等[13]通过修改弹性系数矩阵、约束相应的相对位移并引入各类构件的简化假定提出了实体退化单元模型,降低了计算规模。鉴于此,本文基于隔离非线性理论建立了线性的四面体以及六面体等参单元分析模型。以单元应变分布模式为基础,给出了实体单元非线性应变插值点选择规则,以及相应的非线性应变插值函数。将直接积分格式[14]的6积分点替代线性六面体等参单元的8高斯点作为单元的非线性应变插值点,可在保证数值精度的情况下提高单元状态判定效率。控制方程采用Woodbury公式与组合近似法联合求解,避免了Schur补矩阵的分解求解以及常矩阵的预处理计算,整个求解过程只有矩阵回代以及矩阵与向量的乘积,在提高求解效率的同时也具有足够的计算精度。基于时间复杂度理论[15]定量分析了Woodbury公式与组合近似法联合求解控制方程的计算效率,表明:随着结点自由度数目增加,本方法相对于传统变刚度法具有更高的计算效率。编制了隔离非线性法实体单元模型的计算程序,数值算例表明本文所建立的实体单元模型正确,实际计算时间相对于传统变刚度法更少,显示本文求解方法正确且高效。

1 隔离非线性实体单元模型

1.1 应变分解

隔离非线性法运用应变分解思想,将单元应变分解为线弹性与非线性两部分,图1给出了增量形式的应力-应变关系分解方法。当材料进入非线性状态后,应变增量Δε总是可以分解为线弹性应变增量Δε′与非线性应变增量Δε′′,应力增量Δσ可通过弹性矩阵D表示:

应力增量Δσ又可通过弹塑性矩阵Dt表示成线性化增量的形式,即有:

图1 应力-应变分解
Fig.1 Stress-strain decomposition

1.2 实体单元模型

1.2.1 线性六面体等参单元

图2为典型的线性六面体等参单元,相较于传统有限元单元模型,隔离非线性法额外要求在单元中设置若干个非线性应变插值点,将结点位移增量Δa和引入插值点的非线性应变增量同时作为未知变量。单元的位移场增量Δu以及应变场增量Δε可通过Δa表示,单元的非线性应变场增量 Δε′′以插值的形式表示,即:

图2 线性六面体等参单元
Fig.2 Linear hexahedron isoparametric element

式中:N为位移插值矩阵;B为应变矩阵;C为非线性应变插值矩阵;为非线性应变增量;Ci分别为插值点i(i=1,2,···,n)的非线性应变插值函数与非线性应变增量。Ci应具备完备性性质,其中δij为Kronecker delta函数[11]I6为6阶单位矩阵。

由式(2)和式(5),可得插值点i的非线性应变:

式中,Δεi、ΔσiBi分别为插值点i处的应变增量、应力增量与应变矩阵。

线性六面体等参单元可将2×2×2个高斯点作为非线性应变插值点,单元内的非线性应变分布采用基于物理坐标的多项式函数,该多项式函数的项次分别取此时非线性应变插值函数为:

线性四面体单元为常应力-常应变单元,单元内任意一点均可作为非线性应变插值点,此时插值矩阵C应为6阶单位矩阵,即:

结合式(5)的单元应变场以及式(6)的单元非线性应变场,可将式(2)和式(3)的单元应力场分别改写为:

1.2.2 直接积分的6点插值模型

隔离非线性法的单元状态判定以非线性应变插值点为基准点,插值点数量将直接影响单元状态判定效率,因而有必要对线性六面体等参单元的非线性应变插值点数量以及布置位置进行研究,使得其能够在保证单元计算精度的同时提高单元状态判定效率。

基于立方体单元的直接积分规则[14],将其6个直接积分点替代线性六面体等参单元的2×2×2个高斯点作为非线性应变插值点,此时单元状态判定数目降低1/4。6积分点位置有两种布置方案,其在物理坐标中的位置分别见表1和表2所示。

方案1:将表1中的6个直接积分点作为非线性应变插值点,以作为非线性应变插值多项式函数的项次,此时非线性应变插值函数为:

表1 6点插值布置方案1
Table 1 Scheme 1 of 6 interpolation point layout

n ξ η ζ 权重Hi 阶次2 a b -c 3 -a -b c 1 a -b -c 4/3 3阶4 -a b c 5 -2a0 -c 6 2a 0 c abc===111,,623

表2 6点插值布置方案2
Table 2 Scheme 2 of 6 interpolation point layout

n ξ η ζ 权重Hi 阶次1 1 0 0 2 -1 0 0 3 0 1 0 4 0 -1 0 5 0 0 1 6 0 0 -1 4/3 3阶

方案2:将表2中的6个直接积分点作为非线性应变插值点,采用多项式函数构造单元的非线性应变插值函数将导致系数矩阵奇异,故采用构造法[11]直接构造插值函数Ci。首先构造点1、2、3的非线性应变插值函数分别为然后将点4、5、6的非线性应变插值函数C4C5C6表示成以η(或ζ)方向一次与ζ、ξ(或ξη)方向二次的拉格朗日多项式的乘积,最后对进行修正。

1.3 控制方程

隔离非线性法作为一种基于位移有限元的分析方法,采用最小位能原理建立控制方程。对于在应力增量Δσ、体积荷载增量Δf以及边界荷载增量ΔT作用下满足平衡条件的变形体结构,在满足几何协调条件下的虚位移增量δu)以及虚应变增量δε)上的总虚功等于零,则有:

式中:V为体积域;Sσ为边界域。

基于应变分解与非线性应变插值理论,有:

式中:δε′)为虚弹性应变增量;δε′′)为虚非线性应变增量;为插值点处的虚非线性应变增量;δa)为虚结点位移增量。

将式(15)代入虚功方程式(14),并结合式(11)所示的两种应力场关系以及单元的虚位移场可由虚结点位移表示,即:δu)=a),考虑到δa)与的任意性,整理可得隔离非线性法的控制方程。

式中:

式中:K为结构的初始刚度矩阵;K′为结构的系数矩阵,代表非线性应变与结点力之间的关系,两者在分析前集成得到且在整个计算过程中保持不变;K′′为包含结构材料非线性信息的分块对角矩阵,在计算过程中实时变化;ΔF为变形体结构的荷载增量。

对于线性六面体等参单元,单元矩阵的数值积分根据是否采用6点插值模型而选择6点直接积分或8点高斯积分。由于数值积分点与插值点重合,可以证明矩阵是以插值点为单位的列分块矩阵与分块对角矩阵,且集成后的矩阵K′与K′′具有同样的性质。

式中,分别为插值点i(i=1,2,···,nn=6或8)对应的权重、雅克比矩阵、应变矩阵以及弹塑性矩阵。

对于体积为Ve的四面体单元,单元矩阵采用Hammer数值积分,分别为:

2 控制方程求解

2.1 迭代求解法

隔离非线性法在求解时,需要判断所有非线性应变插值点状态,处于线弹性状态的插值点非线性应变增量为零,控制方程中的矩阵K′′中并不存在与之对应的元素,可消去矩阵K′以及K′′中对应元素。对于第m个增量步,控制方程式(16)可改写为:

基于Sherman-Morrison-Woodbury(SMW)[3]公式可将式(21)方程改写为:

式中:KSchur 为Schur补矩阵,计算如下:

式中:N为刚度矩阵的维数;M为隔离非线性法的计算自由度;Kf为结构刚度矩阵的逆矩阵。式(22)的求解原则为从右至左拆分为6个计算步依次进行[16],使得求解过程始终是矩阵与向量之间的运算,即:

隔离非线性法采用Newton-Raphson方法迭代求解,求解过程具有快速收敛特性,迭代格式如图3所示。第m个增量步的第n+1次迭代的不平衡力与位移增量改进值可分别由式(25)和式(26)计算得到,式(22)的求解只需更新矩阵K′′即可得到修正的Kf算式,避免了更新切线刚度矩阵并分解造成计算效率低下的问题。

2.2 组合近似法联合求解控制方程

式(22)的计算难点在于Schur补矩阵KSchur 与向量之间的运算,即式(24)中的第三步。KSchur 为满阵,当其规模较大时,直接分解高秩的Schur补矩阵计算量较大,且矩阵KTK-1K′预处理计算量大并需要较大的存储空间。

组合近似法[6-7,17](combined approximations approach,CA法)将局部近似法(泰勒级数展开)和全局近似法(多项式拟合法)相结合,利用局部效率上的优势和全局近似精度上的优势估计响应。采用CA法求解式(24)中的第三步,即式(27),求解过程主要是矩阵回代以及矩阵与向量之间的运算,避免了矩阵KSchur 的直接分解求解以及矩阵KTK-1K′的预处理计算。

图3 隔离非线性法的迭代格式
Fig.3 Iteration format of IS FEM

该过程可分为5个子步进行,具体计算流程如下。

1)构造初始基向量r1K′′为分块对角的稀疏矩阵,矩阵K′′的逆与向量Δx1的运算可视为各分块矩阵的逆与对应向量的运算,计算量较小。

2)构造其余的基向量ri以及基向量排列矩阵R=[r1 r2rs]。基向量计算过程中共涉及s次的矩阵回代、矩阵K′′的逆与向量运算以及矩阵K′(与KT)与向量的乘积,其中s为选取的基向量个数。基向量之间应保证线性无关须满足式(30),从而使得缩减矩阵KR非奇异,一般设定ε=0.95。

3)构造缩减矩阵KR。式(29)在求解基向量的过程中已获得矩阵K′(K-1(KR)),故在求解缩减矩阵KR时无须再次计算。

4)构造缩减向量FR

5)计算待求向量。

CA法联合求解控制方程,整个求解过程主要有s+2次的矩阵回代,s+1次的矩阵K′(与矩阵KT)与向量的乘积以及s+1次的矩阵K′′的逆与向量的运算。基于本文提出的隔离非线性实体单元模型以及基于CA法联合求解算法,编制了相应的计算程序,整体分析流程如图4所示。

图4 隔离非线性法的分析流程
Fig.4 Analysis flow of IS FEM

3 效率分析

一种算法经程序实现后其实际的计算效率受多种因素影响,如计算机硬件、分析软件以及程序代码的优化程度等。时间复杂度[15]指算法在具体计算时的工作量,是一种仅与算法本身相关而与环境变量无关的效率评价方法。基于时间复杂度理论对本文提出的CA法联合求解算法的工作量进行统计,并与传统变刚度法对比分析。

3.1 CA法联合求解算法的时间复杂度

CA法联合求解控制方程,只需在计算开始前对规模为N×N阶、带宽为m的刚度矩阵K进行一次LDLT分解,一次分解与回代的时间复杂度[15]约为:

式中:分解约为Nm2+6Nm+2N+(m2-1)/2,回代约为2Nm-N-(m2-1)2。

通过对CA法联合求解算法的工作量进行统计,其时间复杂度约为:

式中:系数A0B0C0见表3。表3为式(24)的6个计算步骤的时间复杂度,计算步骤3的时间复杂度见表4。其中,α为单元的结点自由度数目,线性四面体单元α=12,线性六面体等参单元α=24;β为矩阵 ′′K中小分块矩阵的规模,实体单元β=6。

表3 CA法联合求解算法的时间复杂度
Table 3 Time complexity of algorithm combined with CA

计算步骤 计算项目 时间复杂度1 1 Δ=Δ uKF-Nm-N-m-2(1)/2 2 1 2 T xKu(2α-1)M 3 T11 Δ=Δ 1′1 xKKKKx ⊗(见表4)4 3′2 Δx=KΔx(2α-1)M 5 1′′ ′Δ=-Δ--2()1 Δ=Δ uKx-2 2 3 2(1)/2 Nm-N-muuu N合计 2000 6 12 Δ=Δ+Δ TANBMC=++AmssCssmsm Bss 3222 0 0=+--=+-+-+=+++-++-+2(2)1 , [212(173)]/61,2[4(151)/3](91612)3 2αβ β ββ α 2 2 0

表4 计算步骤3的时间复杂度
Table 4 Time complexity of 3rd step

计算子步 计算项目 时间复杂度3.1 1 rKΔx 2(97)M/3 1=′′-1 β+β-2 3.2 1T1---+++-3.3 TT1 R{[()]}rKKKKr--(21)(1)/2[4(913)/3]Ms+′′ ′=-{[()]}mNsms i 1i αβ β 2 β+-3.4 T KRKRKKKR 22(2s2s)Ms′′ ′=--+-+-合计 000 FRx(2M-1)s 3.5 1 R =Δ1 Δx=RKF 32(97)/3(21)2(-)RR ssssM⊗=++DNEMF DmsFssms Ess 322 0 0 2=-=+-+=+++-+-+(21), [212(173)]/6,2[4(151)/3](109)/3 2 2 0 αββ β β

3.2 效率对比分析

一般的,刚度矩阵带宽可取由式(34)可知传统变刚度法一次迭代求解的时间复杂度约为:

由式(35)可知,CA法联合求解算法的时间复杂度与NM均相关,取s=10、α=24、β=6,此时联合求解控制方程一次迭代求解的时间复杂度约为:

T2=T1,可得M理论上限值Mmax(N),通过式(38)可得到非线性自由度的临界比例γmax。图5给出了临界比例γmax的变化规律,当非线性自由度比例γ=M/Nγmax时,CA法联合求解算法的时间复杂度均小于传统变刚度法的时间复杂度,且在极限状态下M一般不超过N的10倍,即γ<10。

图5 非线性自由度临界比例
Fig.5 Critical ratio of nonlinearity DOF

图6 不同γ的时间复杂度
Fig.6 Time complexity with different γ’s

图6为在不同γ值情况下本文方法与传统变刚度法的时间复杂度对比,当N>5000时,本文方法的时间复杂度均小于传统变刚度法,且N越大,本文方法的时间复杂度也就越低,即计算效率越高。

4 数值验证

4.1 悬臂柱静力非线性分析

图7为承受水平荷载作用的悬臂柱,悬臂柱长L=1.0 m,矩形截面边长H=0.1 m。模型采用两种实体单元进行网格划分,六面体单元共计划分10 000个单元以及12 221个结点,选择2×2×2点及两种6点这三种非线性应变插值方案;四面体单元共计划分73 733个单元以及14 619个结点。模型材料为双线性随动强化弹塑性模型,弹性模量E=201 GPa,泊松比v=0.3,切线模量Et=1 GPa,初始屈服应力σs=250 MPa。悬臂柱顶部某边施加横向均布线荷载作用P,采用多点位移控制算法[18],每个荷载步附加位移约束增量d0=0.1 mm,共计500个荷载步。

图7 悬臂柱与有限元模型
Fig.7 Cantilever column and finite element models

图8和图9分别为两种网格模型的悬臂柱顶部横向位移随荷载P的变化曲线,并与有限元分析软件ANSYS(相同单元类型、网格尺寸、单元数量、材料本构、荷载施加情况)的计算结果进行对比,两者吻合良好,顶端最大位移为50 mm时两类模型荷载P的相对误差均约为0.3%。这是由于CA法联合求解算法是一种近似而非精确的求解方法,较小的误差表明两类实体单元模型正确且本文算法具有足够的计算精度。

图8 六面体模型荷载-位移曲线
Fig.8 Load-displacement curves of hexahedron model

图9 四面体模型荷载-位移曲线
Fig.9 Load-displacement curves of tetrahedron model

图10和图11分别为隔离非线性法与传统变刚度法计算自由度与时间复杂度对比,可以看出即便隔离非线性的计算自由度大于传统变刚度法时,本文提出的CA法联合求解算法的时间复杂度也低于传统变刚度法的时间复杂度。

图10 计算自由度曲线
Fig.10 Calculation DOF curves

图11 时间复杂度曲线
Fig.11 Time complexity curves

表5为两类模型具体的计算时间统计,可以看出计算时间的大部分消耗在单元状态判定,方程求解与其他项占据较小部分,两种6点插值模型相较于2×2×2点插值模型单元状态判定效率提升了约20%。本文方法的方程求解时间远低于变刚度法的求解时间(采用同样软件编制变刚度法程序),通过方程求解时间比可以看出本文提出的CA法联合求解算法的计算效率约为传统变刚度法的6倍~8倍。

表5 计算时间统计
Table 5 Calculation time statistics

网格模型 维数 分解/s 回代/s 插值方案 总迭代步 总时间/s状态判定/s其他项/s方程求解/s迭代步数 方程求解/s刚度矩阵K 组合近似法联合求解算法 变刚度法 方程求解时间比/(%)6点方案1 1447 3329 2542(80.11%)12 77512.34 6点方案2 1474 3340 2593(81.72%)12 735 11.71 2×2×2 1410 3997 3173(100.0%)12 812 12.93四面体 43461 4.034 0.045 单点 1469 4169 3441 15 713 1084 4422 16.12六面体 36300 4.429 0.0441403 6276

4.2 拱坝动力非线性分析

式(39)为结构在地震作用下增量形式的运动方程:

式中:M为质量矩阵;C为阻尼矩阵;ΔF(t)为结构的恢复力增量;为地震动加速度增量向量;分别为结构速度增量以及加速度增量。

基于隔离非线性法理论,结构的恢复力:

将式(40)代入式(39),采用Newmark法求解该动力方程,并结合控制方程式(16)的第二项等式,可得隔离非线性法动力分析的控制方程:

式中:为等效荷载向量;为结构采用常阻尼矩阵的初始等效刚度矩阵。

采用上述理论对某混凝土拱坝进行地震作用下的动力反应分析。该拱坝为双曲拱坝,坝高240 m,坝顶厚13 m,坝宽667 m。坝体采用线性六面体实体单元进行网格划分,整个计算域共剖分为58044个单元和69321个结点,底部固定约束,有限元模型如图12所示。坝体的材料参数分别为:弹性模量E=30.0 GPa,泊松比v=0.2,密度ρ=2500 kg/m3,抗拉强度为2.9 MPa,抗压强度为24.1 MPa,采用理想弹塑性的Drucker-Prager模型模拟坝体的非线性效应。该分析模型仅考虑坝体在自重以及地震荷载作用下的反应,不考虑水压力、扬压力以及淤砂压力等对拱坝反应的影响。地震动采用y方向输入El Centro波,地震动强度参数PGA分别调整为0.3 g、0.4 g、0.5 g、0.6 g四种连续工况,阻尼比取0.05,时间间隔为0.02 s,共计40 s。在PGA=0.6 g的工况下,模型采用2×2×2点与两种6点共三种非线性应变插值方案,其余工况均采用2×2×2点插值方案。

图12 拱坝有限元模型
Fig.12 Arch dam finite element model

图13为四种工况下坝顶中点y方向位移时程曲线,本文的计算结果与ANSYS计算结果吻合良好,最大位移处相对误差均在3%左右。这是由于本文采用转移应力解析解法[19]对单元的本构关系进行积分,与ANSYS中处理方式略微不同,使得在大步长分析时两者计算结果存在较小差异。转移应力解析解法对小步长和大步长加载均有良好的收敛性,ANSYS分析时由于荷载步较长不易收敛,分析时需设置较多的子步。

图13 坝顶Y方向位移时程曲线
Fig.13 Y direction displacement-time curves in dam crest

该拱坝模型刚度矩阵维数为201354,刚度矩阵带宽优化后一次分解与回代的计算时间分别约为73.55 s、0.58 s,前者约为后者的127倍。图14、图15分别为隔离非线性与传统变刚度法的计算自由度与时间复杂度的时间对比,在17.36 s时隔离非线性法的计算自由度达到最大的1092792,相当于39.2%的单元进入非线性状态。CA法联合求解算法的时间复杂度约为3.46×109(基向量个数s=5),传统变刚度法的时间复杂度约为1.64×1011,两者之比约为2.1%。

图14 计算自由度曲线
Fig.14 Calculation DOF curves

图15 时间复杂度曲线
Fig.15 Time complexity curves

表6为4种工况下实际的计算时间统计,ANSYS分析时设置较多荷载子步,故在统计变刚度法的迭代步数时认为与隔离非线性法的迭代步数一致,求解时间认为只有在非线性迭代求解时考虑矩阵分解与回代,弹性求解只考虑矩阵回代。两者的方程求解时间均未考虑单元状态判定、其他项以及矩阵组装等。通过方程求解时间比可以看出,本文方法的计算效率约为传统变刚度法的5倍~12倍。在PGA=0.6 g工况下,三种插值方案所得的计算结果基本一致,两种6点插值方案相较于2×2×2点插值方案单元状态判定效率提高约25%。

表6 计算时间统计
Table 6 Calculation time statistics

工况 组合近似法联合求解算法 变刚度法 方程求解时间比/(%)插值方案 总时间/s 状态判定/s 其他项/s总迭代步 弹性步 迭代步 方程求解/s 方程求解/s 6点方案1 15967 13848(73.80%)184 2167 1851 316 19357.85 6点方案2 16022 13893(74.05%)184 2169 1849 320 1945 7.89 2×2×2 20928 18762(100.0%)190 2168 1850 318 1976 8.02 PGA=0.5 2×2×2 19556 17712 178 2105 1903 202 1666 16078 10.36 PGA=0.4 2×2×2 18383 16733 178 2063 1938 125 1472 10390 14.16 PGA=0.3 2×2×2 17580 16072 174 2037 1964 73 1334 6551 20.36 PGA=0.624646

5 结论

本文基于隔离非线性法基本理论,构造了线性四面体单元与线性六面体等参单元两种实体单元模型,其中线性六面体等参单元具有3种非线性应变插值方案。控制方程采用Woodbury公式与组合近似法联合求解,定量分析了该算法的时间复杂度并统计了实际计算时间。得到如下结论:

(1)线性六面体等参单元的两种6点插值方案相较于2×2×2点插值方案,可提高单元状态判定效率20%~25%,且三者计算结果基本一致。

(2)组合近似法联合求解控制方程,高秩的Schur补矩阵修改切线刚度同样具有较高的计算效率,且避免了Schur补矩阵中常矩阵的前处理计算。

(3)当结构的结点自由度数目大于5000时,组合近似法联合求解算法的时间复杂度均小于传统变刚度法,且结点自由度数目越大,计算效率也越高。

参考文献:

[1]Zienkiewicz O C, Taylor R L.The finite element method:Solid mechanics [M].5th ed.Oxford: Butterworth Heinemann Book Company, Inc, 2000: 1-21.

[2]Kirby R, Knepley M, Logg A, et al.Optimizing the evaluation of finite element matrices [J].Siam Journal on Scientific Computing, 2005, 37(3): 741-758.

[3]Akgün M A, Garcelon J H, Haftka R T.Fast exact linear and non-linear structural reanalysis and the Sherman-Morrison-Woodbury formulas [J].International Journal for Numerical Methods in Engineering, 2001, 50(7):1587-1606.

[4]Triantafyllou S P, Koumousis V K.Hysteretic finite

elements for the nonlinear static and dynamic analysis of structures [J].Journal of Engineering Mechanics, 2014,140(6): 04018008-1-04018008-17.

[5]Li ZG, Wu B S.A preconditioned conjugate approach to structural reanalysis fo rgeneral layout modifications [J].International Journal for numerical Methods in engineering, 2007, 70(5): 505-522.

[6]Kirsch U.A unified reanalysis approach for structure analysis, design, and optimization [J].Structural &Multidisciplinary Optimization, 2003, 25(2): 67-85.

[7]Kirsch U.Reanalysis of structures: A unified approach for linear, nonlinear, static and dynamic systems [M].Netherlands Springer, 2008.

[8]Li Gang, Yu Dinghao.Efficient inelasticity-separated finite element method for material nonlinearity analysis[J].Journal of Engineering Mechanics, 2018, 144(4):04018008-1-04018008-11.

[9]Li Gang, Yu Dinghao, Li Hongnan.Seismic response analysis of reinforced concrete frames using inelasticity-separated fiber beam-column model [J].Earthquake Engineering & Structural Dynamics, 2018,47(5): 1291-1308.

[10]李钢, 余丁浩, 李宏男.基于拟力法的纤维梁有限元非线性分析方法 [J].建筑结构学报, 2016, 37(9): 61-68.Li Gang, Yu Dinghao, Li Hongnan.Nonlinear fiber beam element analysis based on force analysis method [J].Journal of Building Structures, 2016, 37(9): 61-68.(in Chinese)

[11]王勖成.有限单元法[M].北京: 清华大学出版社,2003: 117-122.Wang Xucheng.Finite element method [M].Beijing:Tsinghua University Press, 2003: 117-122.(in Chinese)

[12]Ueda Y, Fujikubo M.Generalization of the plastic node method [J].Computer Methods in Applied Mechanics and Engineering, 1991, 92(1): 33-53.

[13]牛辉, 汪劲丰, 张巍, 等.基于实体退化单元的高墩非线性稳定仿真分析[J].浙江大学学报(工学版), 2012,46(6): 1082-1089.Niu Hui, Wang Jinfeng, Zhang Wei, et al.Simulation anslysis of nonlinear stability of high pier based on degenerated solid element [J].Journal of Zhejiang University(Engineering Science), 2012, 46(6): 1082-1089.(in Chinese)

[14]Dhatt G, Touzot G, Lefrancois E.Finite element method[M].London and Hoboken: STE Ltd and John Wiley &Sons, Inc, 2012: 368-369.

[15]Golub G H,Van Loan C F.Matrix computations [M].4th ed.Beijing: Posts & Telecom Press, 2014: 176-185.

[16]李钢, 贾硕, 余丁浩.基于算法复杂度理论的拟力法计算效率评价 [J].计算力学学报, 2018,35(2): 129-137.Li Gang, Jia Shuo, Yu Dinghao.The efficiency evaluation of force analogy method based on the algorithm complexity theory [J].Chinese Journal of Computational Mechanics, 2018, 35(2): 129-137.(in Chinese)

[17]王琥, 种浩, 高国强, 等.重分析方法研究进展及展望[J].工程力学, 2017, 34(5): 1-15.Wang Hu, Chong Hao, Gao Guoqiang, et al.Review of advances and outlook in reanalysis methods [J].Engineering Mechanics, 2017, 34(5): 1-15.(in Chinese)

[18]黄羽立, 陆新征, 叶列平, 等.基于多点位移控制的推覆分析算法 [J].工程力学, 2011, 28(2): 18-23.Huang Yuli, Lu Xinzheng, Ye Lieping, et al.A pushover analysis algorithm based on multiple point constraints[J].Engineering Mechanics, 2011, 28(2): 18-23.(in Chinese)

[19]杨强, 杨晓君, 陈新.基于D-P准则的理想弹塑性本构关系积分研究 [J].工程力学, 2005, 22(4): 15-19.Yang Qiang, Yang Xiaojun, Chen Xin.On integration algorihms for perfect plasticity based on Drucker-Prager criterion [J].Engineering Mechanics, 2005, 22(4): 15-19.(in Chinese)

THE INELASTICITY-SEPARATED SOLID ELEMENT MODEL AND COMPUTATIONAL EFFICIENCY ANALYSIS

LI Jia-long1 , LI Gang1 , LI Hong-nan1,2
(1.State Key Laboratory of Costal and Offshore Engineering, Dalian University of Technology, Dalian 116024, China;2.College of Civil Engineering, Shenyang Jianzhu University, Shenyang 110168, China)

Abstract: Due to the large number of elements required in the calculation of solid finite element models, large computational resource is consumed in the element state determination process and the factorization of the global tangent stiffness matrix with large dimensions, thus resulting in low efficiency.In this paper, linear tetrahedron and hexahedron isoparametric elements are established based on the inelasticity-separated finite element method.Six direct integration points are considered for hexahedron elements as the nonlinear strain interpolation points instead of eight gauss integration points, of which the computational accuracy is stable and the efficiency is improved.In addition, the main solving process of the governing equation is only the back substitution of the initial stiffness matrix and the matrix-vector multiplication by using the Woodbury formula and combined approximation approach.Therefore, the efficiency is significantly improved.Finally, the computational efficiency of the proposed method based on the time complexity theory indicates that, with the increase of the number of nodal degrees of freedom, the computational efficiency of the proposed method is significantly improved as compared with the traditional variable stiffness method.The numerical examples verify the correctness of the proposed solid element model and the high efficiency of the proposed algorithm.

Key words: inelasticity-separated finite element method; solid element; direct integration; Woodbury formula;combined approximations approach; efficiency analysis

中图分类号:TU375.6;TU311.4

文献标志码:A

doi: 10.6052/j.issn.1000-4750.2018.07.0383

文章编号:1000-4750(2019)09-0040-10

收稿日期:2018-07-13;修改日期:2019-01-04

基金项目:国家自然科学基金项目(51878112);中央高校基本科研业务费专项资金项目(DUT17ZD220);大连市高层次人才创新支持计划项目(2015R044,2017RD04)

通讯作者:李 钢(1979―),男,辽宁葫芦岛人,教授,博士,博导,主要从事结构工程抗震等研究(E-mail: gli@dlut.edu.cn).

作者简介:

李佳龙(1991―),男,四川遂宁人,博士生,主要从事结构非线性分析等研究(E-mail: wekalee@163.com);

李宏男(1957―),男,辽宁沈阳人,教授,博士,博导,主要从事结构工程抗震和结构健康检测等研究(E-mail: hnli@dlut.edu.cn).