動(dòng)態(tài)松弛方法中虛密度與時(shí)間步長(zhǎng)的選擇

陳新民，劉傳奇，薛世峰

(1.中國(guó)石油大學(xué)儲(chǔ)運(yùn)與建筑工程學(xué)院，山東青島 266580;2.勝利油田濱南采油廠，山東濱州 256600)

動(dòng)態(tài)松弛法是將線(xiàn)性或非線(xiàn)性靜力問(wèn)題轉(zhuǎn)化為動(dòng)力問(wèn)題，通過(guò)設(shè)置阻尼消耗能量，當(dāng)位移收斂到平衡位置(即勢(shì)能最小位置)時(shí)，迭代結(jié)束的顯示求解方法［1］。Otter［2］在分析預(yù)應(yīng)力混凝土核容器外殼時(shí)首次提出該方法;Cassell等［3］引入虛質(zhì)量的概念;Underwood［4］對(duì)如何設(shè)定最有效的阻尼矩陣、質(zhì)量矩陣進(jìn)行了研究;David等［5］采用動(dòng)態(tài)松弛法進(jìn)行非線(xiàn)性分析，研究了超彈性結(jié)構(gòu);Michael［6］采用此方法進(jìn)行了受拉結(jié)構(gòu)的找形分析;Li等［7］應(yīng)用動(dòng)態(tài)松弛法求解可變形離散元中的節(jié)點(diǎn)運(yùn)動(dòng)方程;Rezaiee等［8］對(duì)虛質(zhì)量與阻尼的具體算法進(jìn)行了系統(tǒng)論述;李士軍等［9］對(duì)Rayleigh阻尼參數(shù)取值進(jìn)行了分析。當(dāng)求解包含應(yīng)力集中或其他局部網(wǎng)格加密的問(wèn)題時(shí)，可通過(guò)設(shè)置虛密度，提高收斂速度［10］。對(duì)于虛密度影響收斂時(shí)步與計(jì)算結(jié)果的問(wèn)題，相關(guān)文獻(xiàn)研究較少。劉波等［11］從應(yīng)力波傳播角度討論了不同單元如何設(shè)定密度與步長(zhǎng)，但未能給出虛密度影響下時(shí)步設(shè)定的具體表達(dá)式。筆者以常應(yīng)變?nèi)切螁卧獮槔?，探討虛密度與時(shí)步的設(shè)定問(wèn)題。

1 密度對(duì)收斂過(guò)程與結(jié)果的影響

從節(jié)點(diǎn)運(yùn)動(dòng)角度，可變形離散單元保持收斂的必要條件是單元Jacobi行列式恒大于零。如圖1所示，在一個(gè)時(shí)步內(nèi)，若節(jié)點(diǎn)k從原位置移動(dòng)到k1位置，網(wǎng)格未出現(xiàn)畸變，保證在平衡位置附近振動(dòng)，若移動(dòng)到k2(Jacobi行列式為零)，k3(Jacobi行列式小于零)位置，則網(wǎng)格畸變，引起發(fā)散。

圖1 單元節(jié)點(diǎn)位置示意圖Fig.1 Schematic diagram of nodes position

遍歷所有網(wǎng)格，初始狀態(tài)頂點(diǎn)到對(duì)邊最小距離記為Hmin，則計(jì)算收斂的必要條件是:Δt時(shí)間內(nèi)，內(nèi)含集中質(zhì)量的節(jié)點(diǎn)位移小于Hmin。有

式中，F(xiàn)ext為k節(jié)點(diǎn)所受外力，N;mk為k節(jié)點(diǎn)集中質(zhì)量(在三角形常應(yīng)變單元中，mk=ρA/3)，kg;A為單元面積，m2;ρ為材料密度，kg/m3。

將mk帶入式(1)中得

由式(2)可以明顯看出，若單元密度增加，則收斂時(shí)間步長(zhǎng)可相應(yīng)增加。

從信息傳遞角度，迭代收斂的必要條件是當(dāng)前計(jì)算速度必須大于信息傳播的最大速度［11］，由此推導(dǎo)的最大收斂時(shí)間步長(zhǎng)亦正比于密度的0.5次方。

針對(duì)密度變化影響最終平衡位置改變的問(wèn)題，考慮單自由度含阻尼強(qiáng)迫振動(dòng)，其運(yùn)動(dòng)方程為

其中

式中，ε為阻尼比;p為固有頻率，rad/s;m為節(jié)點(diǎn)質(zhì)量，kg;c為阻尼系數(shù)，N·s/m;k為彈簧剛度，N/m;F為所受外力，N。

根據(jù)杜哈梅積分公式［12］，初始條件為零的情況下，方程(3)的全解為

利用分部積分得到

若按式(2)設(shè)定步長(zhǎng)，將t=nΔt代入式(5)，則有

由式(6)可以看出，x(n)與密度無(wú)關(guān)，即按最大收斂時(shí)步進(jìn)行迭代，收斂速度相同。因而，加快收斂的核心是設(shè)定虛密度和步長(zhǎng)，使不同單元擁有相同的收斂時(shí)步，并且步長(zhǎng)盡可能接近最大收斂時(shí)步。

