随着人类对海洋石油天然气的需求越来越高,近海水域的油气储量逐渐减少,海洋油气开采逐步走向深海。钢悬链式立管(SteelCatenaryRiser,SCR)具有成本低、与平台连接灵活和耐高温高压等优点,在海洋工程领域应用广泛。钢悬链式立管分为悬垂段和流线段,悬垂段与流线段的过渡区域称为触地区域。在波浪等外载荷作用下,浮式平台发生周期性运动,带动悬链式立管运动,针对平台和立管的耦合运动,国内学者提出了整体分析模型[1]和二维动力学模型[2]。由于立管的运动,使得触地区域与海床频繁地发生相互作用,对立管的动力响应和疲劳寿命影响显著。因此,准确模拟悬链线立管管土作用过程是国内外研究热点之一。

国内外已有许多学者对管土作用问题进行了数值模拟研究。在众多数值模拟方法中,线性模型相对简单并且易于实现数值模拟。王坤鹏等[3]将管土作用简化为三种不同的线性模型:线性刚度、线性截断、线性截断并考虑海床吸力和刚度衰减,并且通过数值模拟,认为考虑了吸力效应后立管会有较多的疲劳应力循环。由于管土作用过程是一个非线性过程,利用线性模型进行模拟与实际存在一定误差。Aubeny等[4]以及Randolph等[5]分别提出了两类非线性管土作用模型,以下分别简称为A-B模型和R-Q模型。目前国内对管土作用的研究多采用A-B模型[6,7,8,9,10]。相较于A-B模型,R-Q模型的管土作用模式更丰富,其P-y曲线的参数更多,能够更加真实地模拟管土作用过程,目前国外学者对于SCR管土相互作用的研究主要采用R-Q模型。Hejazi等[11]提出了“等效线性土体刚度”的概念,通过设置线性土体刚度使立管的应力水平与利用R-Q模型计算结果相当,从而达到简化非线性管土模型的作用;Shiri[12]利用R-Q模型描述非线性海床与立管的相互作用,研究了顶端垂荡激励对SCR疲劳性能以及应力分布的影响;Muraleedharan等[13]和Dong等[14]分别利用R-Q模型研究了立管的疲劳损伤和立管的节点响应,均得出R-Q模型较为保守的结论。

本文利用ABAQUSSubroutine中的用户自定义单元(UserDefinedElement,UEL)基于R-Q模型开发悬链式立管触地单元,通过与模型试验结果比较分析,验证开发的数值分析模型的合理性,通过与线模型模拟结果对比,讨论非线性R-Q模型体现出的优势,最后进行参数敏感性分析,探讨立管顶端垂荡幅值及周期对沟槽形状、立管位移及弯矩响应等的影响。

1、数值模型

1.1 管土作用模型

R-Q模型是一种非线性管土作用模型,其管土作用过程可分为7种模式:分别为①初始嵌入;②阻力减小的提升;③吸力产生的提升;④吸力衰减的提升;⑤吸力释放的再次嵌入;⑥继续再次嵌入;⑦土体分离后再次嵌入。土体反力P与嵌入深度y的关系曲线表达式如式(1)～(7)所示,图1为该模型的P-y曲线。

极限抗力:

Pu(y)=NcD(S0+Sgy)         (1)Ρu(y)=ΝcD(S0+Sgy)         (1)

式中:Nc为与y有关的无因次剪切因子;D为立管的外直径;S0为海床表面的剪切刚度;Sg为海床的剪切刚度梯度。

极限吸力:

Pu−suc(y)=−fsucPu(y)         (2)Ρu-suc(y)=-fsucΡu(y)         (2)

式中,fsuc为吸力因子。

初始嵌入模式:

PIP(y)=ζ1+ζPu(y)         (3)ΡΙΡ(y)=ζ1+ζΡu(y)         (3)

式中,ζ为与y有关的无因次嵌入因子。

提升模式,取式(4)和(5)中的较大值

