刘文武，陆念力，胡长胜

（哈尔滨工业大学机电工程学院，黑龙江哈尔滨150001）

等厚度三维单元在薄板-土动力接触中的应用

刘文武，陆念力，胡长胜

（哈尔滨工业大学机电工程学院，黑龙江哈尔滨150001）

针对薄板与土体的动力接触问题，根据土体小变形和基于薄板理论，对土体和薄板分别采用一阶基函数和零阶基函数进行建模，建立了一种基于一阶数值流形方法的等厚度三维接触单元。采用虚位移原理和刚塑性本构模型，推导了接触单元的刚度矩阵。通过引入阻尼单元和采用等价线性模型，建立了接触系统的二阶动力微分方程，并利用四阶Runge-Kutta显式方法编程求解得到了薄板与土体动力接触过程中的位移与应力值。相比于Goodman单元，在系统单元数相同条件下，当接触单元厚度为其宽度的0．06倍时，文中结果更符合试验值。通过和原位试验结果比较分析表明，文中建立的接触单元不仅建模方便，而且精度更高。

薄板-土相互作用；数值流形方法；接触单元；等厚度；原位试验；刚度矩阵；阻尼单元；本构模型

众所周知，薄板结构已广泛应用于力学和土木工程等领域，薄板-土相互作用分析在各种结构基础设计中也越来越重要，并成为近年来国内外探讨的热门问题［1-2］。薄板与土体之间接触应力和位移的分布在其动力相互作用过程中至关重要。因此，研究薄板-土动力相互作用中的接触关系具有非常重要的理论和工程价值。

接触问题的非线性特征，使得解析方法难于处理应用力学中的土-结构边值问题。有限元［3-5］和边界元［6-8］等数值分析方法由于通用性而在接触问题中被普遍采用，但高计算消耗是这些数值方法的共同缺点。R．E．Goodman等［9］首先提出了无厚度的接触单元对平面问题进行仿真。苗雨等［10］提出了用改进薄层单元来解决桩与土三维接触问题。齐良锋等［11］提出了在桩与土界面建立接触单元来模拟其相互间的共同作用，雷晓燕等［12］在其基础上进行了改进，通过将节点接触力作为基本未知量，并利用虚位移原理推导了等价接触单元刚度矩阵。为改善接触单元精度，Desai C S等［13］在Goodman单元基础上提出薄层单元，考虑了厚度对接触分析的影响。L．Gaul等［14］提出了一种片段之间的接触单元来模拟多自由度的复杂系统。上述基于有限元等方法所建立的接触单元，均无法在保持单元数不变的条件下增加接触分析的精度。K．Alexander等［15-16］指出阶数提高对描述接触问题能够带来更高的精度。系统单元数保持不变情况下，数值流形方法由于采用有限覆盖技术，只须提高基函数的阶数则能提高求解精度，且对于拥有不规则边界的实体进行建模更加方便［17-21］。

鉴于现有方法或建立的单元很难兼顾求解精度和计算耗费，本文提出基于一阶数值流形方法的等厚度三维接触单元。

1 接触单元结构

在薄板-土接触模型中建立长度和宽度分别为a和b的等厚度接触单元如图1所示，接触单元的厚度为h，有1～8共八个节点。其中，节点1～4为八节点土体等参单元的表面节点，节点5～8为薄板单元的表面节点，节点9～12分别为薄板单元的中面节点。板单元每个节点有3个转动自由度和3个移动自由度，土体等参单元每个节点具有3个移动自由度。

图1 接触模型和接触单元Fig．1 Contact model and contact element

图1中OXYZ为接触系统全局坐标系，O1ζηξ为接触单元和土体单元局部坐标系，O2xyz为薄板单元局部坐标系，薄板单元的厚度为2t。各个局部坐标系原点均位于相应单元的几何中心。

2 接触单元刚度矩阵

土体和薄板单元分别基于一阶和零阶数值流形方法进行建模，则土体单元每个节点位移可表示为

式中:Si采用数值流形方法中的一阶基函数，具体方程表示为

式中:Ni为每个数学覆盖的权函数，ξ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单位矩阵。其中，

式中:ζ、η分别为单元长度和宽度方向的物理坐标。

薄板单元每个节点的位移方程表示为

式中:Sj矩阵能将薄板单元转动自由度转化为移动自由度，具体如方程表示为

式中:Nj为接触单元每个数学覆盖的权函数，di=｛dij｝（i=9，10，11，12；j=1，2，3，4，5，6）。

因此，接触单元每个节点的相对位移方程:

式中:Se=［S1S2S3S4-S5-S6-S7-S8］为接触单元的相对应变矩阵。

由于接触单元厚度较小，故假设沿厚度方向的应变分量没有变化，根据虚位移原理:

则等厚度接触单元的刚度矩阵方程表示为

式中:h和D分别为接触单元厚度和本构关系矩阵。

式中:W=（1+kn）（1+kt）-μ2，kn、kt分别为法向和切向刚度，取值可参考文献［22］。E、μ和Cs分别为接触材料的弹性模量、泊松比和引入参数。

则整个接触系统的刚度矩阵方程表示为

3 动力接触模型

采用等价线性模型作为土体动力本构关系模型，则系统的动力方程为

