等厚度三維單元在薄板-土動力接觸中的應(yīng)用

劉文武，陸念力，胡長勝

（哈爾濱工業(yè)大學(xué)機電工程學(xué)院，黑龍江哈爾濱150001）

等厚度三維單元在薄板-土動力接觸中的應(yīng)用

劉文武，陸念力，胡長勝

（哈爾濱工業(yè)大學(xué)機電工程學(xué)院，黑龍江哈爾濱150001）

針對薄板與土體的動力接觸問題，根據(jù)土體小變形和基于薄板理論，對土體和薄板分別采用一階基函數(shù)和零階基函數(shù)進行建模，建立了一種基于一階數(shù)值流形方法的等厚度三維接觸單元。采用虛位移原理和剛塑性本構(gòu)模型，推導(dǎo)了接觸單元的剛度矩陣。通過引入阻尼單元和采用等價線性模型，建立了接觸系統(tǒng)的二階動力微分方程，并利用四階Runge-Kutta顯式方法編程求解得到了薄板與土體動力接觸過程中的位移與應(yīng)力值。相比于Goodman單元，在系統(tǒng)單元數(shù)相同條件下，當接觸單元厚度為其寬度的0．06倍時，文中結(jié)果更符合試驗值。通過和原位試驗結(jié)果比較分析表明，文中建立的接觸單元不僅建模方便，而且精度更高。

薄板-土相互作用；數(shù)值流形方法；接觸單元；等厚度；原位試驗；剛度矩陣；阻尼單元；本構(gòu)模型

眾所周知，薄板結(jié)構(gòu)已廣泛應(yīng)用于力學(xué)和土木工程等領(lǐng)域，薄板-土相互作用分析在各種結(jié)構(gòu)基礎(chǔ)設(shè)計中也越來越重要，并成為近年來國內(nèi)外探討的熱門問題［1-2］。薄板與土體之間接觸應(yīng)力和位移的分布在其動力相互作用過程中至關(guān)重要。因此，研究薄板-土動力相互作用中的接觸關(guān)系具有非常重要的理論和工程價值。

接觸問題的非線性特征，使得解析方法難于處理應(yīng)用力學(xué)中的土-結(jié)構(gòu)邊值問題。有限元［3-5］和邊界元［6-8］等數(shù)值分析方法由于通用性而在接觸問題中被普遍采用，但高計算消耗是這些數(shù)值方法的共同缺點。R．E．Goodman等［9］首先提出了無厚度的接觸單元對平面問題進行仿真。苗雨等［10］提出了用改進薄層單元來解決樁與土三維接觸問題。齊良鋒等［11］提出了在樁與土界面建立接觸單元來模擬其相互間的共同作用，雷曉燕等［12］在其基礎(chǔ)上進行了改進，通過將節(jié)點接觸力作為基本未知量，并利用虛位移原理推導(dǎo)了等價接觸單元剛度矩陣。為改善接觸單元精度，Desai C S等［13］在Goodman單元基礎(chǔ)上提出薄層單元，考慮了厚度對接觸分析的影響。L．Gaul等［14］提出了一種片段之間的接觸單元來模擬多自由度的復(fù)雜系統(tǒng)。上述基于有限元等方法所建立的接觸單元，均無法在保持單元數(shù)不變的條件下增加接觸分析的精度。K．Alexander等［15-16］指出階數(shù)提高對描述接觸問題能夠帶來更高的精度。系統(tǒng)單元數(shù)保持不變情況下，數(shù)值流形方法由于采用有限覆蓋技術(shù)，只須提高基函數(shù)的階數(shù)則能提高求解精度，且對于擁有不規(guī)則邊界的實體進行建模更加方便［17-21］。

鑒于現(xiàn)有方法或建立的單元很難兼顧求解精度和計算耗費，本文提出基于一階數(shù)值流形方法的等厚度三維接觸單元。

1 接觸單元結(jié)構(gòu)