PUL(y)=P0−ζ0−ζAUL(y)+(ζ0−ζ)(P0−Pu−suc(y))         (4)Pmin(y)=EUL(y)Pu−suc(y)         (5)ΡUL(y)=Ρ0-ζ0-ζAUL(y)+(ζ0-ζ)(Ρ0-Ρu-suc(y))         (4)Ρmin(y)=EUL(y)Ρu-suc(y)         (5)

式中:P0与ζ0为模式开始时的土体反力与嵌入因子;AUL与EUL为无因次参数。

再次嵌入模式,取式(6)和(7)中的较小值

PRP(y)=P0+ζ−ζ0ARP(y)+(ζ−ζ0)(Pu(y)−P0)         (6)Pmax(y)=ERP(y)PIP(y)         (7)ΡRΡ(y)=Ρ0+ζ-ζ0ARΡ(y)+(ζ-ζ0)(Ρu(y)-Ρ0)         (6)Ρmax(y)=ERΡ(y)ΡΙΡ(y)         (7)

式中,ARP与ERP为无因次参数。

本文同时选取线性模型与非线性R-Q模型进行对比。线性模型的P-y曲线如图2所示,其管土作用模式可分为初始嵌入、线性提升、管土部分分离、管土分离和再次嵌入5个部分。

线性模型中的骨干曲线如式(1)所示,y1、P1由骨干曲线确定,P2、y3、y2的值由式(8)～(10)确定

P2=−fsucP1         (8)y1−y3=μsepD         (9)(y2−y3)/(y1−y3)=λsuc         (10)Ρ2=-fsucΡ1         (8)y1-y3=μsepD         (9)(y2-y3)/(y1-y3)=λsuc         (10)

图1R-Q模型P-y曲线

图2线性模型P-y曲线

式中:μsep为管土分离的最大距离因子;λsuc为最大吸力出现位置因子。

得到(y1,P1)、(y2,P2)、(y3,P3)三个关键点之后,通过线性连接即可得到相应的P-y曲线。

1.2 基于UEL的管土接触单元开发

本文基于通用有限元软件平台ABAQUS模拟管土相互作用,采用梁单元模拟立管,自定义触地单元模拟海床。触地单元实质为非线性弹簧,弹簧刚度根据P-y曲线得到。建模示意图如图3所示,在梁单元节点处定义触地单元。借助ABAQUSSubroutine功能进行二次开发,使用FORTRAN语言编写用户定义单元子程序(UEL)。在编写UEL程序时,需对右手端矢量RHS和刚度矩阵AMATRX两个变量进行明确定义,通过SVARS动态数组保存用户定义单元求解过程中依赖的状态变量值。

图3立管模型与触地单元连接示意图

2、试验结果验证

本文选取Wang等[15]在大型3D土壤/结构相互作用测试系统开展的模型试验研究进行数值分析方法合理性验证。试验立管参数如表1所示。

该试验中,立管顶端连接致动器,控制其垂向运动,立管底端自由。垂向激励幅值为48mm(0.3D),频率0.2Hz,循环总数200次。试验记录和比较了不同监测点的位移和弯矩响应,位移监测点DW1～DW4分别距立管顶端3m、4m、5m、6m,弯矩监测点BM1～BM4分别距立管顶端1.5m、2.5m、3.5m、4.5m。

图4为4个位移监测点200个周期下的位移变化曲线,中部为试验结果,左侧和右侧分别为R-Q模型和线性模型模拟结果。结果表明,R-Q模型能够给出较为合理的数值结果,并且相较线性模型与试验结果更为吻合。观察R-Q模型给出的数值结果可以发现,4个位移监测点在200个循环加载下嵌入深度逐渐加深,这与试验结果具有相同的规律,而线性模型模拟下,4个位移监测点在循环次数大于10次后嵌入深度便不再加深。出现这一现象的根本原因在于线性模型的P-y曲线仅取决于最大深度点的P、y值(即图2中的y1,P1),若立管在提升过程中发生反向运动,P-y曲线将按原路返回,所以数值模拟中易发生P-y曲线重叠现象,使得立管不再嵌入;而对于R-Q模型,P-y曲线不仅取决于y1,P1值,当立管发生反向运动时,需记录此时的P、y值,作为下一段P-y曲线的重要参数,所以在数值模拟过程中P-y曲线几乎不可能重叠,从而使立管能够以合理的速度继续嵌入土体,这也体现了R-Q模型相较线性模型的优势。

