鱼雷锚上拔承载力的物质点法数值分析

程 钰，邱长林

(天津大学 建筑工程学院，天津 300072)

鱼雷锚是一种用于系泊深海设施的基础，形似鱼雷，其外表为锥形的圆柱形钢管，内用混凝土和废金属填充。安装时将鱼雷锚在一定水深处释放，其依靠自身重力下落获得的动能上拔海床，安装快速简便，是传统海上锚固件经济性的替代方案。2001年Medeiros C J[1]概述了鱼雷锚在巴西近海坎普斯湾(Campos Basin)进行鱼雷锚现场原位试验，水深范围200～1 000 m，并分析鱼雷锚的贯入深度与其上拔承载力，为鱼雷锚的工程施工规范和理论分析发展奠定了基础。国内目前一般海洋工程项目基础的安全承载力设计主要参考API规范[2]。但由于理论估算受参数影响大，结果存在不稳定性，近年来越来越多的学者采用模型试验和数值分析方法研究鱼雷锚上拔承载力。

Wenkai Wang等[3]进行了240组实验室试验，将11种不同形状的鱼雷锚从不同类型、不同埋深的粘性土体中垂直拉出，提出鱼雷锚不排水单调抗拔承载力估算公式。Hossain M S[4]等，通过200 g的离心机模型试验探究黏土和淤泥中鱼雷锚的单向抗拔承载力，试验结果表明，承载力与安装后固结时间、锚埋深、土体不排水的抗剪强度呈正相关。相比于现场试验或模型试验，数值分析方法可更细致直观研究鱼雷锚的承载能力。目前鱼雷锚承载力数值分析方法主要为有限元法，包括常规小变形有限元方法和可以考虑大变形的几何非线性有限元方法，如CEL、ALE等。Chen[5]通过ABAQUS有限元方法对动力锚在15°～80°不同拉拔倾角下承载力，并分析其影响因素。瑜璐等[6]应用ABAQUS软件分析特定埋深下鱼雷锚的抗拔承载力，探讨了锚型、土体类型、拉拔荷载倾角、拉拔荷载水平分量与锚翼夹角等多种因素对拉拔承载力的影响并得出规律性结论，最终提出预测鱼雷锚承载力的归一化V-H包络公式。Raie[7]通过CEL方法分析研究鱼雷锚安装成功后鱼雷锚周围地基土体二次固结情况，结果直观反映了超孔隙水压力的消散速率以及土体强度恢复情况。关于接触问题，CEL方法无法直接模拟正常固结土中锚-土之间的摩擦作用，大多数学者采用True提出的牛顿第二定律来模拟锚-土之间的相互作用[8]，但该方法无法真正反映土体的率效应和应变软化效应对摩擦力的影响。

鱼雷锚安装后的承载力对海上工程设计施工至关重要，目前研究表明，鱼雷锚的承载力还缺乏有效的数值模拟方法，这导致鱼雷锚的研究还很不成熟。物质点易于追踪物质运动及边界情况，规避了网格大变形问题，在解决鱼雷锚上拔过程土体大变形问题上具有十分明显的优点。针对上述问题，在已有研究基础上，通过三维显式物质点法开源程序[9]模拟鱼雷锚在土体中上拔过程，采用接触界面法向量与API相结合的接触算法模拟锚-土之间摩擦力，并且通过对现场原位试验结果进行可靠性分析，与现场试验结果相互验证补充，为鱼雷锚上拔承载力研究提供一种新的数值算法，此外，还对上拔过程周围土体变形规律和上拔承载力影响因素进行分析，为工程实际提供指导性建议。

1 物质点法方法基本理论

1.1 物质点离散

图1 物质点法示意图[10]Fig.1 Schematic diagram of the material point method[10]

在物质点法(MPM)中，连续体和其运动空间分别用物质点和网格两套不同的系统表示，如图1所示。MPM先将连续体离散为一组质点，物体材料区域内的物质信息均由质点携带，并由质点代表其所在区域，在每个计算时间步内，质点和计算背景网格不发生相对运动，质点会根据具体受力情况在背景网格中运动，背景网格只用于求解动量方程和计算空间导数，计算网格结点与质点间通过单值映射函数进行运动量信息传递。

