基于牛頓-拉夫遜的拉壓不同模量問題的數(shù)值求解

*，，，

(沈陽工業(yè)大學 機械工程學院，沈陽 110178)

1 引 言

經(jīng)典彈性理論的主要研究對象是具有相同拉伸和壓縮彈性行為的問題。然而，在實際工程中，許多材料與結構都在不同程度上表現(xiàn)出拉壓不同的彈性性質。材料方面，如混凝土、石墨、玻璃鋼和一些高分子聚合物[1-4]等，均具有拉壓不同模量；除了材料本身具有拉壓不同模量的特性外，某些結構自身的構型特征決定了即便使用拉壓相同模量的材料，也會具有不同拉壓性質，如已經(jīng)在土建及航天航空等領域得到了廣泛應用的索和膜結構[5-7]。當這類問題的拉壓性質相差大到一定程度時，若仍采用經(jīng)典的彈性理論和有限元方法進行分析計算，會產(chǎn)生誤差甚至是錯誤。綜上所述，針對這類典型問題進行準確計算分析及其精密設計時，迫切需要發(fā)展基于有限元的高效數(shù)值求解方法。

Ambartsumyan[8]將拉壓不同模量材料的本構關系總結為雙直線模型，系統(tǒng)地介紹了該問題的初始假定與基本方程，在此基礎上提出了拉壓不同模量的彈性理論[8,9]。近年來，由于工程問題求解以及理論發(fā)展的需要，不僅解析法[1,10-13]成功地應用到某些典型問題中，而且基于有限元法的數(shù)值求解技術也得到了相應的發(fā)展。其中，楊海天等[14，15]提出了基于初應力技術的有限元迭代格式，并建立了基于連續(xù)域蟻群算法的二維拉壓不同模量反問題的數(shù)值求解模型；He等[16]利用一個平面應力數(shù)值算例，測試了一種針對拉壓雙模量問題的加速收斂方法；張洪武等[5，17]基于參變量變分原理與共旋法，建立一種高效的求解策略，用于雙模量桁架結構與張拉整體結構幾何非線性行為的分析。Ding等[6,7]同樣基于參變量變分原理，分析了薄膜起褶這一現(xiàn)象。

拉壓不同模量問題究其根本為材料本構方程的不連續(xù)性，從而導致此問題具有強非線性。牛頓-拉夫遜格式因其具有收斂速度快、收斂穩(wěn)定及編程簡單等優(yōu)點，已在ANSYS和ABQUS等大型有限元軟件中成功地應用于很多非線性問題的求解，無論在理論研究還是工程應用方面，均取得很好的效果。因此，本文為了快速準確地求解拉壓不同模量問題，借鑒文獻[18]的策略，采用改進的Heaviside函數(shù)將材料本構方程連續(xù)光滑化，并推導連續(xù)化后的切線剛度矩陣，最后基于牛頓-拉夫遜格式進行求解，并將相應的算法應用于大規(guī)模計算問題的分析與設計中。

2 拉壓不同模量問題的數(shù)值求解

2.1 拉壓不同模量問題的介紹

拉壓不同彈性模量問題中，主應力和主應變的關系矩陣可表示為

(1)

根據(jù)拉壓不同模量的基本原理，柔度矩陣A的元素可表示

(2)

式中Et,μt，Ec和μc分別代表拉壓彈性模量和相應的泊松比，下標t代表拉，c代表壓，且假設μt/Et=μc/Ec。不難看出,拉壓不同模量問題的本構方程實質是一個非連續(xù)階梯函數(shù),為了便于采用牛頓-拉夫遜格式進行求解，本文采用如下改進的Heaviside函數(shù)來光滑連續(xù)化這一階梯函數(shù)。

(3)

(4)

圖1 采用改進Heaviside函數(shù)近似本構方程

Fig.1 Approximated material behavior using the modified Heaviside function

式中 全應力矩陣σM可表示為

(5)

式中 相應的特征向量{l1,l2,l3}T,{n1,n2,n3}T和{m1,m2,m3}T分別為三個主方向角{α,β,γ}與初始笛卡爾直角坐標軸{x,y,z}的余弦。

與應力矩陣相對應的全應力向量為

σ= {σx x,σy y,σz z,σx y,σy z,σz x} =