图5为立管嵌入深度对比情况,图中散点表示试验中4个位移监测点在特定时刻的嵌入深度,曲线为数值模拟给出的立管形状,实线和虚线分别为R-Q模型和线性模型结果。由图5可以更直观的看出,线性模型结果中,立管在N=10之后嵌入深度不再增加,R-Q模型的立管嵌入深度和嵌入速度均更为合理,与试验结果相比也更为吻合。

图4两种模型的位移变化曲线与试验结果对比

图6为4个弯矩监测点的弯矩随循环次数的变化曲线,中部为试验结果,左侧和右侧分别为R-Q模型和线性模型模拟结果。无论是线性模型还是R-Q模型,在弯矩的量值方面与试验值相差不大,但是在变化趋势上,试验中4个监测点的弯矩均值在200个循环中有缓慢增加的趋势,R-Q模型给出的结果能够很好地模拟出这个趋势,而线性模型则无法体现出这一变化过程。

图5立管嵌入深度比较

总体而言,无论从位移还是弯矩角度,R-Q模型模拟管土作用的效果均较为合理,并且比线性模型更接近于实测结果。

3、顶端垂荡对立管动力响应的影响

3.1 计算工况

本文选取的立管参数和土体参数如表2所示。

顶端运动条件对SCR动力响应影响显著。本节针对顶端垂荡幅值A和周期T开展参数敏感性分析,A分别取2m、4m和6m,T分别取5s、10s、15s、20s、25s和30s,共计18个工况。有限元模型中以海床锚点作为坐标原点。

图6两种模型的弯矩变化曲线与试验结果对比

3.2 立管不同位置处动力响应

在动力分析过程中,立管触地点并非一个固定的点,而是会随着时间不断变化。以A=6m,T=30s工况为例,TDP的位置如图7所示。65s时TDP位于管长730m处,而80s时TDP则位于管长764m处。

对于相同的顶端垂荡激励,不同的立管位置可能会有截然不同的动力响应。根据其不同的响应特性,立管节点可分为三类,第一类是锚点至TDP之间的节点,始终紧贴着海床沟槽,第二类是TDP附近的节点,时而紧贴沟槽时而与沟槽分离,第三类是TDP至立管顶端之间的节点,从不与海床沟槽接触。就A=6m,T=30s工况而言,管长720m、750m以及780m处的节点分别对应上述三类节点,三个节点的位移时历和P-y曲线比较如图8和图9所示,其中780m处的节点由于未与海床接触故不存在P-y曲线。

图7不同时刻TDP位置

图8位移时历比较

图9P-y曲线比较

由图8和图9可知,由于720m处的节点始终紧贴着沟槽,位移受到约束,所以位移响应的幅值较小,其对应的P-y曲线范围也较小,没有出现土体吸力,同时其嵌入土体的速度比较缓慢;而750m处的节点处于悬垂段到触地点之间,运动相对自由,所以位移响应的幅值较大,同时由于土体反力的存在,其位移波形的谷底较平坦,其对应的P-y曲线完整,包括了阻力减小、吸力产生和土体分离等部分,其嵌入土体的速度相较720m处的节点更快;对于780m处的节点,由于未受到管土作用的影响,其位移响应遵循标准的正弦变化曲线。

由于第二类节点的管土作用模式最为复杂,嵌入土体深度较深且嵌入速度也较快,因此下文将主要讨论顶端垂荡激励对第二类节点的影响。