在薄板-土接觸模型中建立長度和寬度分別為a和b的等厚度接觸單元如圖1所示，接觸單元的厚度為h，有1～8共八個節(jié)點。其中，節(jié)點1～4為八節(jié)點土體等參單元的表面節(jié)點，節(jié)點5～8為薄板單元的表面節(jié)點，節(jié)點9～12分別為薄板單元的中面節(jié)點。板單元每個節(jié)點有3個轉(zhuǎn)動自由度和3個移動自由度，土體等參單元每個節(jié)點具有3個移動自由度。

圖1 接觸模型和接觸單元Fig．1 Contact model and contact element

圖1中OXYZ為接觸系統(tǒng)全局坐標系，O1ζηξ為接觸單元和土體單元局部坐標系，O2xyz為薄板單元局部坐標系，薄板單元的厚度為2t。各個局部坐標系原點均位于相應(yīng)單元的幾何中心。

2 接觸單元剛度矩陣

土體和薄板單元分別基于一階和零階數(shù)值流形方法進行建模，則土體單元每個節(jié)點位移可表示為

式中:Si采用數(shù)值流形方法中的一階基函數(shù)，具體方程表示為

式中:Ni為每個數(shù)學(xué)覆蓋的權(quán)函數(shù)，ξi、ηi和ζi為中間自由度。di=｛dij｝（i=1，2，3，4；j=1，2，3，4，5，6，7，8，9，10，11，12），I為3×3單位矩陣。其中，

式中:ζ、η分別為單元長度和寬度方向的物理坐標。

薄板單元每個節(jié)點的位移方程表示為

式中:Sj矩陣能將薄板單元轉(zhuǎn)動自由度轉(zhuǎn)化為移動自由度，具體如方程表示為

式中:Nj為接觸單元每個數(shù)學(xué)覆蓋的權(quán)函數(shù)，di=｛dij｝（i=9，10，11，12；j=1，2，3，4，5，6）。

因此，接觸單元每個節(jié)點的相對位移方程:

式中:Se=［S1S2S3S4-S5-S6-S7-S8］為接觸單元的相對應(yīng)變矩陣。

由于接觸單元厚度較小，故假設(shè)沿厚度方向的應(yīng)變分量沒有變化，根據(jù)虛位移原理:

則等厚度接觸單元的剛度矩陣方程表示為

式中:h和D分別為接觸單元厚度和本構(gòu)關(guān)系矩陣。

式中:W=（1+kn）（1+kt）-μ2，kn、kt分別為法向和切向剛度，取值可參考文獻［22］。E、μ和Cs分別為接觸材料的彈性模量、泊松比和引入?yún)?shù)。

則整個接觸系統(tǒng)的剛度矩陣方程表示為

3 動力接觸模型

采用等價線性模型作為土體動力本構(gòu)關(guān)系模型，則系統(tǒng)的動力方程為

式中:

其中，K、C和M分別為接觸單元的剛度、阻尼和質(zhì)量矩陣，F(xiàn)（t）和U分別為動力接觸系統(tǒng)的力向量和位移向量，Ds和ω分別為土體的阻尼比和固有頻率。

將接觸單元阻尼矩陣組裝到系統(tǒng)阻尼矩陣當中，假設(shè)在短時間間隔內(nèi)系統(tǒng)的平均加速度為定值，采用四階Runge-Kutta顯示方法求解動力方程并編寫動力接觸關(guān)系程序則可求得動力作用下，薄板與土體之間的接觸應(yīng)力和接觸位移。

4 試驗驗證

在哈爾濱市王崗地區(qū)對薄板-土體動力接觸關(guān)系進行現(xiàn)場試驗，試驗裝置如圖2（a）所示。液壓系統(tǒng)壓力為12 MPa，試驗裝置的力學(xué)分析過程見文獻［23］?？妆谕馏w中測試點的布置如圖2（b）所示。薄板的長度為hp=0．7 m，外半徑R=0．75 m，中心角φ=50°，厚度為tp=0．01 m，彈性模量E和泊松比νp分別為2 100 GPa和0．3。試驗中樁孔的半徑R= 1．35 m，深度h=2．0 m，土體分析模型的半徑R0= 6．75 m，深度為H=10．0 m，表1列出了土體的其他參數(shù)。

圖2 試驗裝置和測試點布置圖Fig．2 Experimental device and site layout of test points