式中:

其中，K、C和M分别为接触单元的刚度、阻尼和质量矩阵，F（t）和U分别为动力接触系统的力向量和位移向量，Ds和ω分别为土体的阻尼比和固有频率。

将接触单元阻尼矩阵组装到系统阻尼矩阵当中，假设在短时间间隔内系统的平均加速度为定值，采用四阶Runge-Kutta显示方法求解动力方程并编写动力接触关系程序则可求得动力作用下，薄板与土体之间的接触应力和接触位移。

4 试验验证

在哈尔滨市王岗地区对薄板-土体动力接触关系进行现场试验，试验装置如图2（a）所示。液压系统压力为12 MPa，试验装置的力学分析过程见文献［23］。孔壁土体中测试点的布置如图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列出了土体的其他参数。

图2 试验装置和测试点布置图Fig．2 Experimental device and site layout of test points

表1 试验区土体参数Table 1 Parameters of the soil of the test area

4．1 数值分析模型

六面体作为系统的数学覆盖，数值模型建模过程如图3（a）所示。由数值流形方法建模原理，得到系统的物理覆盖如图3（b）所示。最终系统的数值分析模型如图4所示，其中虚线区域表示接触单元。考虑剪切波速和土体-薄板的相互作用频率较小，模型水平和竖直方向单元尺寸分别为50 mm和100 mm。接触区域单元尺寸为50×50 mm2。最终薄板结构由64个流形元组成，每个流形元有4个数学覆盖，土体由23 000个流形元组成，每个流形元有8个数学覆盖，薄板与土体之间共有64个接触单元。土体模型四周为粘弹性边界，避免产生过大反射波［24］，土体底部为固定约束。

4．2 动力接触分析

在激振器作用下，由文献［23］中的力学分析可知，F2（t）=230sin（0．1πt）。数值仿真过程中，讨论厚度h分别等于0．01、0．02…0．1倍接触单元宽度b时的本文解，并与试验值和Goodman单元解相比较，Goodman单元法向和切向刚度按文献［9］均取为20 MPa／m。最终P1、P3的接触位移和P2、P4的接触应力如图5～8所示，其中循环时间为60 s，采样频率为1 Hz。

图3 数值模型的建模过程Fig．3 Modeling process of the numerical model

图4 数值模型轴向剖视图Fig．4 Axial cross-sectional view of the numerical model

由图5和图6分析可知，在不同接触单元的厚度下，本文解的结果相差较大。当接触单元厚度h为其宽度b的0．06倍时，本文解相比于Goodman单元解更加接近试验值。在整个循环时间内P1点接触位移的本文解相对于试验值的最小误差为5．57％，P2点接触应力的本文解相对于试验值的最小误差为5．68％。结果表明，在系统单元数相同条件下，增加覆盖函数的阶次能够获得更高的求解精度。

图5 测试点P1的接触位移Fig．5 Contact displacement of point P1

图6 测试点P2的接触应力Fig．6 Contact stress of point P2

图7 测试点P3的接触位移Fig．7 Contact displacement of point P3

图8 测试点P4的接触应力Fig．8 Contact stress of point P4

由图7和图8分析可知，不同接触单元的厚度下，本文解的值也存在较大差异。同样当接触单元厚度为其宽度的0．06倍时，本文解相比于Goodman单元解更加接近试验值。试验过程中激振器将激振力传递给连杆，由于测试点P3和P4位于薄板轴向中间线上，在轴向连杆组的共同作用下，P3点的接触位移和P4点的接触应力都将产生拍振。在整个循环时间内P3点接触位移的本文解相对于试验值的最小误差为6．58％，P4点接触应力的本文解相对于试验值的最小误差为6．78％。对比结果也表明，在系统单元数相同条件下，增加覆盖函数的阶次能够获得更高的求解精度。

5 结论

本文提出一种基于一阶数值流形方法的等厚度三维接触单元来模拟薄板与土体之间的动力接触关系。采用刚塑性本构关系模型和根据虚位移原理，推导了接触单元的刚度矩阵。引入Rayleigh阻尼反应薄板-土体系统在动力接触过程中的能量耗散。通过现场试验验证了本文所建立单元的有效性和精确性。具体结论如下:

1）采用数值流形方法通过规则的六面体对不规则的接触系统进行建模，简化了建模过程。

2）针对薄板-土体接触模型，接触单元厚度为其宽度的0．06倍时，仿真效果更符合试验值，并提出了厚度公式，其正确性需后续进一步验证。

3）通过和现场试验结果比较，表明在单元数相同条件时所建立接触单元比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，FUENMAYOR 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函数在时域匹配边界元中的应用研究［J］．哈尔滨工程大学学报，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］．岩土力学，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］．岩土力学，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，FU 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］．岩土工程学报，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

刘文武，陆念力，胡长胜．等厚度三维单元在薄板-土动力接触中的应用［J］．哈尔滨工程大学学报，2016，37（12）:1698-1703．

2015-10-12．

2016-09-28．

国家国际科技合作专项项目（2013DFA71120，2015DFA70100）；中央高校基础研发基金项目（HIT．NSRIF．2013049）．

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

陆念力（1955-），男，教授，博士生导师；

胡长胜（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．