3.3 顶端垂荡激励幅值的影响

选取垂荡激励周期T=30s,激励幅值A=2m、A=4m、A=6m三个工况,进行对比分析。

3.3.1 沟槽形状及其发展速度

图10为不同工况下沟槽发展过程,N为循环加载次数,N=0表示静力模拟的结果。顶端垂荡激励幅值越大,在相同的循环次数时形成的沟槽越宽越深。同时,激励幅值A=2m、A=4m、A=6m的工况对应的沟槽最深点的位置分别为距锚点742m、741m和740m处,即沟槽位置略有向锚点靠近的趋势。

定义沟槽深度增量Δyi=yi-yi-1,其中yi表示循环次数N=i时刻沟槽最大深度,y0则表示静力模拟的沟槽最大深度。定义沟槽相对发展速度αe=Δyi/yi-1,αe值越大则说明一个循环下的沟槽深度的相对增量越大,从而沟槽发展得越快。图11和图12分别为不同激励幅值下沟槽深度和沟槽发展速度对比情况。由图11可以看出随着循环次数增加沟槽的深度不断增加,而图12则表明沟槽发展速度随循环次数增加而逐渐减小。在相同的循环次数下,激励幅值A越大,沟槽发展速度也越快。

图10不同激励幅值沟槽形状比较

3.3.2 位移和弯矩

图13和图14分别为不同激励幅值下740m处节点位移时历图和P-y曲线图。顶端垂荡激励幅值A越大,立管节点的垂向位移幅值也越大,从而节点P-y曲线的范围也越大。A=2m与A=4m工况仅有极少时间立管受到海床吸力作用,而对于A=6m工况,由于P-y曲线的范围大,立管受到海床吸力作用的时间也较长,使立管嵌入深度更深,这也导致顶端大幅值激励下沟槽较深且发展速度较快。

图11沟槽最大深度比较

图12沟槽相对发展速度比较

图13740m处节点的位移时历比较

图14740m处节点的P-y曲线比较

图15为不同激励幅值下立管弯矩分布。顶端垂荡激励幅值越大,立管最大弯矩的绝对值越小。由于沟槽区域相当于立管悬垂段与拖地段之间的缓冲区,结合上述沟槽深度与激励幅值的关系,可以得出这样的规律:顶端垂荡激励幅值越大,形成的沟槽就越深,从而立管的缓冲区越长,立管曲率变化就相对缓和,因此立管弯矩也就有所减小。

图15不同激励幅值弯矩比较

总的来说,顶端垂荡激励幅值越大,则节点的位移幅值也越大,从而形成的沟槽就越深,沟槽发展速度也越快,同时最大弯矩的绝对值越小。

3.4 顶端垂荡激励周期的影响

选取激励幅值A=2m,激励周期T=5s、T=10s、T=15s、T=20s、T=25s、T=30s六个工况。

(1)沟槽形状及其发展速度

图16为不同工况下沟槽发展过程。激励周期T越小,沟槽的深度越深,T=5s工况的沟槽深度变化最大。T=5s工况沟槽最深点出现在距锚点744m处,其余工况的沟槽最深点均出现在距锚点742m处。由此可知,激励周期对沟槽最深点的位置影响不大。

图17和图18分别为不同激励周期下沟槽深度和沟槽发展速度对比情况。可以看出,激励周期T越小,沟槽越深且发展速度越快。这一规律在周期较小时(即T<15s)表现得尤为显著。

(2)位移和弯矩

图19和图20分别为不同激励周期下742m处节点的位移随循环次数的变化曲线图和P-y曲线图。顶端垂荡激励的周期越小,立管节点的垂向位移幅值越大,特别是T=5s的工况,由于顶端垂荡激励周期很小(激励程度剧烈),使得在相同的顶端幅值条件下,各立管节点的运动速度更快,所以在相同的时间内可以嵌入更深的深度。

