CN114330072A

CN114330072A - 一种考虑土体空间变异性的管道撞击损伤分析方法、设备

Info

Publication number: CN114330072A
Application number: CN202111669661.XA
Authority: CN
Inventors: 姜逢源; 赵恩金
Original assignee: China University of Geosciences
Current assignee: China University of Geosciences
Priority date: 2021-12-31
Filing date: 2021-12-31
Publication date: 2022-04-12

Abstract

本发明提供了一种考虑土体空间变异性的管道撞击损伤分析方法、设备，首先，根据海床土体尺寸和海床土体强度参数统计特征，将土体离散成一系列欧拉实体单元；通过Python脚本读取每个欧拉实体单元的形心坐标；然后生成三维土体强度平稳随机场，根据初始温度场进行叠加，进而得到三维土体强度非平稳随机场，通过子程序技术与大变形有限元数值模型进行耦合分析，得到N组管道凹痕变形值，然后进行分布拟合，得到管道损伤分析结果。本发明的有益效果是：针对大变形有限元数值模型中欧拉单元在空间保持固定不变的特点，通过设置温度场及子程序技术有效捕捉土体在计算过程中的位置，实现土体强度随空间变化的特征，分析得到的管道损伤结果与实际情况更接近。

Description

一种考虑土体空间变异性的管道撞击损伤分析方法、设备

技术领域

本发明涉及海底管道损伤领域，尤其涉及一种考虑土体空间变异性的管道撞击损伤分析方法、设备。

背景技术

海底管道是海洋油气开发产业的“生命线”，其作为连接海底设施和陆上终端的介质，承担着输送资源的重要任务，其结构的安全性是保障油气输运体系正常运行的关键。由于海底管道所处的特殊环境，其检修维护成本极大增加，海底管道一旦失效，将带来严重的经济损失和环境污染。其中，由人类第三方活动引发的坠物落水事故频有发生，使得结构物撞击附近海域的服役管道，造成管道的冲击损伤，进而引发管道的屈曲破坏，甚至直接导致管道破裂。

如何抵御坠物的冲击作用直接关系到管道的运营安全。实际工程中，常对管道进行埋深处理以保障其安全，已有研究对于埋置管道的撞击损伤分析尚处于初级阶段：DNV-RP-F107中仅考虑到碎石层对冲击能量的吸收作用，对于埋深并未说明。美国船级社规范中仅给出埋深的建议值。巴西、日本等国家对此提出的标准也不全面。我国海洋工程领域对于挖沟埋深的标准没有明确规定，设计工作主要凭经验完成。同时，土体材料由于其内部物质组成，应力条件、沉积条件、埋藏条件等因素，具有明显空间变异性特征，影响土体承载能力。但目前国内外的研究对象局限在均质土体阶段，忽略了土体的空间变异性，无法真实反映管道结构损伤响应。在近年来的海洋岩土工程研究中，土体的空间变异性因素已在桩基承载力、边坡稳定性等领域有了充分考虑，研究结果表明土体的空间变异性会影响土体的失效机理，明显降低土体的承载力，若忽略该因素得到的评估结果偏危险。对于管道受坠物冲击的损伤分析，有限元数值分析为最为有效、经济的研究手段。由于坠物侵入土体过程中会引发土体的大变形行为，传统的有限元分析方法会由于变形过大而出现精度低或不收敛等问题，需要利用大变形有限元分析技术建立相应的数值模型。因此，如何将土体三维空间变异性因素引入坠物-管道-土体大变形有限元数值模型中，捕捉管道结构的损伤行为是保障管道安全的关键技术。

综上所述，针对具有一定埋深的海底管道的受冲击损伤，提出一种基于三维随机场大变形有限元分析技术的分析方法，得到的评估结果更加符合实际情况，降低了工程设计中的风险。首先基于卡南-洛伊夫展开生成三维平稳随机场，随后通过初始温度场及子程序技术完成三维土体强度非平稳随机场的建立并将其映射至大变形有限元数值模型中，最后经过蒙特卡洛法抽样模拟得到管道在土体三维变异性条件下的统计响应特征。这一方法可有效完成大变形有限元模型的三维土体强度非平稳随机场的建立，分析考虑土体三维空间变异性因素下管道的损伤响应特征，为海底管道防冲击损伤提供有效理论依据。

