基于S-R和分解定理的無網格Galerkin法求解幾何非線性問題

摘 要：針對幾何非線性問題，基于S-R和分解定理與更新拖帶坐標描述法，從運動變換的角度出發，由虛功率原理建立全局弱式積分方程.采用無網格Galerkin法(EFG)求得該問題的數值解.對于非線性問題的數值計算，通常采用增量和迭代法.S-R和分解定理給出了一種新的應變張量，該應變張量能合理地反映大位移大轉動情況下物體的應變狀態及轉動情況，而且有利于求解幾何非線性問題.數值算例驗證了該方法的合理性和有效性.

關鍵詞：S-R和分解定理；無網格Galerkin法；幾何非線性；更新拖帶坐標法

中圖分類號：O343.5 文獻標識碼：A

Element Free Galerkin Method for Geometrically Nonlinear Problems Based on the S-R Decomposition Theorem

CHEN Fang-zu， LUO Dan

(State Key Laboratory of Advanced Design and Manufacturing for Vehicle Body， Hunan Univ， Changsha， Hunan 410082，China)

Abstract:Aiming at geometrically nonlinear problems and using movement transformation as a point of departure， a global weak form integral equation was established from the virtual power principle on the basis of the S-R decomposition theorem and updated co-moving coordinate formulation. Numerical solutions were obtained by using element free Galerkin method (EFG). Usually， incremental and iterative methods are used to solve the nonlinear problems of numerical calculation. S-R decomposition theorem brings out a new strain tensor， which can reasonably reflect the state of stain and rotation on the condition of large displacement and large rotation， and are helpful in solving geometrically nonlinear problems. Numerical examples have shown that this method is reasonable and effective.

Key words:S-R decomposition theorem; element free Galerkin method; geometrical nonlinearity; updated co-moving coordinate formulation

無網格法是一種直接通過節點信息來離散問題域的數值計算方法.與傳統的有限元法相比，它節省了前期為離散問題域而劃分單元的成本，具有很好的自適應性，其計算結果通常是光滑連續的，因此一般不需要應力均勻等后處理.無網格的近似函數是通過節點直接構建的，不依賴于網格，抗畸變能力強，在求解幾何非線性問題時具有很大的優勢.無網格Galerkin法[1-2]是一種常用的全局弱形式無網格法，它使用移動最小二乘近似來構造形函數，通過Galerkin弱勢來建立離散系統方程，得到的系統矩陣是帶狀稀疏矩陣且具有對稱性.這種方法在精確性、收斂性以及效率方面更優于傳統有限元.目前，無網格Galerkin法在解決幾何非線性、材料非線性和板殼結構等問題上得到了廣泛的應用.

對于幾何非線性問題，傳統的幾何非線性理論[3-4]都是基于格林應變張量和Kirchhoff應力張量建立起來的，格林應變張量只是能夠反映物體運動的相對位移，而無法反映物質點在形變過程中的相對轉動情況，在數學上講格林張量是不完備的，而且其本身也存在一些缺陷.

S-R和分解定理[5]從運動變換的角度出發，推導出新的應變張量，這種張量既反映物體的相對位移，又反映物體內各個物質點的轉動情況，是一種合理而又精確的應變描述.更新拖帶坐標[6-7]是指在t時刻的所有物理量都以t時刻的初始拖帶系的度規張量表示，參考位形隨著增量步的不斷變化而更新.S-R和分解定理在有限元分析中已經得到很好的應用[8].本文基于新的應變張量，從虛功率原理[9]出發建立更新拖帶坐標法(U.C)下全局弱式控制方程，采用無網格Galerkin法離散系統方程，最后采用增量法和Newton-Raphson迭代結合進行數值求解.

1 S-R和分解定理

S-R和分解定理：“給定物體一個可能的位移函數，此函數在變形點集域內是單值連續的，處處具有一階導數，則此運動變換總可以唯一分解為正交與對稱兩個子變換的直和.正交變換表現為點集的轉動，對稱變換表現點集的形變.”［5］即

Fij=δij+uij=Sij+Rij，(1)