{σ1,σ2,σ3,σ4,σ5,σ6}

(6)

σ=TσI

(7)

相應地，全應變向量可表示為

ε= {εx x,εy y,εz z,εx y,εy z,εz x} =

{ε1,ε2,ε3,ε4,ε5,ε6}

(8)

εI=TTε

(9)

因此，將式(1,9)代入式(7)可得，全應力向量和全應變向量之間的關系為

σ=TD0TTε=D(σI)ε

(10)

式(7,9)的轉軸矩陣T定義為

(11)

式中l(wèi)i,mi和ni(i=1,2,3)分別為式(4)特征向量的元素。

2.2 切向剛度矩陣

拉壓不同模量情況下，應力的計算式(10)的張量形式為

(12)

式中 自由指標i=1,2,…,6，啞指標j=1,2,…,6，k,l=1,2,3。推導切線剛度張量，即推導式(12)關于應變張量的導數(shù)，結果如下，

(13)

第一項?D0k l/?εm可表示為

(14)

第二項?Tli/?εm可表示為

(15)

(i= 1,2,…,6) (16)

2.3 牛頓-拉夫遜求解格式

基于有限元方法的基本原理，拉壓不同模量問題的非線性方程組可表示為

K(U)U=F

(17)

式中K(U)為結構整體非線性剛度矩陣，U為位移向量，F(xiàn)為外載荷向量。

根據(jù)牛頓-拉夫遜迭代格式，方程(17)的求解格式可表示為

ΔU=-K-1t(P-F)

(18)

式中 ΔU為節(jié)點位移增量向量，P為結構整體節(jié)點內力向量，Kt為結構整體的切線剛度矩陣。整體切線剛度陣與整體節(jié)點內力向量分別由單元切線剛度陣與單元節(jié)點內力矢量組集而成，具體計算列式如下，

(19)

式中Ket為單元切線剛度矩陣，Pe表示單元節(jié)點內力矢量，B為單元應變矩陣，ΔD為當前材料切線剛度矩陣。ΔD可依據(jù)式(13～16)進行計算，D為當前材料剛度矩陣，D可依據(jù)式(10)進行計算。

牛頓-拉夫遜迭代格式為

(2) 計算結構整體剛度矩陣K，形成結構外載荷向量F，求解KU=F。

(4) 判斷‖ΔU‖ <κ是否成立，成立則進行步驟(5)；否則，U=U+ΔU，回到步驟(3)。

3 算 例

算例1受水平載荷的平面應力問題

選取文獻[17]的經(jīng)典算例以驗證本文方法的有效性。該算例屬于平面應力問題，且本文均采用二維平面問題作為算例。本文所有計算公式均以三維空間問題為研究對象，只需將計算式中三維問題的矩陣與向量更換為二維問題的矩陣與向量即可。

為了形成鮮明的對比，本文采用與文獻[16]相同的幾何和物理參數(shù)以及有限元網(wǎng)格參數(shù)設置等，如圖2所示，抗拉模量Et=2.2 GPa，抗壓模量Ec=3.32 GPa，相應的泊松比分別為μt=0.22與μc=0.332，以滿足μt/Et=-μc/Ec。

圖2 受水平載荷的平面應力問題

Fig.2 Plane -stress problem under the horizontal load

分別采用文獻[17]與本文方法求解圖2所示有限元模型，兩個方向的位移和各單元主應力的計算結果分別列入表3和表4。可以看出，本文所得解與文獻[17]吻合，從而證明了本文方法的正確性。表5 給出了在選取不同控制誤差時，文獻[17]中所述的求解策略和本文方法所需的迭代次數(shù)?？梢钥闯觯疚姆椒ň哂辛己玫氖諗啃?，且效率較高。雖然采用本文方法求解此算例所需迭代終止步數(shù)稍多于文獻[17]的方法，但本文方法每一迭代步只需求解線性方程組即可，而文獻[17]需要求解計算規(guī)模較大的線性互補問題，因此本文方法更加適用于大規(guī)模問題的求解。

ξ-迭代次數(shù)ux(本文)ux(文獻[17])1×10-12-0.283981×1021×10-23-0.324754×1021×10-33-0.376413×102-0.425152×1021×10-44-0.406229×1021×10-55-0.415338×1021×10-66-0.425147×102