发明内容

为了解决上述问题，本发明提供了一种考虑土体空间变异性的管道撞击损伤分析方法、设备，一种考虑土体空间变异性的管道撞击损伤分析方法，主要包括以下步骤：

S1：获取海床土体的尺寸和海床土体强度参数的统计特征；

S2：基于Abaqus平台，根据海床土体尺寸建立相应欧拉实体，并进行网格划分，得到一系列欧拉实体单元；

S3：通过Python脚本读取每个欧拉实体单元的形心坐标；

S4：基于海床土体强度参数统计特征，生成三维土体强度平稳随机场；

S5：将所述形心坐标代入至所述三维土体强度平稳随机场中，得到对应于三维土体强度平稳随机场；

S6：根据海床土体强度随深度增加的梯度，得到初始温度场；

S7：定义土体、管道、坠物的材料属性；在Abaqus平台中根据管道、坠物的尺寸建立相应的管道实体、坠物实体，并为管道实体、坠物实体划分网格；定义土体材料的本构模型为摩尔-库伦模型，管道的本构模型为各向同性硬化模型；定义分析步为显示动力学分析步；定义土体、坠物、管道彼此之间的接触条件为通用接触；定义坠物的初始速度；定义海床土体的边界条件为底面固定，四周侧向约束，定义管道两端的边界条件为固定约束，以上述信息作为输入量，建立相应的大变形有限元数值模型；

S8：根据所述三维土体强度平稳随机场和初始温度场进行叠加，得到三维土体强度非平稳随机场；

S9：以同样的方法得到N组三维土体强度平稳随机场以及N组三维土体强度非平稳随机场；

S10：通过子程序VUSDFLD将N组三维土体强度非平稳随机场与大变形有限元数值模型进行耦合并分析，输出得到N组管道凹痕变形值；

S11：对所述N组管道凹痕变形值进行分布拟合，得到海底管道受坠物撞击损伤的分析结果。

进一步地，，生成三维土体强度平稳随机场的方法采用的是卡南-洛伊夫展开法。

进一步地，所述三维土体强度平稳随机场为：

式中，S_u0(x,y,z,θ)为三维土体强度平稳随机场；λ_ln为对数正态分布的均值；ξ_ln为对数正态分布的标准差；M为卡南-洛伊夫展开的项数；i为第i项卡南-洛伊夫展开的下标；ξ_i(θ)为相互独立的标准正态随机变量；λ_i为自相关函数的特征向量；f_i(x,y,z)为自相关函数的特征函数。

进一步地，

其中，

为不排水抗剪强度S_u的均值，

为变异系数。

进一步地，所述三维土体强度平稳随机场为：

其中，(x_c,i,y_c,i,z_c,i)为欧拉实体单元的形心坐标；i为第i个欧拉实体单元的下标；λ_ln为对数正态分布的均值；ξ_ln为对数正态分布的标准差；M为卡南-洛伊夫展开的项数；j为第j项卡南-洛伊夫展开的下标；ξ_j(θ)为相互独立的标准正态随机变量；λ_j为自相关函数的特征向量；f_j(x_c,i,y_c,i,z_c,i)为对应于形心坐标的特征函数。

进一步地，所述初始温度场为：

其中，S_u0为海床表面z＝0处的土体不排水抗剪强度，

为土体强度随深度增加的梯度，z为海床土体深度。

进一步地，所述三维土体强度非平稳随机场为：

其中，S_u0为海床表面z＝0处的土体不排水抗剪强度，Temp_i表示第i个欧拉土体单元的初始温度场，S_u0(x_c,i,y_c,i,z_c,i,θ)为三维土体强度平稳随机场。