Sij=12(ui|j+ui|Tj)-LikLkj(1-cosθ)，(2)

Rij=δij+Lijsin θ+LikLkj(1－cos θ)， (3)

θ=±arcsin ［－ωijωji］1/2，(4)

Lij=ωij/sin θ， (5)

ωij=(ui|j-ui|Tj)/2.(6)

式中：Fij為運動變換張量分量；Sij為應變張量分量；Rij為轉動張量分量；θ表示物質點的平均整旋角；Lij為轉軸方向余弦.S-R和分解定理表明物體在運動過程中轉動和形變是同時進行的，不能把它看作轉動和形變的簡單疊加.該定理的唯一性、存在性和客觀性已得到證明.在微小變形情況下，θ的值很小，此時應變張量Sij也就是在小變形情況下常用的Cauchy應變張量.

2 更新拖帶坐標法下的增量變分方程

在非線性大變形過程中，應力與外力呈非線性關系.對于非線性問題一般采用增量方程求解，通常采用時間參量t來描述問題的外載荷和應力狀態的關系.

假設從0到t時刻的所有時間點的力學量已經求得，那么下一步要求t+Δt的各個力學量.由虛功率變分原理，t+Δt時刻勢能率變分為：

t+ΔtδJ=∫t+ΔtΩ(t+Δt)σijδt+ΔtVj‖idΩ-t+ΔtQ.(7)

其中

t+ΔtQ=∫Γt(t+Δt)piδt+ΔtVidΓ+

∫t+ΔtΩ(t+Δt)ρfiδt+ΔtVidΩ.(8)

式中:t+ΔtVj‖i表示在t+Δt時刻的速度t+ΔtVj關于t+Δt時刻的拖帶系協變基矢t+Δtgi的協變導數.假設在t時刻各個物理量已經得到，那么對于在t+Δt時刻的位形，由泰勒展開式，得到

t+ΔtVk‖i=(δkl-Δuk‖j)tVj‖i，(9)

t+Δtσij=tσij+Δtij.(10)



將式(9)和(10)代入變分方程(7)，化簡得到

Δt∫tΩ(t)ijδtVj‖idΩ-Δt∫tΩ(t)σijVj‖kδtVk‖idΩ=

t+ΔtQ-∫tΩ(t)ijδtVj‖idΩ.(11)

式(11)右端最后項作為虛擬荷載，在求解過程中與外載荷一起進行平衡迭代.為求解方程(11)，采用t時刻初始拖帶系t0gi內的速度分量t0Vi作為基本插值量.取t時刻初始拖帶系為標準直角坐標系，把t時刻瞬時拖帶系tgi上的tVi‖j轉換到初始拖帶系t0gi上：

tVi‖j=t0Vk|lxitxktxlxj.(12)

對于t0Vk|l，由廣義歐拉運動公式可表示為

t0Vi|j=t0ij+t0ij.(13)

t0ij表現為物體的應變速率，t0ij表現為物體的轉動角速度與方位角余弦的乘積.即

t0ij=12(t0Vi|j+t0Vi|Tj)，(14)

t0ij=Lij＝12(t0Vi|j-t0Vi|Tj).(15)

由速率型物性方程，在t時刻的瞬時拖帶系tgi上，

tij=tDik tjllk.(16)

由張量的性質，應力速率滿足坐標變換.在t時刻的初始拖帶系t0gi上，

t0ij=t0Dikjlt0lk.(17)

其中，

t0Dikjl=tDmrnsxntxjtxixmxstxltxkxr.(18)

對于各向同性材料，

t0Dikjl=tDmrns.(19)

把式(17)代入式(11)，得到