表1 試驗區(qū)土體參數(shù)Table 1 Parameters of the soil of the test area

4．1 數(shù)值分析模型

六面體作為系統(tǒng)的數(shù)學(xué)覆蓋，數(shù)值模型建模過程如圖3（a）所示。由數(shù)值流形方法建模原理，得到系統(tǒng)的物理覆蓋如圖3（b）所示。最終系統(tǒng)的數(shù)值分析模型如圖4所示，其中虛線區(qū)域表示接觸單元?？紤]剪切波速和土體-薄板的相互作用頻率較小，模型水平和豎直方向單元尺寸分別為50 mm和100 mm。接觸區(qū)域單元尺寸為50×50 mm2。最終薄板結(jié)構(gòu)由64個流形元組成，每個流形元有4個數(shù)學(xué)覆蓋，土體由23 000個流形元組成，每個流形元有8個數(shù)學(xué)覆蓋，薄板與土體之間共有64個接觸單元。土體模型四周為粘彈性邊界，避免產(chǎn)生過大反射波［24］，土體底部為固定約束。

4．2 動力接觸分析

在激振器作用下，由文獻［23］中的力學(xué)分析可知，F(xiàn)2（t）=230sin（0．1πt）。數(shù)值仿真過程中，討論厚度h分別等于0．01、0．02…0．1倍接觸單元寬度b時的本文解，并與試驗值和Goodman單元解相比較，Goodman單元法向和切向剛度按文獻［9］均取為20 MPa／m。最終P1、P3的接觸位移和P2、P4的接觸應(yīng)力如圖5～8所示，其中循環(huán)時間為60 s，采樣頻率為1 Hz。

圖3 數(shù)值模型的建模過程Fig．3 Modeling process of the numerical model

圖4 數(shù)值模型軸向剖視圖Fig．4 Axial cross-sectional view of the numerical model

由圖5和圖6分析可知，在不同接觸單元的厚度下，本文解的結(jié)果相差較大。當接觸單元厚度h為其寬度b的0．06倍時，本文解相比于Goodman單元解更加接近試驗值。在整個循環(huán)時間內(nèi)P1點接觸位移的本文解相對于試驗值的最小誤差為5．57％，P2點接觸應(yīng)力的本文解相對于試驗值的最小誤差為5．68％。結(jié)果表明，在系統(tǒng)單元數(shù)相同條件下，增加覆蓋函數(shù)的階次能夠獲得更高的求解精度。

圖5 測試點P1的接觸位移Fig．5 Contact displacement of point P1

圖6 測試點P2的接觸應(yīng)力Fig．6 Contact stress of point P2

圖7 測試點P3的接觸位移Fig．7 Contact displacement of point P3

圖8 測試點P4的接觸應(yīng)力Fig．8 Contact stress of point P4

由圖7和圖8分析可知，不同接觸單元的厚度下，本文解的值也存在較大差異。同樣當接觸單元厚度為其寬度的0．06倍時，本文解相比于Goodman單元解更加接近試驗值。試驗過程中激振器將激振力傳遞給連桿，由于測試點P3和P4位于薄板軸向中間線上，在軸向連桿組的共同作用下，P3點的接觸位移和P4點的接觸應(yīng)力都將產(chǎn)生拍振。在整個循環(huán)時間內(nèi)P3點接觸位移的本文解相對于試驗值的最小誤差為6．58％，P4點接觸應(yīng)力的本文解相對于試驗值的最小誤差為6．78％。對比結(jié)果也表明，在系統(tǒng)單元數(shù)相同條件下，增加覆蓋函數(shù)的階次能夠獲得更高的求解精度。

5 結(jié)論

本文提出一種基于一階數(shù)值流形方法的等厚度三維接觸單元來模擬薄板與土體之間的動力接觸關(guān)系。采用剛塑性本構(gòu)關(guān)系模型和根據(jù)虛位移原理，推導(dǎo)了接觸單元的剛度矩陣。引入Rayleigh阻尼反應(yīng)薄板-土體系統(tǒng)在動力接觸過程中的能量耗散。通過現(xiàn)場試驗驗證了本文所建立單元的有效性和精確性。具體結(jié)論如下:

1）采用數(shù)值流形方法通過規(guī)則的六面體對不規(guī)則的接觸系統(tǒng)進行建模，簡化了建模過程。

2）針對薄板-土體接觸模型，接觸單元厚度為其寬度的0．06倍時，仿真效果更符合試驗值，并提出了厚度公式，其正確性需后續(xù)進一步驗證。

3）通過和現(xiàn)場試驗結(jié)果比較，表明在單元數(shù)相同條件時所建立接觸單元比Goodman單元精度更高。

［1］MOONEY M，BOUTON C．Vibratory plate loading of compacted and instrumented field soil beds［J］．Geotechnical testing journal，2005，28（3）:221-230．

［2］AUERSCH L．Response to harmonic wave excitation of finite or infinite elastic plates on a homogeneous or layered halfspace［J］．Computers and geotechnics，2013，51:50-59．

［3］MAO Jianqiang．A finite element approach to solve contact problems in geotechnical engineering［J］．International journal for numerical and analytical methods in geomechanics，2005，29（5）:525-550．

［4］GINER E，SUKUMAR N，F(xiàn)UENMAYOR F J，et al．Singularity enrichment for complete sliding contact using the partition of unity finite element method［J］．International journal for numerical methods in engineering，2008，76（9）:1402-1418．

［5］SEGURA J M，CAROL I．On zero-thickness interface elements for diffusion problems［J］．International journal for numerical and analytical methods in geomechanics，2004，28（9）:947-962．

［6］GUN H．Boundary element analysis of 3-D elasto-plastic con-tact problems with friction［J］．Computers＆structures，2004，82（7／8）:555-566．

［7］PADRóN L A，MYLONAKIS G，BESKOS D E．Simple superposition approach for dynamic analysis of piled embedded footings［J］．International journal for numerical and analytical methods in geomechanics，2012，36（12）:1523-1534．

［8］李志富，任慧龍，石玉云，等．Laguerre函數(shù)在時域匹配邊界元中的應(yīng)用研究［J］．哈爾濱工程大學(xué)學(xué)報，2016，37（5）:629-633．LI Zhifu，REN Huilong，SHI Yuyun，et al．Research on the application of Laguerre function in the time domain Rankine-Green panel method［J］．Journal of Harbin Engineering University，2016，37（5）:629-633．

［9］GOODMAN R E，TAYLOR R L，BREKKE T L．A model for the mechanics of jointed rock［J］．Journal of the soil mechanics and foundations division，1968，94（3）:637-660．

［10］苗雨，李威，鄭俊杰，等．改進的薄層單元在樁-土動力相互作用中的應(yīng)用［J］．巖土力學(xué)，2015，36（11）: 3223-3228．MIAO Yu，LI Wei，ZHENG Junjie，et al．Application of modified thin layer element to the analysis of dynamic pilesoil interaction［J］．Rock and soil mechanics，2015，36（11）:3223-3228．

［11］齊良鋒，簡浩，唐麗云．引入接觸單元模擬樁土共同作用［J］．巖土力學(xué)，2005，26（1）:127-130．QI Liangfeng，JIAN Hao，TANG Liyun．The computing model of adopting contact element to simulate the pile-soil’s reciprocity［J］．Rock and soil mechanics，2005，26（1）: 127-130．

［12］LEI Xiaoyan．Contact friction analysis with a simple interface element［J］．Computer methods in applied mechanics and engineering，2001，190（15-17）:1955-1965．

［13］DESAI C S，ZAMAN M M，LIGHTNER J G，et al．Thinlayer element for interfaces and joints［J］．International journal for numerical and analytical methods in geomechanics，1984，8（1）:19-43．

［14］MAYER M H，GAUL L．Segment-to-segment contact elements for modelling joint interfaces in finite element analysis［J］．Mechanical systems and signal processing，2007，21（2）:724-734．

［15］KONYUKHOV A，SCHWEIZERHOF K．Incorporation of contact for high-order finite elements in covariant form［J］．Computer methods in applied mechanics and engineering，2009，198（13／14）:1213-1223．