一种考虑土体空间变异性的管道撞击损伤分析设备，其特征在于：包括：处理器及存储设备；所述处理器加载并执行所述存储设备中存储的指令及数据用于实现所述的基于三维随机场大变形有限元分析技术的海底管道受坠物撞击损伤分析方法。

本发明提供的技术方案带来的有益效果是：针对大变形有限元数值模型中欧拉单元网格在空间保持固定不变的特点，通过设置温度场及子程序技术有效捕捉了土体在计算过程中的位置，实现了土体强度随空间变化的特征，且分析得到的管道损伤结果与实际情况更加接近。

附图说明

下面将结合附图及实施例对本发明作进一步说明，附图中：

图1是本发明实施例中一种考虑土体空间变异性的管道撞击损伤分析方法的流程图。

图2是本发明实施例中土体强度三维土体强度非平稳随机场实现的示意图。

图3是本发明实施例中管道凹痕统计特征的示意图，其中，图(a)为管道凹痕深度概率密度示意图，图(b)为管道凹痕深度累积概率示意图。

图4是本发明实施例中硬件设备工作的示意图。

具体实施方式

为了对本发明的技术特征、目的和效果有更加清楚的理解，现对照附图详细说明本发明的具体实施方式。

本发明的实施例提供了一种考虑土体空间变异性的管道撞击损伤分析方法、设备。

请参考图1，图1是本发明实施例中一种考虑土体空间变异性的管道撞击损伤分析方法的流程图，是针对于土体随机场建模绘制而成的，实施步骤可以分为三个部分，分别为“三维土体强度平稳随机场生成”、“三维土体强度非平稳随机场生成及大变形有限元数值模型建立”、“蒙特卡洛抽样模拟”。本实施例中有限元数值计算的建模及求解基于大型通用有限元程序Abaqus平台来完成。

各部分可细分为以下步骤：

1、三维平稳随机场生成

1.1确定海床土体的尺寸，包括长度L、宽度B和高度H；

1.2确定海床土体强度参数的统计特征，包括不排水抗剪强度S_u的均值

变异系数

和土体强度参数的自相关距离η_x，η_y，η_z。

1.3基于Abaqus平台，根据海床土体尺寸(L,B,H)建立相应欧拉实体，并采用结构化网格划分技术进行网格划分，将土体离散成一系列欧拉实体单元，单元种类为三维八节点缩减积分欧拉实体单元。

1.4通过Python脚本读取模型数据库中每个欧拉实体单元的形心坐标(x_c,i,y_c,i,z_c,i)，该模型数据库为网格划分之后生成的对应的欧拉实体单元组成的数据库，可以从中获取每个欧拉实体单元的形心。划分网格的目的就是为了获取每个欧拉实体单元的形心坐标。

1.5根据土体强度参数统计特征

及自相关距离(η_x，η_y，η_z)，基于卡南-洛伊夫展开法生成三维土体强度平稳随机场，一般情况下，土体强度参数遵循对数正态分布，其三维土体强度平稳随机场的表达式为如式(1)所示：

式中，S_u0(x,y,z,θ)为三维土体强度平稳随机场；λ_ln为对数正态分布的均值；ξ_ln为对数正态分布的标准差；M为卡南-洛伊夫展开的项数；ξ_i(θ)为相互独立的标准正态随机变量；λ_i为自相关函数的特征向量；f_i(x,y,z)为自相关函数的特征函数。

上述用于生成三维土体强度平稳随机场使用的方法还可以采用(1)协方差矩阵分解法、(2)移动平均法、(3)谱表现法、(4)转动带法等。

步骤1.5中生成三维土体强度平稳随机场的具体过程如下：

1.5.1确定三维自相关函数ρ_ln(x,y,z)的形式为指数形式，如式(4)所示，由于其为可分离形，三维自相关函数的特征向量及特征函数可由一维自相关函数对应的特征向量及特征函数相乘得到。