MPM将连续体离散后，其密度为

(1)

式中：mp是质点p的质量；np为离散质点总数；δ(x)是Dirac Delta函数；xip是质点p的坐标。

1.2 动量方程及其求解

在不考虑热量交换情况下，连续体更新拉格朗日的动量方程和边界条件分别为

(2)

(3)

为求解动量方程，取虚位移δuj∈R0和权函数R0={δuj|δuj∈C0，δuj|Γu=0}，结合给定面力边界条件(3)式，得出动量方程(2)式等效积分弱形式为

(4)

将(1)式代入(4)式中，得到连续体离散后动量方程等效积分弱形式

(5)

由于应用物质点法求解动量方程时，质点和背景网格结点之间信息传递是通过在背景网格结点上建立的有限元函数NI(Xi)来实现的。因此，质点p的位移uip和虚位移δuip为

uip=NIpuiI

(6)

δuip=NIpδuiI

(7)

式中：uiI为结点位移；下标I、J、K表示网格结点。将(6)、(7)式代入(5)式中，得到背景网格结点的运动方程

(8)

(9)

(10)

网格结点I在i方向的动量piI为

(11)

式中：网格质量矩阵mIJ为

(12)

应用集中质量阵(ng为结点总数)

(13)

可以将结点运动方程(8)式简化为

(14)

之后，利用中心差分法对(14)式进行求解，得到当前时间步背景网格结点的位移增量，再借助形函数插值映射到质点，得到质点相关的其他物理量。运动方程求解完成，计算步结束，抛弃变形后的网格。

1.3 接触算法

解决无滑移黏着接触问题，标准MPM不需采取额外的处理，当需要考虑对象间的相对滑动和分离时，需要将接触算法引入到MPM中。采用接触界面法向量计算方法[10]与API方法结合模拟锚-土间摩擦作用力。

通常来说，当两个物体r和s间满足接触判定条件式(15)时

(15)

考虑到两物体法向单位向量不完全共线时，会导致动量不守恒和界面穿透情况，故两物体外法向单位向量选取至关重要。当物体r比物体s更硬时，应选取物体r来计算接触面的公法线方向，即取

(16)

(17)

(18)

(19)

(20)

1.4 材料模型

物质点法中，本构方程在质点上计算。由于钢的强度和模量远大于土体，鱼雷锚的本构模型采用线弹性模型，其应力应变关系为

(21)

(22)

土体本构模型采用理想弹塑性模型，其破坏形式包括剪切破坏和拉伸破坏。土体剪切失效模型采用Drucker-Prager模型，其屈服函数为

(23)

式中：J2为偏应力张量的第二不变量；kΦ、qΦ与材料的黏聚力c和摩擦角φ的关系由下式确定

(24)

(25)

式中：加号表示D-P屈服面在π平面上内接Mohr-Coulomb屈服面，减号则表示外接，计算采用内接圆来进行计算[11]。土体拉伸失效屈服函数为

ft=σm-σt

(26)

式中：σt为材料抗拉强度。剪切失效和拉伸失效相应的流动法则均采用关联流动法则。

为了保证应力率不受刚体转动的影响，每个物质点的应力率都采用焦曼应力率进行应力控制。剪切失效和拉伸失效相应的流动法则均采用关联流动法则。

图2 鱼雷锚上拔示意图Fig.2 Schematic diagram of torpedo anchor pulling up

2 鱼雷锚上拔过程物质点法分析

以巴西近海坎普斯湾鱼雷锚上拔现场原位试验[1]为研究对象，通过物质点法建立三维模型对其进行模拟，从而研究无翼鱼雷锚在饱和黏土中上拔承载力。

2.1 分析模型

巴西近海坎普斯湾现场试验[1]将在正常固结黏土中的泥线角θ0= 0°的无翼鱼雷锚从20 m深处拉出，但并未介绍锚链拉出角θa，拉出角θa定义为负载方向和水平方向之间的角度，考虑到现场试验实际情况，与竖向拉出情况类似，故取尾翼拉出角θa=90°(θ0，θa如图2所示)。