［16］NGUYEN D T，RAUCHS G，PONTHOT J P．The impact of surface higher order differentiability in two-dimensional contact elements［J］．Journal of computational and applied mathematics，2013，246:195-205．

［17］PABISEK E．Self-learning FEM／NMM approach to identification of equivalent material models for plane stress problem［J］．Computer assisted mechanics and engineering sciences，2008，15（1）:67-78．

［18］ZHANG H W，ZHOU L．Numerical manifold method for dynamic nonlinear analysis of saturated porous media［J］．International journal for numerical and analytical methods in geomechanics，2006，30（9）:927-951．

［19］MA Guowei，AN Xinmei，HE Lei．The numerical manifold method:a review［J］．International journal of computational methods，2010，7（1）:1-32．

［20］QU X L，F(xiàn)U G Y，MA G W．An explicit time integration scheme of numerical manifold method［J］．Engineering analysis with boundary elements，2014，48:53-62．

［21］ZHENG Hong，XU Dongdong．New strategies for some issues of numerical manifold method in simulation of crack propagation［J］．International journal for numerical methods in engineering，2014，97（13）:986-1010．

［22］劉晶波，谷音，杜義欣．一致粘彈性人工邊界及粘彈性邊界單元［J］．巖土工程學(xué)報，2006，28（9）:1070-1075．LIU Jingbo，GU Yin，DU Yixin．Consistent viscous-spring artificial boundaries and viscous-spring boundary elements［J］．Chinese journal of geotechnical engineering，2006，28（9）:1070-1075．

［23］LIU Wenwu，HU Changsheng，LU Nianli．Bearing board’s force analysis of the rodless drilling set and its prototype test［J］．Applied mechanics and materials，2014，446-447:526-531．

［24］TAO Guilan，ZHANG Li．Dynamic response analysis of docking chamber structure considering viscoelastic boundary condition［J］．Applied mechanics and materials，2013，405-408:1939-1944．

Application of a 3D constant-thickness element in thin-plate-soil dynamic contact interaction

LIU Wenwu，LU Nianli，HU Changsheng
（School of Mechatronics Engineering，Harbin Institute of Technology，Harbin 150001，China）

To solve the dynamic contact problem between a thin plate and soil，a new 3D constant-thickness contact element based on the first-order numerical manifold method is proposed．In this method，the soil is modeled by a firstorder basic function and the thin plate is modeled by a zero-order basic function．Using the principle of virtual displacement and a rigid-plastic constitutive model，the stiffness matrix of the new contact element was deduced．By adding a damping component and equivalent linear model，the second-order differential dynamic equation of the contact system was built．Then，the displacement and stress values of the dynamic contact process were solved using the four-order Runge-Kutta explicit method．Compared with the Goodman results，the new contact element agrees better with the experimental results when the thickness of the contact element is 0．06 times its width．Comparison with the field test shows that the new contact element is not only convenient for modeling but also has higher accuracy．

thin-plate-soil interaction；numerical manifold method；contact element；constant thickness；field test；stiffness matrix；damping component；constitutive model

10．11990／jheu．201510020

http：／／www．cnki．net／kcms／detail／23．1390．u．20160928．0936．034．html

TU67

A

1006-7043（2016）12-1698-06

劉文武，陸念力，胡長勝．等厚度三維單元在薄板-土動力接觸中的應(yīng)用［J］．哈爾濱工程大學(xué)學(xué)報，2016，37（12）:1698-1703．

2015-10-12．

2016-09-28．

國家國際科技合作專項項目（2013DFA71120，2015DFA70100）；中央高校基礎(chǔ)研發(fā)基金項目（HIT．NSRIF．2013049）．

劉文武（1985-），男，博士研究生；

陸念力（1955-），男，教授，博士生導(dǎo)師；

胡長勝（1971-），男，副教授．

胡長勝，E-mail:hit_liuwenwu＠163．com．

LIU Wenwu，LU Nianli，HU Changsheng．Application of a 3D constant-thickness element in thin-plate-soil dynamic contact interaction［J］．Journal of Harbin Engineering University，2016，37（12）:1698-1703．