2 虛密度與步長(zhǎng)設(shè)定方法

單元節(jié)點(diǎn)運(yùn)動(dòng)受彈簧作用力影響，偏離平衡位置位移越大，所受彈簧力越大。彈簧力影響節(jié)點(diǎn)加速度，導(dǎo)致節(jié)點(diǎn)速度發(fā)生變化，單個(gè)時(shí)步的位移受到影響，反過(guò)來(lái)，影響下一時(shí)步的彈簧力。各因素影響關(guān)系如圖2所示。

圖2 迭代參數(shù)相互影響示意圖Fig.2 Schematic diagram of iterative parameters mutual influence

圖中，aij、vij、sij分別表示第 i步 j節(jié)點(diǎn)的加速度、速度、位移。假設(shè)j節(jié)點(diǎn)的最大收斂時(shí)間步長(zhǎng)為Δtj，則有

采用三次迭代后，節(jié)點(diǎn)速度至少衰減為初速度一半的標(biāo)準(zhǔn)，確定時(shí)間步長(zhǎng)，即:v3j≤0.5v1j，帶入各自表達(dá)式，取等號(hào)進(jìn)行方程運(yùn)算，則有

此方程為一元五次方程，共有5個(gè)不同的解，分別為

整個(gè)過(guò)程中未考慮阻尼影響，計(jì)算最大收斂時(shí)步需添加阻尼影響系數(shù)c，則j節(jié)點(diǎn)最大收斂時(shí)步為

式中，ρj為所設(shè)虛密度，kg/m3;Aj為單元面積，m2;k為彈簧剛度，N/m;c0為比例系數(shù)。

系統(tǒng)的收斂時(shí)步為ΔT=min{Δtj}，在實(shí)際應(yīng)用中，首先計(jì)算最小單元面積，記為Amin，假定其密度，將單元j的密度設(shè)定為

此時(shí)，各單元收斂時(shí)步相同，再按照

設(shè)定步長(zhǎng)，其中c0通過(guò)反復(fù)試算，取為0.45。

3 編程驗(yàn)證

基于可變形離散單元法進(jìn)行編程驗(yàn)證。以平面受拉薄板為例，驗(yàn)證理論推導(dǎo)結(jié)論以及虛密度設(shè)定方案的有效性。計(jì)算模型如圖3所示，由于薄板受集中力作用，采用不均勻網(wǎng)格(圖3(b))。原來(lái)連續(xù)的節(jié)點(diǎn)分裂為夾有彈簧的兩個(gè)角點(diǎn)，圖中每個(gè)網(wǎng)格線(xiàn)間均有一對(duì)彈簧，方向分別為界面法向與切向。

為減弱由彈簧變形引起的位移不連續(xù)，將彈簧剛度設(shè)為單元彈性模量的50倍［7］，阻尼采用局部阻尼系數(shù)0.8［13］，彈性模量為20 GPa，泊松比為 0.2，密度為1.8×103kg/m3，彈簧剛度為1 000 GPa，阻尼系數(shù)為0.8。

圖3 計(jì)算模型示意圖Fig.3 Schematic diagram of calculation model

采用C++語(yǔ)言編寫(xiě)可變形離散元程序。計(jì)算過(guò)程為:單元與單元之間通過(guò)彈簧相連，根據(jù)單元節(jié)點(diǎn)的相對(duì)位移計(jì)算彈簧作用力，采用高斯定理和Wilkins的大變形［14］有限差分法計(jì)算單元變形力，采用局部無(wú)黏性阻尼計(jì)算阻尼力，按照動(dòng)態(tài)松弛法不平衡力—加速度—速度—位移—不平衡力的順序進(jìn)行迭代，當(dāng)不平衡力數(shù)值小于預(yù)定值時(shí)，迭代完成。需要指出，速度與位移的求解是按照中心差分格式進(jìn)行的。

圖4為可變形離散元法和有限元軟件ANSYS的計(jì)算結(jié)果，Y向最大位移分別為4.6×10－7m和4.7×10－7m，兩種結(jié)果的相對(duì)誤差為2.1%。通過(guò)改變單元密度進(jìn)行計(jì)算，結(jié)果并未發(fā)生變化。

圖4 可變形離散元法和ANSYS軟件的計(jì)算結(jié)果Fig.4 Computed result of deformable discrete element method and ANSYS software

選取頂部左節(jié)點(diǎn)的內(nèi)力作為參照，分別對(duì)各單元具有相同密度、不同密度的情況進(jìn)行討論，以驗(yàn)證推導(dǎo)結(jié)論及虛密度、步長(zhǎng)設(shè)定方案的有效性。

3.1 不同單元具有相同密度

對(duì)比3種相同密度下的最大收斂時(shí)步與公式(11)的估計(jì)時(shí)步，如表1所示。其中最大收斂時(shí)步通過(guò)反復(fù)試算求得。

表1 收斂時(shí)步與估計(jì)時(shí)步對(duì)照Table 1 Comparison of convergence and estimated time step