Δt(∫(Ω(t)t0Dikjlt0lkδt0jidΩ+∫Ω(t)t0ijδt0jidΩ)-

Δt(∫(Ω(t)t0σijt0jkδt0kidΩ+Δt∫Ω(t)t0σijt0jkδt0kidΩ+

∫Ω（t)t0σijt0jkδt0kidΩ+∫Ω(t)t0σijkiδt0jkdΩ)=

t+ΔtQ-∫Ω(t)t0σijδt0jkdΩ.(20)

對于式(20)左端的第2項，由應力張量的對稱性可以證明在t時刻初始拖帶系下，

t0ij-t0σikt0kj+t0ikt0σkj=0.(21)

將式(21)代入(20)得到

Δt∫Ω(t)t0Dikjlt0lkδt0jidΩ-Δt(∫Ω(t)t0σijt0jkδt0kidΩ+

∫Ω(t)t0σijt0jkδt0kidΩ+∫Ω(t)t0σijt0jkδt0kidΩ+

∫Ω（t)t0σijkiδt0jkdΩ)=t+ΔtQ-∫Ω(t)t0σijδt0jidΩ.(22)

此式作為無網格法數值計算的基礎積分方程.

3 更新拖帶坐標法無網格Galerkin法格式

對于節點位移增量ui和速度Vi，采用移動最小二乘插值形函數進行近似:

ui=∑nk=1φkuik，Vi=∑nk=1φkVik. (23)

式中:φk為節點k的移動最小二乘法(MLS)的形函數; uik為節點位移; Vik為節點速度; ui為計算點的近似位移，Vi為計算點近似速度.

把式(23)代入式(22)，離散后寫成矩陣形式，有

［tKL-tKN+Kα］ΔU=t+ΔtQ-tF+Fα.(24)

各對應矩陣表達式：

tKL=∫Ω(t)BsTDBsdΩ，(25)

tKN=∫Ω(t)BsTσ1BsdΩ+

∫Ω(t)BrTσ2BrdΩ+∫Ω(t)BsTσ3BrdΩ+

∫Ω(t)BrTσ3TBsdΩ，(26)

t+ΔtQ=∫Γt(t+Δt)ΦTpdΓ+∫Ω(t+Δt)ΦTρfdΩ. (27)

tF=∫Ω(t)BsTσdΩ，(28)

Kα=－∫ΓuΦTαΦdΓ，

Fα=－∫ΓuΦTαdΓ.(29)

式中：Kα和Fα分別為總體懲罰剛度矩陣和總體懲罰力向量；Γu表示問題的位移邊界；t0D為物性矩陣.對于二維平面幾何非線性問題， 各矩陣和向量的具體表達式為:

Bs={Bs1，Bs2，…，Bsn}，(30)

Br={Br1，Br2，…，Brn}，(31)



Bsk=Φkx1+Γ111ΦkΓ112ΦkΓ221ΦkΦkx2+Γ222Φk

Φkx2+(Γ121+Γ211)ΦkΦkx1+(Γ122+Γ212)Φk，(32)

Brk=Φkx2+(Γ121－Γ211)ΦkΦkx1+(Γ122－Γ212)Φk，(33)

σ=［σ11，σ22，σ21］T，(34)

σ1=σ110σ1220σ22σ122σ122σ122σ11+σ224，(35)

σ2=－σ11+σ224，(36)

σ3=－σ122，σ122，σ11－σ224T.(37)

式中：Γijk為第二類Christoffel符號，當取初始坐標系與直線直角坐標系同胚，Γijk=0Γijk=0，因此在更新拖帶坐標法求解幾何非線性問題中的Bs與求解線性問題中的應變矩陣B相同.

4 數值算例

