劉文武，陸念力，胡長(zhǎng)勝

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

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

劉文武，陸念力，胡長(zhǎng)勝

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

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

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

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

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

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

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

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

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

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

2 接觸單元?jiǎng)偠染仃?/h2>

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

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

式中:Ni為每個(gè)數(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單位矩陣。其中，

式中:ζ、η分別為單元長(zhǎng)度和寬度方向的物理坐標(biāo)。

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

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

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

因此，接觸單元每個(gè)節(jié)點(diǎn)的相對(duì)位移方程:

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

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

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

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

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

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

3 動(dòng)力接觸模型

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

式中:

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

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

4 試驗(yàn)驗(yàn)證

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

圖2 試驗(yàn)裝置和測(cè)試點(diǎn)布置圖Fig．2 Experimental device and site layout of test points

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

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

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

4．2 動(dòng)力接觸分析

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

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

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

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

圖5 測(cè)試點(diǎn)P1的接觸位移Fig．5 Contact displacement of point P1

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

圖7 測(cè)試點(diǎn)P3的接觸位移Fig．7 Contact displacement of point P3

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

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

5 結(jié)論

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

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

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

3）通過(guò)和現(xiàn)場(chǎng)試驗(yàn)結(jié)果比較，表明在單元數(shù)相同條件時(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ù)在時(shí)域匹配邊界元中的應(yīng)用研究［J］．哈爾濱工程大學(xué)學(xué)報(bào)，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］苗雨，李威，鄭俊杰，等．改進(jìn)的薄層單元在樁-土動(dòng)力相互作用中的應(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］齊良鋒，簡(jiǎn)浩，唐麗云．引入接觸單元模擬樁土共同作用［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é)報(bào)，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

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

2015-10-12．

2016-09-28．

國(guó)家國(guó)際科技合作專項(xiàng)項(xiàng)目（2013DFA71120，2015DFA70100）；中央高?；A(chǔ)研發(fā)基金項(xiàng)目（HIT．NSRIF．2013049）．

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

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

胡長(zhǎng)勝（1971-），男，副教授．

胡長(zhǎng)勝，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．