(1)密度相同(1.8×103kg/m3)、時(shí)步不同時(shí)，收斂曲線(xiàn)如圖5(a)所示。由圖5(a)可以看出，時(shí)步增加，收斂速度加快，這是由于密度不變，最大收斂時(shí)步不變，若時(shí)步增加，則時(shí)步比增加，增大了收斂速度。

圖5 相同密度下不同時(shí)步對(duì)照收斂曲線(xiàn)Fig.5 Comparison of convergence curve in different time step under same density

(2)時(shí)步相同(2.0×10－6s)，密度不同時(shí)，收斂曲線(xiàn)如圖5(b)所示。由圖5(b)可以看出，密度增加，收斂速度減緩。這是由于密度增加，則最大收斂時(shí)步增加，而時(shí)步恒定，則時(shí)步比減小，因而放緩了收斂速度。

(3)時(shí)步比均為1，密度不同時(shí)，收斂曲線(xiàn)如圖5(c)所示。由圖5(c)可以看出，時(shí)步比相同時(shí)，迭代速度相差不大，驗(yàn)證了按最大收斂時(shí)步進(jìn)行迭代，收斂速度相同的結(jié)論。

3.2 不同單元具有不同密度

虛密度加速收斂的核心在于不同單元設(shè)定不同的密度，按照公式(10)進(jìn)行密度設(shè)定，公式(11)進(jìn)行時(shí)步設(shè)定，對(duì)比真實(shí)密度下按照最大收斂時(shí)步進(jìn)行迭代的收斂曲線(xiàn)，如圖6所示。由圖6可以看出，按照本文所述方法進(jìn)行密度設(shè)定，時(shí)步設(shè)定，加速收斂效果十分明顯，驗(yàn)證了此方法的有效性。

圖6 真實(shí)密度與虛密度對(duì)照收斂曲線(xiàn)Fig.6 Comparison of convergence curve in actual density and fictitious density

4 結(jié)論

(1)虛密度增大，可增大收斂時(shí)步。

(2)收斂速度與時(shí)步比有關(guān)，時(shí)步比越大，收斂速度越快。

(3)虛密度加速收斂的核心在于不同單元設(shè)置不同密度，使得不同單元擁有相同的最大收斂時(shí)步，按照所提出的密度、時(shí)步設(shè)定方案進(jìn)行參數(shù)設(shè)定，加速收斂效果很明顯。

［1］ TONG P.An adaptive dynamic relaxation method for static problems［J］.Computational Mechanics，1986(1):127-140.

［2］ OTTER J R H.Computations for prestressed concrete reactor pressure vessels using dynamic relaxation［J］.Nuclear Structural Engineering，1965，1(1):61-75.

［3］ CASSELL A C，KINSEY P J，SEFTON D J.Cylindrical shell analysis by dynamic relaxation［J］.ICE Proceedings，1968，39(1):75-84.

［4］ UNDERWOOD P.Dynamic relaxation［R］.Lockheed Report LMSC-D678265，1979.

［5］ DAVID R Oakley，NORMAN F Knight Jr.Adaptive dynamic relaxation algorithm for non-linear hyperelastic structures(Ⅰ):formulation［J］.Computer Methods in Applied Mechanics and Engineering，1995，126(1/2):67-89.

［6］ MICHAEL R Barnes.Form finding and analysis of tension structures by dynamic relaxation［J］.International Journal of Space Structures，1999，14(2):89-104.

［7］ LI Shi-hai，LIU Xiao-yu，LIU Tian-ping，et al.Continuum-based discrete element method and its applications［C］.Beijing:UK-China Summer School/International Symposium on DEM，2008:24-30.

［8］ REZAIEE-PAJAND M，ALAMATIAN J.The dynamic relaxation method using new formulation for fictitious mass and damping［J］.Structural Engineering and Mechanics，2010，34(1):109-133.

［9］ 李士軍，馬大為，朱孫科.動(dòng)態(tài)松弛方法中Rayleigh阻尼參數(shù)取值分析［J］.計(jì)算力學(xué)學(xué)報(bào)，2010，27(1):169-172.LI Shi-jun，MA Da-wei，ZHU Sun-ke.Analysis of rayleigh danp parameters in dynamic relaxation method［J］.Chinese Journal of Computational Mechanics，2010，27(1):169-172.

［10］ 王泳嘉，邢紀(jì)波.離散單元法及其在巖土力學(xué)中的應(yīng)用［M］.沈陽(yáng):東北工學(xué)院出版社，1991:46-47.

［11］ 劉波，韓彥輝.FLAC原理、實(shí)例與應(yīng)用指南［M］.北京:人民交通出版社，2005:11-12.

［12］ 仇偉德.機(jī)械振動(dòng)［M］.東營(yíng):石油大學(xué)出版社，2001:45-46.

［13］ PAPADRAKAKIS M.A method for the automatic evaluation of the dynamic relaxation parameters［J］.Computer Methods in Applied Mechanics and Engineering，1981，25(1):35-48.

［14］ LANRU Jing，SSTEPHANSSON O.Fundamentals of discrete element methods for rock engineering:theory and applications［M］.Amsterdam:Elsevier，2007.