式中，x₁，x₂为x方向上任意两点的坐标；y₁，y₂为y方向上任意两点的坐标；z₁，z₂为z方向上任意两点的坐标；η_x表示x方向上的自相关距离；η_y表示y方向上的自相关距离；η_z表示z方向上的自相关距离

1.5.2通过理论解析法，求解一维随机场自相关函数ρ_ln(x)的特征向量及特征函数。通过求解式(6)所示的超静定方程，可得到有穷多个正实根w_x。将w_x代入式(7)、式(8)，即可求得ρ_ln(x)对应的特征向量λ_x,i和特征函数f_i(x)。

(η_x ²w_x ²-1)sin(w_xL_x)＝2η_xw_xcos(w_xL_x) (6)

f_i(x)＝a_icos(w_x,ix)+b_isin(w_x,ix) (8)

a_i＝η_xw_x,ib_i (9)

式中，ρ_ln(x)为一维自相关函数；L_x为一维随机场的区域长度；a_i为待定系数；b_i为待定系数；w_x,i为式(6)中第i个正实根。

1.5.3根据一维自相关函数的特征向量及特征函数，合成三维自相关函数的特征向量及特征函数。按照一维自相关函数的特征值乘积

进行任意组合，按特征值乘积降序排列，取其前M项，从而确定各一维特征向量对应的组合次序(i_x,j_y,k_z)。按照该组合次序，按式(11)-(12)进行相乘，得到三维自相关函数的特征向量及特征函数。

式中，

为在x方向上的一维随机场特征值；

为在y方向上的一维随机场特征值；

为在x方向上的一维随机场特征值；f_i(x,y,z)表示特征函数，

为x方向上的一维特征函数；

为y方向上的一维特征函数；

分别表示为z方向上的一维特征函数。

1.5.4将1.5.3步骤中得到的自相关函数特征向量λ_i及特征函数f_i(x,y,z)代入式(1)，可得到三维土体强度平稳随机场的表达式。

1.6将步骤1.4得到的欧拉单元形心坐标(x_c,i,y_c,i,z_c,i)代入式(1)中，即可得到对应于三维土体强度平稳随机场。

2、三维土体强度非平稳随机场的生成及大变形有限元数值模型建立

由于本发明针对管道受坠物损伤的分析方法中数值模型不光包括土体，还包括管道和坠物，故此处添加了关于管道以及坠物模型的建模描述，具体描述如下：

2.1确定土体本构模型，定义土体材料的强度参数，将土体的不排水抗剪强度S_u定义为以场变量FV为自变量的函数S_u(FV)。

2.2根据管道和坠物的几何信息、材料信息生成对应的实体，并对管道模型及坠物模型划分基于拉格朗日算法的网格。

2.3将管道、土体、坠物实体映射至装配空间。

2.4定义分析步，为显示动力学分析步，设定计算终止时间t_t。

2.5定义接触条件，为通用接触；

2.5定义边界条件，管道两端设置固定约束，土体底面为固定约束，四周表明为侧向位移约束；

2.6对土体施加初始地应力场，对坠物施加初始速度；

以上述信息作为输入量，建立相应的大变形有限元数值模型，其输出量为管道凹痕变形值，如果直接利用该大变形有限元数值模型，输出量就是没有考虑土体变异性的凹痕变形值，为了得到考虑土体变异性的管道凹痕变形值，就需要将大变形有限元数值模型与三维土体强度非平稳随机场进行耦合。

2.7确定土体强度随深度增加的梯度

对初始温度场采用解析式法将其定义为温度随深度z以梯度

线性增加的形式，如式(14)所示

式中，S_u0为海床表面z＝0处的土体不排水抗剪强度。