图3 鱼雷锚几何尺寸Fig.3 Torpedo anchor geometry

巴西近海坎普斯湾现场试验[1]的T-40鱼雷锚，直径D=0.76 m，长L=12 m，锚干重为240 kN，无尾翼，鱼雷锚形态及其参数如图3所示，下端部为圆锥形状，在模型中其离散间距为0.3 m的物质点。考虑到鱼雷锚上拔影响范围，土体模型尺寸取10 m×10 m×60 m，离散间距为0.5 m的物质点。整个模型总共55 390个物质点。模型计算区域为10 m×10 m×60 m，网格单元类型为八节点六面体，网格间距为0.5 m。鱼雷锚和土体物质点离散模型如图4所示。

图4 鱼雷锚-土物质点离散模型Fig.4 Discrete model of torpedo anchor and soils

鱼雷锚材料为钢，其本构参数见表1。结合已有的研究结果[12]，土体为软黏土，其内摩擦角为0°，不排水抗剪强度Su取(5+2z)kPa，其中z为地基深度，弹性模量E0取500Su。土体本构参数见表2。

表1 鱼雷锚本构参数Tab.1 Torpedo anchor constitutive parameters

表2 土体本构参数Tab.2 Soil constitutive parameter

模型初始状态为鱼雷锚埋入土体中，鱼雷锚的初始埋深H为20 m[12]。边界条件施加在计算网格上，4个侧面和底面均为固定。

外部荷载方面，鱼雷锚上拔模拟是通过施加竖直向上恒定锚链拉力将鱼雷锚拔出，阻力仅考虑侧摩阻力。接触选取鱼雷锚作为主体、主平面，取鱼雷锚的表面法线为公法线，摩擦系数μ取0，取API公式计算由土体粘聚力产生的侧摩阻力作为切向强度。

API规范[2]中规定对于黏土中嵌入式桩，沿桩长向侧摩阻力Fs为

Fs=f(z)As=αsuAs

(27)

式中：不排水抗剪强度Su取锚-土接触点埋深处对应值，界面摩擦比α取灵敏度的倒数1/St，仍采用已有的研究结果St取3，即α取1/3[12]。

算法采用显示积分格式，计算时间步长0.167E-04，计算时间为2.0 s，共60 010步。

图5 鱼雷锚竖向上拔阻力-归一化位移关系曲线Fig.5 Vertical pull-up resistance-normalized displacement relationship curve of torpedo anchor

2.2 鱼雷锚上拔承载力验证

在巴西近海坎普斯湾进行的无翼鱼雷锚现场试验中，平均贯入深度为20 m。如图3所示鱼雷锚的设计承载力为1 400 kN[1]，现场原位试验实测极限承载力在900～1 100 kN范围内变化，平均极限承载力为1 000 kN。根据现场原位试验数据，对埋深H=20 m的无翼鱼雷锚进行物质点法上拔数值模拟，得到鱼雷锚竖向上拔阻力与归一化位移关系曲线如图5所示。通过计算发现，锚归一化位移(锚尖位移与锚尖初始埋深比值)U/H为0.08～0.16范围内竖向上拔阻力达到最大值，故取位移达0.08H～0.16H所对应的承载力平均值为竖向抗拔承载力，即最终竖向上拔承载力为1 187 kN，处于现场试验测得的抗拔承载力范围内，与现场平均极限承载力相差18.7%，物质点法数模结果与现场测试数据的合理一致性证实了该物质点法数模在评估鱼雷锚在黏土中的抗拔承载力的可靠性。

2.3 上拔过程地基变形特性

为便于观察研究结果，其中计算输出区域为锚两侧各6D，锚在土中运动时深度向取26D，锚拔出土体时深度向取32D。