图21为不同激励周期下立管弯矩分布情况。可以看出,顶端垂荡激励周期越小,立管最大弯矩的绝对值也越小。这与不同幅值工况的结果类似,由于激励周期减小,沟槽的深度增加,从而立管曲率变化相对缓和,立管弯矩有所减小。

图16不同激励周期沟槽形状比较

图17沟槽最大深度比较

图18沟槽相对发展速度比较

总的来说,顶端垂荡激励的周期越小,在同样的循环次数下,管土作用形成的沟槽就越深,沟槽发展速度也越快,而立管最大弯矩的绝对值越小。

图19742m处节点的位移变化曲线比较

图20742m处节点的P-y曲线比较

4、结论

本文以ABAQUS梁单元为基础,利用用户自定义单元(UEL)实现了基于R-Q模型的触地单元的数值建模,从而进行了非线性管土作用过程的时域模拟。通过与试验结果进行比较,验证了建立的数值分析模型的有效性,并讨论了其相较线性模型的优势。最后探讨了立管顶端垂荡幅值及周期对沟槽发展、立管位移及弯矩响应等的影响,得到以下结论:

图21不同激励周期的立管弯矩比较

(1)论从沟槽发展还是位移及弯矩响应角度,采用非线性R-Q模型模拟管土相互作用均能得出较为合理的预报结果,并且相较传统线性模型更接近于试验观测情况,优势明显。

(2)管土相互作用过程中,立管的不同位置处可能会有截然不同的动力响应,其中触地点附近节点的管土作用模式最为复杂,响应的非线性特征最为显著,立管嵌入土体的深度也较深;

(3)端垂荡幅值越大或周期越小(即顶端激励越剧烈),立管节点的位移响应幅值越大,形成的沟槽越长也越深,沟槽发展速度越快。顶端垂荡幅值对管土作用的影响效应相较垂荡周期更为显著。

(4)槽的发展速度在循环一开始时较快,并且随着加载次数增加逐渐减小,沟槽逐步趋于稳定。

(5)着顶端垂荡幅值的增加,沟槽最深点出现的位置略有向锚点靠近的趋势,这一趋势在顶端垂荡周期变化时表现得不明显。

(6)槽是立管与海床之间的缓冲区域,沟槽的形成有助于减小立管弯矩,在顶端垂荡幅值和周期变化时,沟槽越深立管最大弯矩的绝对值越小。

参考文献：

[1]孟庆飞,黄维平,刘建军.深水钢悬链式立管与浮式平台整体分析方法研究[J].振动与冲击,2013,32(17):19-23.

[2]孟丹.钢悬链线输流立管顶部浮体激励非线性响应研究[J].振动与冲击,2013,32(4):96-101.

[3]王坤鹏,薛鸿祥,唐文勇.基于海床吸力和刚度衰减模型的深海钢悬链线立管动力响应分析[J].上海交通大学学报,2011,45(4):585-589.

[6]白兴兰,黄维平,谢永和,等.基于非线性海床刚度模型的钢悬链线立管动力响应分析[J].振动与冲击,2015,34(21):82-88.

[7]周阳,杨超凡,黄维平.海床土刚度非线性的钢悬链式立管响应分析[J].哈尔滨工程大学学报,2017,38(3):356-362.

[8]陈振新,李捍平,李世强,等.考虑非线性管-土接触模型的钢悬链线立管触地区动态曲率分析[J].海洋工程,2018,36(6):77-83.

[9]常爽,黄维平,杨超凡.非线性管土作用下钢悬链式立管动力响应分析[J].中国海洋大学学报(自然科学版),2018,48(5):111-118.

郑孟添,袁昱超,薛鸿祥,唐文勇.浮体垂荡激励下深海悬链线立管触地端动力响应及沟槽发展研究[J].振动与冲击,2020,39(21):34-41+54.

基金:国家自然科学基金(51579146;51979166;51490674).