2.8定义子程序VUSDFLD函数，VUSDFLD函数为Abaqus中预留的子程序接口，可通过对其进行编写完成计算数据的传递。在计算初始时刻，首先读取步骤1.6中生成的三维土体强度平稳随机场对应于每个欧拉实体单元积分点的值，也即大变形有限元数值模型的三维非平稳土体强度随机场S_u0(x_c,i,y_c,i,z_c,i,θ)，随后读取初始温度场中对应于每个欧拉土体单元积分点的值Temp_i，将二者按式(15)进行叠加，生成三维土体强度非平稳随机场。如式(16)-(17)所示，根据三维土体强度非平稳随机场变量的值FV_i，更新每个欧拉单元积分点处土体的不排水抗剪强度S_u,i。

式中，

为海床表面z＝0处的土体不排水抗剪强度均值，Temp_i表示第i个欧拉土体单元的初始温度场，S_u0(x_c,i,y_c,i,z_c,i,θ)为三维土体强度平稳随机场。

FV_i(t＝0)＝S_u(x_c,i,y_c,i,z_c,i,θ) (16)

S_u,i＝S_u(FV_i(t＝0)) (17)

2.9在计算分析中，如式(18)-(19)所示，在每一个增量步中，通过状态变量SDV_i将前一增量步中的场变量值FV_i传递至当前增量步，完成三维土体强度非平稳随机场的信息传递。

FV_i(t)＝FV_i(t-dt)＝SDV_i(t) (18)

SDV_i(t)＝FV_i(t) (19)

2.10采用中心差分算法对步骤2.1-2.7中建立的大变形有限元数值模型进行求解，同时将步骤2.8-2.9中定义的子程序VUSDFLD嵌入求解器中，完成预定任务。

2.11后处理，计算完成后，通过编写Python脚本，读取计算完成后得到的结果文件中管道受坠物撞击中心位置处的凹痕变值δ并进行记录。

3、蒙特卡洛抽样模拟

3.1将步骤1.5中生成三维土体强度平稳随机场的方法重复N次，得到对应N组三维土体强度平稳随机场。

3.2将3.1步骤中生成的N组三维土体强度平稳随机场按照步骤2中的方法生成对应的N组三维土体强度非平稳随机场。

3.3将生成的N组三维土体强度非平稳随机场通过步骤2中描述的子程序技术进行大变形有限元分析，得到N组管道凹痕变形值δ_i,(i＝1,2,…,N)。

3.4对步骤3.3中得到的N组管道受冲击损伤后的凹痕变形值进行分布拟合，得到其统计特性，包括其均值μ_δ，标准差σ_δ。凹痕变形值为判断管道是否发生失效，评估管道损伤情况的重要指标，故分析所得的凹痕损伤值的统计特性可有效评估管道损伤情况。

实例计算：

首先通过物理模型试验对本发明步骤1-2中所建立的三维土体强度平稳随机场和三维土体强度非平稳随机场进行验证。该物理模型试验中，坠物为质量为m＝100kg的无杆海军锚，释放高度h＝0.4-1.0m。管道的屈服强度σ_y＝235MPa，管径D＝0.074m，壁厚t＝0.0028m，长度L＝4.9m。试验土体为粘土，管道试样埋置于土体中一定深度处，埋深为e＝0-0.2m，土体的不排水抗剪强度S_u通过扁铲测胀仪获取，其具体参数列于表1中。根据该物理模型试验的参数，采用大变形有限元数值模型进行模拟分析，得到相应工况下对应的管道凹痕值δ及相对误差绝对值err，如表2所示，总体上二者的结果吻合较好，平均相对误差绝对值约为10.8％，满足工程的设计要求，这证明了本发明建立的大变形有限元数值模型的合理性。

表1不同深度处对应的不排水抗剪强度

表2管道凹痕数值模拟及物理模型试验对比

以某工程实例为研究对象，对本发明提出的分析方法的应用进行阐述。管道的几何及材料参数见表3，坠物为质量m＝3000kg的球状坠物，其密度为ρ_o＝7850kg/m³，撞击速度为v_t＝8.873m/s；海床为粘土海床，其不排水抗剪强度的均值为μS_u＝5kPa，变异系数为cov_u＝0.3，自相关距离为η_x＝50.7m，η_y＝50.7m，η_z＝3.8m不排水抗剪强度随深度增加的梯度为k_Su＝1kPa/m。