表2 迭代參數(shù)

Tab.2 Iteration parameters

ξ(初始)ξ-Z11×10-620

表3 節(jié)點位移Tab.3 Displacements of nodes

表4 單元主應力和主方向Tab.4 Principal stress and principal direction of elements

表5 不同控制誤差下求解所需迭代次數(shù)

Tab.5 Iterative numbers for different control error

κ迭代次數(shù)文獻[17]本文κ迭代次數(shù)文獻[17]本文1×10-1361×10-711155×10-2461×10-813151×10-2581×10-914165×10-36101×10-1016161×10-36111×10-1117181×10-59141×10-1218191×10-610141×10-132022

算例2拓撲優(yōu)化設計(單純受壓)

以計算規(guī)模最大拓撲優(yōu)化設計問題為例，驗證本文方法計算效率高和求解速率快的特性。

如圖3所示，60 m×71 m的矩形結構作為初始設計區(qū)域，4個角點均約束固定。根據(jù)結構幾何形式、載荷分布情況以及邊界條件的對稱性，本文設計時只取右半邊作為設計域，中間一層單元作為非設計區(qū)域，用于承受均布的壓力載荷1 kN/m。設計目標函數(shù)采用結構柔順性最小化，以獲得剛度最大化的設計，材料用量不超過整個設計域15%；抗壓模量Ec=210 GPa，泊松比μc=0.3，為了避免數(shù)值奇異性，令Ec=Et/106，因此μc=μt/106，且迭代求解參數(shù)仍選用表2的數(shù)據(jù)。

如圖4所示，單純受壓結構的拓撲在設計過程中的演化，僅需45步即可得到非常清晰的拓撲形式?？梢钥闯?，由于采用了單純受壓的材料模型，在整個設計過程中承載面之上的區(qū)域材料基本為空，幾乎所有的材料都布置在承載面之下來承受壓力。最優(yōu)的拓撲形式如圖5所示，設計結果為經(jīng)典的下承式拱橋結構，與文獻[20]中設計結果的拓撲形式基本相同。且根據(jù)表6的數(shù)據(jù)對比可知，本文所提出的方法很大程度上提高了計算效率與求解精度。

圖3 初始設計區(qū)域

Fig.3 Initial design domain

圖4 單純受壓結構拓撲隨迭代步演化過程

Fig.4 Iteration histories of structural topology with compression-only material

圖5 單純受壓結構最優(yōu)拓撲

Fig.5 Optimal topologies of the structure with compression-only material obtained by using different methods

表6 不同方法所獲得的計算結果

Tab.6 Results obtained by different method

方法迭代次數(shù)單步運算時間/s目標函數(shù)本文452816.84e-2文獻604267.27e-2

4 結 論

本文基于牛頓-拉夫遜求解格式，提出了拉壓不同模量有限元求解方法，主要策略在于光滑連續(xù)化了本構方程，在此基礎上利用特征值與特征向量的求導方法推導了切線剛度矩陣的列式。數(shù)值算例證明了本文方法的正確性，并采用拓撲優(yōu)化設計算例進一步驗證了本文方法在穩(wěn)定性和效率上的優(yōu)勢。

參考文獻(References):

[1] Yang Q,Zheng B L,Zhang K,et al.Elastic solutions of a functionally graded cantilever beam with different modulus in tension and compression under bending loads[J].AppliedMathematicalModelling,2014,38(4):1403-1416.

[2] Khan A H,Patel B P.Nonlinear forced vibration response of bimodular laminated composite plates[J].CompositeStructures,2014,108:524-537.

[3] 魏 靖，石多奇，孫燕濤，等.拉壓不同模量的縫合三明治夾芯結構梁彎曲性能[J].復合材料學報,2015,32(1):160-166.(WEI Jing,SHI Duo -qi,SUN Yan-tao,et al.Flexural properties of stitched sandwich structure beam with different modulus in tension and compression[J].ActaMateriaeCompositaeSinica,2015,32(1):160-166.(in Chinese))

[4] Iraola B,Cabrero J M.An algorithm to model wood accounting for different tension and compression elastic and failure behaviors[J].EngineeringStructures,2016,117:332-343.