物质点法可以通过土体质点位移变化，得出土体颗粒行动方向，进而得出某一时刻土体具体变形情况。取锚拔出4个典型时刻土体竖向位移进行分析，如图6所示(图中位移正值表示位移向上)。由图6可看出，锚上拔过程地基土体竖向位移以向上为主，锚尾部附近地基土体变形大于锚侧附近地基土体，地基土体产生向上的竖向位移的范围为以初始锚尖埋深为长轴、12D为短轴的半椭圆形范围；地基土体产生向下的竖向位移的范围集中在初始锚尖埋深下方6D，且距锚轴线位置2D～4D范围，范围几乎不变。由于锚附近地基土体受到挤压作用，从径向向外扩张、流动，距锚轴线越远，受影响程度越小，故远锚侧土体竖向位移也就越小。此外，锚上拔过程中带动地基土体运动，运动轨迹产生的空腔在锚完全拔出后未被土体完全填充。由图6-c～图6-d可看出，锚虽然在不断上拔，但从锚体部分离开土体开始，地基土体扰动区域几乎没有发生改变，且竖向位移数值基本一致。由图6可看出，锚上拔完成变形土体区域为以初始锚尖埋深下方6D为高，12D为底边的等腰三角形范围。

6-a H=18 m6-b H=12 m6-c H=4 m6-d H=0 m图6 地基土体竖向位移等值线云图(单位:cm)Fig.6 Contour cloud diagram of vertical displacement of foundation soil

上拔过程地基土体累计等效塑性应变分布如图7所示。为便于观察锚周土体塑性区分布，土体范围取锚水平向左右6D，深度向取海床表面以下40D范围。黏土室内三轴试验累计轴向应变达到15%即认为其产生破坏，有研究指出冲击荷载作用下软黏土室内三轴试验累计轴向应变达到2%时土体产生屈服[13]，故认为累计等效塑性应变大于10%为土体塑性变形显著区域。由图7可看出，随着锚不断上拔，累计等效塑性应变区域由锚身开始逐渐向周围扩散，且变形范围不断增加，部分锚体拔出土体后变形区域基本不变。此外，随着鱼雷锚不断上拔，除了锚上拔路径周围土体产生变形外，海床表面两侧土体同样产生较大范围变形，鱼雷锚回收利用时要考虑土体变形范围，以防锚轨迹空洞被两侧土体逐渐填补出现土体塌陷，从而对其周围设施产生影响。

对于变形区域，由图7可看出，锚侧附近地基土体向两侧扩散，径向距离锚侧越近变形越大，直到部分锚体离开土体，锚上拔路径空腔区域出现缩小情况(图7-c)，主要原因是空腔由于深度过高，土体在重力作用下有向腔内垮塌，导致腔体侧壁中下部累计等效塑性应变增大。由图7-c、图7-d中可以看出，鱼雷锚拔出土体后土体塑性区无明显变化。图7结果表明，锚上拔显著塑性变形区为以锚尖初始埋深下5～6倍锚径位置的地基埋深为顶点、以地基表面12倍锚径为底边的等腰三角形。

7-a H=18 m7-b H=12 m7-c H=4 m7-d H=0 m图7 累计等效塑性应变云图Fig.7 Cumulative equivalent plastic strain cloud diagram

3 竖向上拔承载力影响因素分析

鱼雷锚贯入地基后的承载力与锚的质量和几何形状、土体强度分布和埋深有关，又考虑到土体强度分布和埋深对鱼雷锚安装后承载力起到关键性的作用，通过物质点法数值模拟，得出埋深和土体强度对竖向上拔承载力的影响。

3.1 埋深对竖向上拔承载力的影响

保持鱼雷锚及土体参数与现场原位试验数据[1]相同，现场原位试验中平均贯入深度为20 m，在典型的深水黏土沉积物中，锚埋深可达到锚长的3倍，竖向单调上拔承载力通常小于锚干重的5倍，此外有数值研究[4]表明，锚的埋深为锚长度的1.5～2.6倍，锚竖向上拔承载力在锚干重的2.4～4.1倍，故埋深分别选取20 m、25 m、30 m、35 m，分别对应锚长的1.7倍、2.1倍、2.5倍、2.9倍。为深入研究埋深对竖向上拔承载力的影响，模拟鱼雷锚在正常固结土(Su=5+2zkPa)中上拔，得到上拔过程中不同埋深下鱼雷锚归一化上拔阻力随归一化位移(锚埋深与参考埋深H0=20 m之比)变化曲线如图8所示。由图8可看出，随埋深增加，鱼雷锚上拔阻力稳定所需归一化位移略有增加，说明埋深增加时，土体强度越大，土体对鱼雷锚所产生的阻力相应较大。