表3某工程管道的几何参数及材料参数

图2为一次典型的三维土体强度非平稳随机场实现，可以观察到土体强度表现出明显的沿空间随机分布的特征。由于水平自相关距离η_x、η_y较大，因而土体强度在水平方向的变异性要小于竖直方向的变异性，同时考虑到土体强度随深度线性增加的特性，土体强度在竖直方向的强度将沿着直线S_u(z)＝S_u0+k_Suz变化并伴有一定幅度的波动。上述分析中展示出的特性与相关研究中非平稳随机场表现出的特性一致。

对于本工程，取蒙特卡洛法抽样的次数为N＝400次，对结果进行统计分析。作为对比分析，另进行一次定性模拟，土体强度随沿深度线性增加，不考虑其随机性，其余条件保持不变。图3给出了所有抽样工况下管道凹痕的频率直方图及累积概率分布情况，经过柯尔莫可洛夫-斯米洛夫检验，其管道凹痕的分布服从正态分布。其管道凹痕的均值为μ_δ＝85mm，标准差σ_δ＝5.57mm。管道凹痕均值μ_δ要大于确定性分析工况下的管道凹痕值δ_d＝82mm。可以看出，土体强度空间随机分布特征会降低土体的抗冲击承载力，对管道造成更大的损伤，如果不考虑该因素，得到的损伤结果偏小，会引入潜在的风险。本发明提出的分析方法可从统计意义上得到相应的特征参数，进一步可结合可靠度理论进行风险分析，对于海底管道的安全设计具有重要意义。

请参见图4，图4是本发明实施例的硬件设备工作示意图，所述硬件设备具体包括：一种考虑土体空间变异性的管道撞击损伤分析设备401、处理器402及存储设备403。

一种考虑土体空间变异性的管道撞击损伤分析设备401：所述一种考虑土体空间变异性的管道撞击损伤分析设备401实现所述一种考虑土体空间变异性的管道撞击损伤分析方法。

处理器402：所述处理器402加载并执行所述存储设备403中的指令及数据用于实现所述一种考虑土体空间变异性的管道撞击损伤分析方法。

存储设备403：所述存储设备403存储指令及数据；所述存储设备403用于实现所述一种考虑土体空间变异性的管道撞击损伤分析方法。

本发明的关键点为：

(1)实现三维土体强度非平稳随机场与大变形有限元数值模型的耦合：首先，步骤1中根据数值模型中欧拉实体单元的形心位置建立三维土体强度平稳随机场；随后，步骤2中通过设置初始温度场与步骤1中生成的三维土体强度平稳随机场叠加，得到三维土体强度非平稳随机场，并通过VUSDFLD子程序技术将随机场信息导入到完整的数值模型中，从而完成随机场与大变形有限元数值模型的耦合。

(2)形成考虑土体强度空间变异性的管道受坠物损伤分析方法，这里的变异性是指土体强度随空间发生变化。

本发明的有益效果是：

1、本发明建立了三维土体强度非平稳随机场大变形有限元数值模型

岩土工程中常需要大变形有限元技术来完成分析，针对三维大变形有限元数值模型中欧拉单元网格在空间保持固定不变的特点，本发明通过设置温度场及子程序技术，有效捕捉了土体在计算过程中的位置，实现了土体强度随空间的变化特征。

2、本发明提出了考虑土体强度三维空间变异性的管道受坠物冲击损伤分析方法

基于建立的三维土体强度非平稳随机场大变形有限元数值模型，该模型由三维土体强度非平稳随机场与大变形有限元数值模型耦合得到，结合蒙特卡洛抽样法，提出了考虑土体强度三维空间变异性的管道受坠物冲击损伤分析方法。相比于传统的确定性分析方法，该方法得到的管道损伤结果与实际情况更加接近，规避了忽略土体强度变异性而产生的风险。

以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。