算例1 分析一大變形的懸臂梁(如圖1所示），取單位厚度，梁的尺寸為25 m×2 m.在梁的自由截面端的中點處施加一集中載荷，材料參數為ν=0.3，E=3.0×104 N/cm2.這里假設材料為線彈性的.在用無單元Galerkin法計算時，采用31×5個規則分布的節點來描述該梁的初始域.移動最小二乘(MLS)形函數插值的基函數選取6個單項式，采用三次樣條函數作為權函數，取權函數影響域尺寸半徑r=2.5 cm.對于整個加載過程采用20個增量步.對于每一步通過標準Newton-Raphson迭代法求解，計算收斂很快，其中每個荷載增量步的迭代次數都不超過5.該懸臂梁端點的撓度計算結果如圖2所示.由圖2可以看出在大位移大變形情況下，懸臂梁的撓度相比線性分析呈現剛化現象，而且無網格Galerkin法的數值解與有限元的結果也吻合得很好.

算例2 圖3所示的受均布荷載的簡支梁，取單位厚度，L=6 m， μ=0.3， E=3×107Pa，梁的上表面受到均布荷載P=2.0×105 N/m，假設材料為線彈性，不計體力.在用無網格Galerkin法求解時，采用6×16個規則節點來離散問題域.用移動最小二乘法構建形函數，采用四次樣條函數作為權函數，權函數的半徑r=0.8 m.在分布力邊界上采用4個高斯積分點.載荷始終垂直于梁的頂面和底面.采用20個增量步，在每個增量步內采用Newton-Raphson法迭代求解.得到簡支梁的變形情況如圖4所示，各節點的位移見表1，與ANSYS的計算結果比較可見，無網格Galerkin法的計算y結果是合理而有效的.

5 結 論

本文將基于S-R和分解定理的更新拖帶坐標法運用于求解幾何非線性問題，通過數值算例證明更新拖帶坐標理論運用在求解大位移大轉動問題上是合理而有效的，即使在大轉動情況下仍可采用大的增量步.由于該理論所推導出的控制方程是以張量表述的，因此可以進行張量變換，滿足廣義量綱原理，具有一定的普遍意義.這一理論同樣可以進一步運用于求解材料非線性、板殼結構問題，理論上講應該能夠得到更為精確的結果.無網格法在求解非線性問題上體現了較之有限元更高的精度和收斂性，同樣也存在一些缺點，比如計算成本較高等問題，這也是需要進一步解決的.

參考文獻

［1］ BELYTSCHKO T，LU Y Y，GU L．Element free Galerkin methods[J]．International Journal for Numerical Methods in Engineering，1994，37(2)：229-256．

[2] LIU G R， GU Y T. An introduction to meshfree methods and their programming[M]. Berlin: Springer， 2005: 161-236.

[3] 熊淵博，龍述堯，胡德安. 幾何非線性問題的無網格法分析[J]. 機械強度，2006，28(1)：083-087.

XIONG Yuan-bo， LONG Shu-yao， HU De-an. Analysis of the geometrically nonlinear problem by meshless method [J]. Journal of Mechanical Strength，2006，28(1)：083-087. (In Chinese)

[4] 王勖成. 有限單元法[M]. 北京：清華大學出版社， 2003: 618-629.

WANG Xu-cheng. Finite element method[M]. Beijing:Tsinghua University Press， 2003: 618-629.(In Chinese)

[5] 陳至達. 理性力學[M].重慶：重慶出版社，1999:128-134.

CHEN Zhi-da. Rational mechanics[M]. Chongqing: Chongqing Press， 1999:128-134.(In Chinese)

[6] 李平. 非線性大變形有限元分析的更新拖帶坐標法[D].徐州：中國礦業大學，1991.

LI Ping. The updated co-moving coordinate formulation for the nonlinear large deformation finite element analysis and application[D]. Xuzhou: China University of Mining and Technology， 1991.(In Chinese)

[7] LI Ping， CHEN Zhi-da. The updated co-moving coordinate formulation of continuum mechanics based on the S-R decomposition theorem [J]. Computer Methods in Applied Mechanics and Engineering， 1994， 114:21-34.

[8] 高立堂， 宋玉普， 董毓利. 火災下鋼筋混凝土板的熱彈塑性有限元分析——基于S-R 分解原理[J]. 計算力學學報，2007， 24(1): 86-90.

GAO Li-tang， SONG Yu-pu， DONG Yu-li. Thermal-elastic-plastic finite element analysis of reinforced slabs under fire based on S-R decomposition theorem[J]. Chinese Journal of Computational Mechanics， 2007，24(1): 86-90.(In Chinese)

[9] 秦忠. 基于新的大變形理論的非線性有限元及其應用[D].徐州：中國礦業學院，1986.

QIN Zhong. The nonlinear large deformation finite element analysis and application—base on the new large deformation theorem[D]. Xuzhou: China College of Mining and Technology， 1986.(In Chinese)