鱼雷锚归一化竖向上拔承载力随埋深比(初始锚尖埋深与锚长之比)变化曲线如图9所示。由图9可看出，在实验范围内，鱼雷锚归一化竖向上拔承载力与埋深比近似成线性关系，线性比约为8.34，二者线性化公式为

(28)

在实验范围内，鱼雷锚埋深每增加5 m，相应竖向上拔承载力增加锚重的3.1～4.0倍，鱼雷锚竖向上拔承载力约为锚重的5.5～15.5倍。

图8 鱼雷锚上拔阻力-位移关系曲线Fig.8 Torpedo anchor normalized uplift resistance-normalized displacement relationship curve图9 鱼雷锚竖向上拔承载力随埋深比变化曲线Fig. 9 Curve of normalized vertical pull-out bearing capacity of torpedo anchor with depth-to-depth ratio

图10 鱼雷锚归一化上拔阻力-归一化位移关系曲线Fig.10 Torpedo anchor normalized uplift resistance-normalized displacement relationship curve图11 鱼雷锚归一化竖向上拔承载力随k值变化曲线Fig.11 Curve of normalized vertical pull-out bearing capacity of torpedo anchor with k value

3.2 土体强度对竖向上拔承载力的影响

为了研究不同性质黏土对竖向上拔承载力的影响，保持上述鱼雷锚和土体参数不变，埋深H=20 m不变情况下，不排水抗剪强度Su取(5+kz)kPa，其中k分别取0、1、2、3，得出鱼雷锚在不同不排水抗剪强度Su土体中归一化上拔阻力与归一化位移关系曲线如图10所示。

由图10可见，随k值增加，土体强度增加，土体对鱼雷锚所产生的阻力相应较高，且随k值增加，鱼雷锚上拔阻力增加阶段越来越长，即所需位移越来越大。

如图11为不同k值下鱼雷锚归一化竖向上拔承载力随k值变化曲线。由图11可见，归一化竖向上拔承载力与k值呈正相关，二者关系可拟合为

Q=(1.305k+0.697 9)WA

(29)

4 结论

应用物质点法对室内模型试验和巴西近海坎普斯湾[1]原位试验中埋深20 m的无翼鱼雷锚上拔过程进行物质点法分析，得到鱼雷锚竖直上拔承载力，并进一步研究上拔过程地基土体变形规律、塑性变形区范围与竖向上拔承载力影响因素，得出以下结论：

(1)物质点法数模结果与现场试验结果之间的比较表明，应用物质点法研究鱼雷锚上拔过程规律是可靠的、有效的。

(2)锚上拔过程地基土体竖向位移以向上为主，锚尾部附近地基土体变形大于锚侧附近地基土体，竖向位移向上地基土体范围为以初始锚尖埋深为长轴，12D为短轴的半椭圆形范围；竖向位移向下地基土体集中在初始锚尖埋深下方6D，12D为底边的等腰三角形范围，范围几乎不变。

(3)锚上拔显著塑性变形区为以锚尖初始埋深下5～6倍锚径位置的地基埋深为顶点、以地基表面12倍锚径为底边的等腰三角形。因此，在鱼雷锚上拔回收过程中，需要考虑该范围内其对其他结构物的影响。

(4)鱼雷锚竖向上拔承载力影响因素分析结果表明，试验范围内，埋深、土体强度与竖向上拔承载力呈正相关。在实验范围内，不同埋深下，鱼雷锚归一化竖向上拔承载力与深长比近似呈正相关；不同土体强度下，鱼雷锚归一化竖向上拔承载力与k值呈正相关。