本发明涉及深海油气管道铺设安装领域，具体涉及一种深海管道s型铺设动力模拟方法。
背景技术：
深海管道将海底油井产出的石油和天然气运送到陆地，是海洋油气开发重要的集输运手段。s型铺管法将油气管道通过铺管船上张紧器张拉着，经由托管架逐渐下放到海底，具有铺设管径大、速度快等优点，是海洋油气管道铺设最常用的方法。随着海洋油气开发迈向深海，深海复杂海洋环境使得管道铺设难度大、风险高，对管道s型铺设模拟计算方法提出更高要求。
深海s型铺设管道涉及大挠度几何变形、材料性能、边界接触和水动力荷载等非线性，给铺设动力模拟计算带来较大挑战。目前常采用悬链线法和摄动法将动力模拟简化为拟静力分析，无法准确计算海洋环境动力荷载，或采用有限差分法和非线性有限元法计算海洋环境动力荷载，但通常将托管架和海床土简化为刚性，难以准确计算管道的动力响应。
技术实现要素：
有鉴于此，本发明的目的在于提供一种深海管道s型铺设动力模拟方法，同时考虑管道与托管架以及海床土的动态接触作用，分析过程能够随时改变荷载与边界条件，更好模拟动态海洋环境作用下管道受力特征。
为实现上述目的，本发明采用如下技术方案：
一种深海管道s型铺设动力模拟方法，包括以下步骤：
步骤s1:基于向量式有限元建立s型铺设管道质点运动控制方程；
步骤s2：根据得到的s型铺设管道质点运动控制方程，采用显示中央差分法求解获得管道质点位移和转角；
步骤s3:计算铺设管道内力；
步骤s4:计算铺设管道外力；
步骤s5:更新质点位置和转角,循环步骤s2-步骤s4，直到达到计算总时长，输出动力模拟计算结果。
进一步的，所述步骤s1具体为:
步骤s11:预设在整体域坐标系(x,y,z)，管道的初始位置为水平，顶端固定在铺管船张紧器上，底端自由；
步骤s12:基于向量式有限元方法，将管道离散成n 1个由n个弯曲梁单元相连的质点，单元质量平均分于相连接的两端质点，所受内力、外力和弯矩等效施加于两端质点，管道结构的运动和受力都由质点来反映；
步骤s13:管道在重力和浮力作用下自由降落到达海床，底端锚固，形成s形态；
步骤s13:预设初始和终止时间分别为t0和tn，将整个分析历时分成一组时间点t0,…,ti-1,ti,ti 1,…,tn，时段ti-1-ti为一个途径单元，在任一途径单元内管道质点p的运动满足牛顿第二定律，运动控制方程表示为：
其中，mp为质点p的质量矩阵，ip为质点p所在截面的质量惯性矩，xp为质点p的线位移矩阵，βp为质点p的角位移向量，ζ为管道的阻尼系数，
进一步的，所述步骤s2具体为：在途径单元ti-1≤t≤ti内，一个任意连接空间点(o，p)的管道单元，已知ti-1时刻点p的位置和转角向量
中央差分法的显示积分求解控制方程的计算公式为：
其中：c1＝1/(1 ζh/2),c2＝c1(1-ζh/2)。ip是质点p的质量惯性矩向量；h为计算时间步。
进一步的，所述步骤s3具体为：
步骤s31:在主轴坐标中，由节点位置求出单元跨度li、li-1，节点转角转化为主轴坐标分量
步骤s32:计算出ti-1～ti时段内管道单元主轴方向变化向量、单元的变形,并进一步计算出管道单元的变形、扭角和弯角；
步骤s33:根据材料力学的挠曲理论，计算内力增量
步骤s34:建立静平衡方程计算出其它6个力分量，得到管道单元在虚拟位置的节点内力和节点弯矩的增量，并通过转换坐标、叠加ti-1时的内力和正向转动，得到ti时整体域坐标下的节点内力向量；
步骤s35:将节点内力提供给连接空间点做内力集成。
进一步的，所述步骤s32具体为：管道单元的变形和单元的扭角与弯角的计算公式为：
其中，li和li-1分别为单元o-p在ti和ti-1时刻的跨度；
进一步的，所述步骤s33具体为：根据材料力学的挠曲理论，内力增量计算公式为：
其中，e和g为管道的弹性模型与切线模量，a和
进一步的，所述步骤s34具体为：建立静平衡方程计算出其它6个力分量，得到管道单元在虚拟位置的节点内力和节点弯矩的增量
通过转换坐标、叠加ti-1时的内力和正向转动，得到ti时整体域坐标下的节点内力向量
其中，rt是空间转动矩阵，
进一步的，所述管道的外力包括：重力、浮力、水动力荷载、托管架滚轮反力和海床土壤抗力。
进一步的，所述步骤s4具体为：
①管道质点的重力和浮力为：
其中，mp为均分到每个单元的质量；g为重力加速度；ρw为海水的密度；vp为每个单元的体积。
②水动力荷载主要为海流对管道的拖拽力：
其中，fc,n和fc,τ分别为水流对管道质点的法向和切向拖拽力；cn和cτ分别为法向阻力系数和切向拖拽力系数；ρ为海水密度；d为管道外径；vn、vτ为海流沿管道法向和切向的相对速度；
③托管架滚轮反力：
采用多点支撑型托管架来模拟管道与托管架的相互作用。在计算滚轮反力时，首先判定每个管道单元是否与滚轮有接触：
若d0＞r1 r2，管道与滚轮分离，则fn＝0；
若d0≤r1 r2，管道与滚轮接触，则fn＝[1/(1/k1 1/k2)]×[d0-(r1 r2)]
其中，d0为管道距离滚轮的距离；r1为管道半径；r2为滚轮半径；fn为滚轮对管道的支撑力；k1和k2为管道和滚轮的接触刚度；
质点1和质点2所受的外力计算公式为：
④海床土壤抗力：
采用非线性滞回土壤模型计算土壤竖向抗力，管土的竖向作用过程分为无接触、初贯入、抬升和再贯入四种模式，土壤抗力p(z)与贯入深度z有关：
p(z)＝h(z/d)·nc(z/d)·su(z)d(10)
其中，h(z/d)为贯入系数，不同的贯入模式，贯入系数与无量纲贯入系数(z/d)的关系式不同；d是管道的外径；su(z)为海床土的不排水抗剪强度，su(z)＝su0 sug·z，su0为土壤表面泥线位置处的抗剪强度，sug为抗剪强度梯度。
本发明与现有技术相比具有以下有益效果：
1、本发明将管道离散成质点，荷载和重量集中在质点上，采用中央差分法的显示积分求解，不必构成刚度矩阵，避免迭代和不收敛的问题；
2、本发明基于途经单元概念，简化内力计算，每个时间步独立，能够较容易处理不连续行为；同时考虑管道与托管架以及海床土的动态接触作用，分析过程能够随时改变荷载与边界条件，更好模拟动态海洋环境作用下管道受力特征。
附图说明
图1是本发明一实施例中管道s型铺设模拟示意图；
图2是本发明一实施例中托管架滚轮支座反力示意图；
图3是本发明一实施例中模拟流程图。
具体实施方式
下面结合附图及实施例对本发明做进一步说明。
请参照图3，本发明提供一种深海管道s型铺设动力模拟方法，包括以下步骤：
步骤s1:基于向量式有限元建立s型铺设管道质点运动控制方程；
如图1所示，在本实施例中，在整体域坐标系(x,y,z)，管道的初始位置为水平，顶端固定在铺管船张紧器上，底端自由。基于向量式有限元方法，将管道离散成n 1个由n个弯曲梁单元相连的质点，单元质量平均分于相连接的两端质点，所受内力、外力和弯矩等效施加于两端质点，管道结构的运动和受力都由质点来反映。管道在重力和浮力作用下自由降落到达海床，底端锚固，形成“s”形态。
假定分析的初始和终止时间分别为t0和tn，将整个分析历时分成一组时间点t0,…,ti-1,ti,ti 1,…,tn，时段ti-1-ti为一个途径单元，在任一途径单元内管道质点p的运动满足牛顿第二定律，运动控制方程表示为：
其中，mp为质点p的质量矩阵，ip为质点p所在截面的质量惯性矩，xp为质点p的线位移矩阵，βp为质点p的角位移向量，ζ为管道的阻尼系数，
步骤s2：根据得到的s型铺设管道质点运动控制方程，采用显示中央差分法求解获得管道质点位移和转角；
在本实施例中，具体的，在途径单元ti-1≤t≤ti内，一个任意连接空间点(o，p)的管道单元，已知ti-1时刻点p的位置和转角向量
中央差分法的显示积分求解控制方程的计算公式为：
其中：c1＝1/(1 ζh/2),c2＝c1(1-ζh/2)。ip是质点p的质量惯性矩向量；h为计算时间步。
步骤s3:根据材料力学的挠曲理论，计算铺设管道内力；
在主轴坐标中，由节点位置求出单元跨度li、li-1，节点转角转化为主轴坐标分量
管道单元的变形和单元的扭角与弯角的计算公式为：
其中，li和li-1分别为单元o-p在ti和ti-1时刻的跨度；
根据材料力学的挠曲理论，内力增量计算公式为：
其中，e和g为管道的弹性模型与切线模量，a和
建立静平衡方程可计算出其它6个力分量，得到管道单元在虚拟位置的节点内力和节点弯矩的增量。
再通过转换坐标、叠加ti-1时的内力(矩)和正向转动，得到ti时整体域坐标下的节点内力(矩)向量。
其中，rt是空间转动矩阵，
最后，节点内力提供给连接空间点做内力集成
步骤s4:计算铺设管道外力；
在本实施例中，铺设过程管道的外力包括：重力、浮力、水动力荷载、托管架滚轮反力、海床土壤抗力。
①管道质点的重力和浮力为：
其中，mp为均分到每个单元的质量；g为重力加速度；ρw为海水的密度；vp为每个单元的体积。
②水动力荷载主要为海流对管道的拖拽力：
其中，fc,n和fc,τ分别为水流对管道质点的法向和切向拖拽力；cn和cτ分别为法向阻力系数和切向拖拽力系数；ρ为海水密度；d为管道外径；vn、vτ为海流沿管道法向和切向的相对速度。
③托管架滚轮反力：
采用多点支撑型托管架来模拟管道与托管架的相互作用，如图2所示，将托管架考虑为理想圆弧的刚性托管架，滚轮之间的相对位置不会发生改变。在计算滚轮反力时，首先判定每个管道单元是否与滚轮有接触：
若d0＞r1 r2，管道与滚轮分离，则fn＝0；
若d0≤r1 r2，管道与滚轮接触，则fn＝[1/(1/k1 1/k2)]×[d0-(r1 r2)]
其中，d0为管道距离滚轮的距离；r1为管道半径；r2为滚轮半径；fn为滚轮对管道的支撑力；k1和k2为管道和滚轮的接触刚度。
质点1和质点2所受的外力计算公式为：
④海床土壤抗力：
采用非线性滞回土壤模型计算土壤竖向抗力，管土的竖向作用过程分为无接触、初贯入、抬升和再贯入四种模式，土壤抗力p(z)与贯入深度z有关：
p(z)＝h(z/d)·nc(z/d)·su(z)d(10)
式中：h(z/d)为贯入系数，不同的贯入模式，贯入系数与无量纲贯入系数(z/d)的关系式不同；d是管道的外径；su(z)为海床土的不排水抗剪强度，su(z)＝su0 sug·z，su0为土壤表面泥线位置处的抗剪强度，sug为抗剪强度梯度。
步骤s5:更新质点位置和转角,循环步骤s2-步骤s4，直到达到计算总时长，输出动力模拟计算结果。
在本实施例中，采用matlab实现整体方案，进行深海管道s型铺设静动力模拟分析，流程如图3所示。
首先设定铺设管道长度、直径、时长等基本参数，根据步骤s1，基于向量式有限元建立s型铺设管道质点运动控制方程；根据步骤s2采用显示中央差分法求解获得管道质点位移和转角；根据步骤s3计算管道单元的内力和内力矩；再根据步骤s4计算重力、浮力、水动力等外力荷载，作用于管道质点上；最后更新质点的位置和转角，进入下一个时间步的循环，达到计算时长。动力模拟分析是在静力计算管道形成稳定“s”形态基础上，施加波浪流及船舶运动位移荷载，按照上述四个步骤进行动力循环计算，直到达到计算总时长，输出动力模拟计算结果。
以上所述仅为本发明的较佳实施例，凡依本发明申请专利范围所做的均等变化与修饰，皆应属本发明的涵盖范围。
1.一种深海管道s型铺设动力模拟方法，其特征在于，包括以下步骤：
步骤s1:基于向量式有限元建立s型铺设管道质点运动控制方程；
步骤s2：根据得到的s型铺设管道质点运动控制方程，采用显示中央差分法求解获得管道质点位移和转角；
步骤s3:计算铺设管道内力；
步骤s4:计算铺设管道外力；
步骤s5:更新质点位置和转角,循环步骤s2-步骤s4，直到达到计算总时长，输出动力模拟计算结果。
2.根据权利要求1所述的一种深海管道s型铺设动力模拟方法，其特征在于，所述步骤s1具体为:
步骤s11:预设在整体域坐标系(x,y,z)，管道的初始位置为水平，顶端固定在铺管船张紧器上，底端自由；
步骤s12:基于向量式有限元方法，将管道离散成n 1个由n个弯曲梁单元相连的质点，单元质量平均分于相连接的两端质点，所受内力、外力和弯矩等效施加于两端质点，管道结构的运动和受力都由质点来反映；
步骤s13:管道在重力和浮力作用下自由降落到达海床，底端锚固，形成s形态；
步骤s13:预设初始和终止时间分别为t0和tn，将整个分析历时分成一组时间点t0,…,ti-1,ti,ti 1,…,tn，时段ti-1-ti为一个途径单元，在任一途径单元内管道质点p的运动满足牛顿第二定律，运动控制方程表示为：
其中，mp为质点p的质量矩阵，ip为质点p所在截面的质量惯性矩，xp为质点p的线位移矩阵，βp为质点p的角位移向量，ζ为管道的阻尼系数，
3.根据权利要求1所述的一种深海管道s型铺设动力模拟方法，其特征在于，所述步骤s2具体为：在途径单元ti-1≤t≤ti内，一个任意连接空间点(o，p)的管道单元，已知ti-1时刻点p的位置和转角向量
中央差分法的显示积分求解控制方程的计算公式为：
其中：c1＝1/(1 ζh/2),c2＝c1(1-ζh/2)。ip是质点p的质量惯性矩向量；h为计算时间步。
4.根据权利要求1所述的一种深海管道s型铺设动力模拟方法，其特征在于，所述步骤s3具体为：
步骤s31:在主轴坐标中，由节点位置求出单元跨度li、li-1，节点转角转化为主轴坐标分量
步骤s32:计算出ti-1～ti时段内管道单元主轴方向变化向量、单元的变形,并进一步计算出管道单元的变形、扭角和弯角；
步骤s33:根据材料力学的挠曲理论，计算内力增量
步骤s34:建立静平衡方程计算出其它6个力分量，得到管道单元在虚拟位置的节点内力和节点弯矩的增量，并通过转换坐标、叠加ti-1时的内力和正向转动，得到ti时整体域坐标下的节点内力向量；
步骤s35:将节点内力提供给连接空间点做内力集成。
5.根据权利要求4所述的一种深海管道s型铺设动力模拟方法，其特征在于，所述步骤s32具体为：管道单元的变形和单元的扭角与弯角的计算公式为：
其中，li和li-1分别为单元o-p在ti和ti-1时刻的跨度；
6.根据权利要求4所述的所述的一种深海管道s型铺设动力模拟方法，其特征在于，所述步骤s33具体为：根据材料力学的挠曲理论，内力增量计算公式为：
其中，e和g为管道的弹性模型与切线模量，a和
7.根据权利要求4所述的所述的一种深海管道s型铺设动力模拟方法，其特征在于，所述步骤s34具体为：建立静平衡方程计算出其它6个力分量，得到管道单元在虚拟位置的节点内力和节点弯矩的增量
通过转换坐标、叠加ti-1时的内力和正向转动，得到ti时整体域坐标下的节点内力向量
其中，rt是空间转动矩阵，
8.根据权利要求1所述的所述的一种深海管道s型铺设动力模拟方法，其特征在于，所述管道的外力包括：重力、浮力、水动力荷载、托管架滚轮反力和海床土壤抗力。
9.根据权利要求8所述的所述的一种深海管道s型铺设动力模拟方法，其特征在于，所述步骤s4具体为：
①管道质点的重力和浮力为：
其中，mp为均分到每个单元的质量；g为重力加速度；ρw为海水的密度；vp为每个单元的体积。
②水动力荷载主要为海流对管道的拖拽力：
其中，fc,n和fc,τ分别为水流对管道质点的法向和切向拖拽力；cn和cτ分别为法向阻力系数和切向拖拽力系数；ρ为海水密度；d为管道外径；vn、vτ为海流沿管道法向和切向的相对速度；
③托管架滚轮反力：
采用多点支撑型托管架来模拟管道与托管架的相互作用。在计算滚轮反力时，首先判定每个管道单元是否与滚轮有接触：
若d0＞r1 r2，管道与滚轮分离，则fn＝0；
若d0≤r1 r2，管道与滚轮接触，则
其中，d0为管道距离滚轮的距离；r1为管道半径；r2为滚轮半径；fn为滚轮对管道的支撑力；k1和k2为管道和滚轮的接触刚度；
质点1和质点2所受的外力计算公式为：
④海床土壤抗力：
采用非线性滞回土壤模型计算土壤竖向抗力，管土的竖向作用过程分为无接触、初贯入、抬升和再贯入四种模式，土壤抗力p(z)与贯入深度z有关：
p(z)＝h(z/d)·nc(z/d)·su(z)d(10)
其中，h(z/d)为贯入系数，不同的贯入模式，贯入系数与无量纲贯入系数(z/d)的关系式不同；d是管道的外径；su(z)为海床土的不排水抗剪强度，su(z)＝su0 sug·z，su0为土壤表面泥线位置处的抗剪强度，sug为抗剪强度梯度。
技术总结