[5] Zhang L,Gao Q,Zhang H W.An efficient algorithm for mechanical analysis of bimodular truss and tensegrity structures[J].InternationalJournalofMechanicalSciences,2013,70:57-68.

[6] Ding H L,Yang B G,Lou M,et al.New numerical method for two -dimensional partially wrinkled membranes [J].AIAAJournal,2003,41(1)：125-132.

[7] Ding H L,Yang B G.The modeling and numerical analysis of wrinkled membranes [J].InternationalJournalforNumericalMethodsinEngineering,2003,58：1785-1801.

[8] 阿姆巴爾楚米楊.不同模量彈性理論[M].張允真，譯.北京：中國鐵道出版社,1986.(Ambartsumyan S A.DifferentModulusofElasticityTheory[M].ZHANG Yun-zhen,translated.Beijing:China Railway Publishing House,1986.(in Chinese))

[9] Sun J Y,Zhu H Q,Qin S H,et al.A review on the research of mechanical problems with different moduli in tension and compression[J].JournalofMechanicalScienceandTechnology,2010,24(9):1845-1854.

[10] 趙慧玲,葉志明.拉壓不同模量彈性問題[J].上海大學學報(自然科學版)，2014,20(5):550-558.(ZHAO Hui-ling,YE Zhi-ming.Elastic Bodies with different modules in tension and compression[J].JournalofShanghaiUniversity(NaturalScience)，2014,20(5):550-558.(in Chinese))

[11] He X T,Chen S L,Sun J Y.Applying the equivalent section method to solve beam subjected to lateral force and bending-compression column with different moduli [J].InternationalJournalofMechanicalSciences,2007,49(7):919-924.

[12] 吳 曉，楊立軍.拉壓彈性模量不同曲梁的彈性理論解[J].工程力學,2013,30(1):76-80.(WU Xiao,YANG Li-jun.The elastic theory solution for curved beam with difference elastic modulus in tension and compression[J].EngineeringMechanics,2013,30(1):76-80.(in Chinese))

[13] 張良飛,姚文娟.拉壓不同模量矩形板雙向彎曲問題[J].上海大學學報(自然科學版)，2017,23(1):128-137.(ZHANG Liang-fei,YAO Wen-juan.Biaxial bending of rectangular plates with different modulus[J].JournalofShanghaiUniversity(NaturalScience)，2017,23(1):128-137.(in Chinese))

[14] 楊海天,鄔瑞鋒,楊克儉,等.初應力法解拉壓雙彈性模量問題[J].大連理工大學學報,1992,32(1):35-39.(YANG Hai-tian,WU Rui-feng,YANG Ke -jian,et al.Solution to problem of dual extension-compre -ssion elastic modulus with initial stress method[J].JournalofDalianUniversityofTechnology,1992,32(1):35-39.(in Chinese))

[15] 張國慶，楊海天.蟻群算法求解二維拉壓不同模量反問題[J].計算力學學報,2014,31(6):687-693.(ZHANG Guo -qing,YANG Hai-tian.Ant colony algorithm based numerical solution for inverse bimodular problems[J].ChineseJournalofComputationalMechanics,2014,31(6):687-693.(in Chinese))

[16] He X T,Zheng Z L,Sun J Y,et al.Convergence an-alysis of a finite element method based on different moduli in tension and compression[J].InternationalJournalofSolidsandStructures,2009,46(20):3734-3740.

[17] Zhang H W,Zhang L,Gao Q.An efficient computational method for mechanical analysis of bimodular structures based on parametric variational principle[J].Computers&Structures,2011,89(23)：2352-2360.

[18] 楊海天，朱應利.光滑函數(shù)法求解拉壓不同彈性模量問題[J].計算力學學報，2006，23(1)：19-23.(YANG Hai-tian,ZHU Ying-li.Solving elasticity problems with bi-modulus via a smoothing technique[J].ChineseJournalofComputationalMechanics，2006，23(1)：19-23.(in Chinese))

[19] Lee I W,Jung G H.An efficient algebraic method for the computation of natural frequency and mode shape sensitivities(part I):distinct natural frequencies[J].Computers&Structures,1997,62(3):429-435.

[20] Cai K.A simple approach to find optimal topology of a continuum with tension-only or compression-only material[J].StructuralandMultidisciplinaryOptimization,2011,43(6):827-835.