CN105223315A - 一种天然气水合物储层-结构物相互作用的计算预测方法 - Google Patents

Abstract

一种天然气水合物储层-结构物相互作用的计算预测方法，包括水合物参数标定模块、可视化前后处理模块以及求解模块。该方法包括：利用水合物参数标定模块确定水合物本构模型并标定本构参数；利用前后处理模块建立水合物储层-结构物有限元计算模型；利用求解模块计算在水合物不发生分解时水合物储层及结构物在自重作用下的初始应力、应变，进而结构物在受外荷载时模型应力应变；利用自主开发的大型有限元静力、动力计算软件GEODYNA结合应力释放法，计算水合物分解后水合物、土以及结构物应力、变形随时间变化规律，为水合物开采提供风险评估和预测并为施工指导。

Description

一种天然气水合物储层-结构物相互作用的计算预测方法

技术领域

本发明属于天然气水合物安全开采技术领域，涉及一种海底天然气水合物安全开采风险评估预测的分析计算预测方法。

背景技术

天然气水合物作为一种能量密度高的洁净能源，在全世界具有极大的资源量。2009年以来在我国青藏高原地区、南海北部神狐海域、广东沿海珠江口盆地东部海域相继发现了储量巨大的天然气水合物矿藏资源。这些资源的发现展示出了我国天然气水合物巨大的资源前景，有望改善我国的清洁能源结构和比重，更为重要的是能够在一定程度上解决能源依靠进口的国家能源安全问题。

天然气水合物在开采过程中会发生分解，不可避免的对天然气水合物沉积层的结构和强度产生不利影响，从而影响海底结构物的工程安全性，甚至有可能造成海底地质灾害。因此对海底天然气水合物安全开采风险进行评估显的尤为重要。然而目前我国天然气水合物资源开采安全评估关键技术研究仍处于起步阶段。突破相关技术瓶颈，形成具有我国自主知识产权的关键技术，实现天然气水合物高效、安全开采评估预测，是我国未来天然气水合物资源安全高效利用面临的重大课题及挑战。

天然气水合物储层与结构物相互作用的计算预测方法是开采安全评估中的关键技术。由于水合物沉积物的力学性质尤其是本构模型方面的研究成果比较匮乏，对于水合物开采工程中天然气水合物储层与结构物相互作用计算比较少。已有的计算多采用强度折减法，通过改变安全系数调整强度、模量等参数值，并分别进行计算，得出不同安全系数下应力、变形等，从而确定极限平衡状态下的安全系数及相应的变形，应力等。使用这种算法进行计算时需要不断调整参数，耗费大量精力，并且不能得出天然气水合物储层及结构物的应力、变形等随水合物分解的变化过程。因此需要提出一种新的算法，可以通过一套参数完成计算，得到水合物储层及结构物在水合物开采过程中的力学响应。

发明内容

本发明要解决的技术问题是：针对目前天然气水合物储层与结构物相互作用计算方法上存在的不足，提供一种了高效、准确的天然气水合物储层和结构物相互作用的计算预测方法。

本发明所采用的技术方案是：

一种天然气水合物储层-结构物相互作用的计算预测方法，该计算预测方法包含三个模块：水合物参数标定模块、可视化前处理后处理模块和计算模块，通过调用三个模块相互配合，实现预测水合物储与结构物相互作用的目的。水合物参数标定模块又细分为本构模型模块、输入输出模块、图形显示模块和优化与标定模块，它的功能主要是根据水合物三轴试验数据标定和优化非线性本构模型参数；可视化前后处理模块的功能是读入并显示有限元网格，修改单元属性，设置材料应力模型和材料参数，施加和修改约束条件，输出应力、变形；求解模块的功能是对建立好的有限元计算模型求解，得出水合物与结构物的相互作用规律。

计算中根据实际数据调用各个模块，具体步骤下：

a、在参数标定与优化模块中选取合理的本构模型，再将水合物沉积物三轴试验获取不同分解时间下的应力应变关系分别输入到水合物参数标定模块中，自动拟合得到不同分解时间下的本构参数值。根据不同分解时间本构参数值获取随时间变化的本构参数曲线，进而建立参数随时间变化的本构模型；

b、根据实际工程中水合物沉积物层、土层以及结构物范围在可视化前处理模块中建立并划分有限元网格，添加接触面，施加边界约束，设置材料应力模型并设置材料参数，其中水合物选取的是步骤a中确定的变参数本构模型，从而建立了完整的有限元计算模型；

c、利用求解模块计算在水合物不发生分解的条件下水合物储层及结构物在自重作用下的初始应力、应变，进而计算水合物物储层与结构物在受外荷载条件下的应力、应变；

d、在水合物沉积物选用随时间变化参数的本构模型基础上，在求解模块中利用应力释放法,计算在水合物分解条件下水合物储层与结构物的应力、变形等反应。

步骤a中所述的水合物本构模型目前已经实现随时间变化参数的Mohr-Coulomb模型，该模型考虑分解过程中材料属性的变化，其参数(包括密度、剪切模量、体积模量、泊松比、摩擦角和凝聚力)均定义了一条随时间变化的曲线，程序在进行计算时可根据时间步进行插值。为了更好的描述水合物沉积物的力学性质其他随时间变参数的本构模型也在开发当中，例如适用于水合物的邓肯E-B模型。

步骤b中所述的接触面，施加于结构物与水合物储层之间，施加接触面后模型计算的结果能合理的反映水合物储层与结构物之间的相对变形。接触面模型采用的是海洋地基与结构物广义塑性接触面模型，该模型可以描述接触面的剪胀与剪缩，模拟接触面的滑动、张开和闭合，模拟单调和循环荷载，从而利用一套参数就可以完成静、动力分析全过程。海洋地基与结构物广义塑性接触面模型表示为：

$D_s=D_{s0}\frac{1+e}{e}{\[{\left(\frac{{\mathrm{\σ}}_n}{p_a}\right)}^2+{\left(\frac{\mathrm{\τ}}{p_a}\right)}^2\]}^{0.5}$

$D_n=\frac{D_{n0}}{D_{s0}}{D}_s$

$d_g=r_d\mathrm{\α}\left(\right(M_c+k_m\mathrm{\ψ})\sqrt{{\mathrm{\ρ}}_{max}/\mathrm{\ρ}}-\mathrm{\η})\exp (c_0/\mathrm{\η})$

$H=H_0\frac{1}{1+\mathrm{\ψ}}\left(\frac{{\mathrm{\σ}}_n}{p_a}\right)(1-\frac{\mathrm{\ρ}}{{\mathrm{\ρ}}_p}){(1+\mathrm{\ρ})}^{-2}{f}_h$

式中，τ为剪应力，σ_n为法向应力，p_a是标准大气压，e为当前孔隙比，e₀为初始孔隙比，H为塑性模量，D_n0、D_s0为弹性参数，H₀、f_h为塑性模量参数，c₀为很小的常数取0.0001，D_n为法向刚度，D_s为剪切刚度，ψ为剪胀角，M_c为临界状态参数，α、γ_d、k_m为塑性流动方向参数，ρ为反弯点到当前应力状态的距离，ρ_p为峰值应力边界面上映射点到当前应力状态的距离，ρ_max为最大应力历史边界面上映射点到当前应力状态的距离，d_g为剪胀比。

步骤c所述的应力释放法是指同时输入分解前和分解后的材料参数，运用分解前的材料参数和应力应变关系求出各单元的初始状态(应力σ_A、应变ε_A)，然后对单元进行约束，固定其应变，再根据分解后的材料参数和应力应变关系重新计算单元应力。这样，必然会产生不平衡力。把不平衡力F加到土体上，按分解后的材料参数和应力应变关系，则可以求出分解产生的附加应力和变形。开发了应力释放荷载类来控制随时间变化的应力释放计算。这样通过特殊的本构模型和荷载组合实现了应力释放法。其中，不平衡力F用下面的式子表示：

F＝∫_vB^T[σ_A-σ_B(ε_A)]dV

其中，B为应变－位移矩阵；σ_A为分解前的应力；σ_B是与分解前应变ε_A对应的按分解后应力～应变关系计算得到的应力。

所述的水合物储层-结构物相互作用计算模型采用二维或三维模型，采用的计算原则包括，

A、边界条件设定为外边界水平方向固定，竖直方向自由，底边界X、Y、Z方向固定，上部边界水平。

B、计算模型的初始应力场、应变场由水合物储层在重力作用下的自平衡得到，将位移值置零，在此基础上计算结构物在荷载作用下与水合物作用的应力场、应变场，进而计算在水合物分解时结构物以及储层的相应应力与变形。

本发明通过建立随时间变参数的水合物本构模型，水合物三轴试验应力应变关系，标定出不同分解时间的本构参数值，利用结合应力释放法输入不同分解时间的本构参数值，求得不平衡力，再将不平衡力施加于模型，实现了水合物分解条件下水合物储层与结构物的相互作用计算。水合物储层与结构物的接触采用了海洋地基与结构物广义塑性接触面模型，该模型可以描述接触面的剪胀与剪缩，模拟接触面的滑动、张开和闭合，模拟单调和循环荷载，从而利用一套参数就可以完成静、动力分析全过程。利用本发明可以高效、准确的计算二维以及三维情况下水合物储层及结构物随水合物分解的应力和变形。

附图说明

图1是本发明的计算预测方法的模块结构图。

图2是本发明的计算预测方法的数据流程图。

图3是水本发明的计算预测方法中应力释放法示意图。

图4是本发明的计算预测方法中海洋地基广义接触面模型示意图。

具体实施方式

本实施例针对水合物特有力学性质，根据三轴试验实验数据标定出变参数水合物本构模型的参数曲线，利用自主开发的大型动力、静力有限元软件GEODYNA实现了应力释放法，从而实现了水合物储层-结构物相互作用的预测计算。方法包含三个模块，通过调用三个模块相互配合，实现预测水合物储与结构物相互作用的目的；三个模块包含：水合物参数标定模块、可视化前处理后处理模块和计算模块。

水合物参数标定模块又细分为本构模型模块、输入输出模块、图形显示模块和优化与标定模块，它的功能主要是根据水合物三轴试验数据标定和优化非线性本构模型参数；

可视化前后处理模块的功能是读入并显示有限元网格，修改单元属性，设置材料应力模型和材料参数，施加和修改约束条件，输出应力、变形；

求解模块的功能是对建立好的有限元计算模型求解，得出水合物与结构物的相互作用规律；

计算中根据实际工况调用各个模块，具体步骤下：

a、在参数标定与优化模块中选取随时间变化参数的Mohr-Columb本构模型，再将水合物沉积物三轴试验获取的分解时间为12h、24h、36h、48h下的应力应变关系分别输入到水合物参数标定模块中，模块拟合得到各分解时间下弹性模量E、粘聚力c、内摩擦角水合物密度ρ、泊松比ν、剪胀角ψ的值。根据各分解时间本构参数值获取随时间变化的本构参数曲线，各分解时间段内参数值按线性内插取值，进而建立参数随时间变化的本构模型；

b、根据实际工程中水合物沉积物层、土层以及结构物范围在可视化前处理模块中建立并划分有限元网格，添加接触面并选用广义塑性接触面模型，施加边界约束，外边界水平方向固定，竖直方向自由，底边界X、Y、Z方向固定，上部边界自由，选用Mohr-Columb本构模型并输入步骤a中获取参数曲线，至此建立了完整的有限元计算模型；

所述海洋地基与结构物广义塑性接触面模型表示为：

$D_s=D_{s0}\frac{1+e}{e}{\[{\left(\frac{{\mathrm{\σ}}_n}{p_a}\right)}^2+{\left(\frac{\mathrm{\τ}}{p_a}\right)}^2\]}^{0.5}$

$D_n=\frac{D_{n0}}{D_{s0}}{D}_s$

$d_g=r_d\mathrm{\α}\left(\right(M_c+k_m\mathrm{\ψ})\sqrt{{\mathrm{\ρ}}_{max}/\mathrm{\ρ}}-\mathrm{\η})\exp (c_0/\mathrm{\η})$

$H=H_0\frac{1}{1+\mathrm{\ψ}}\left(\frac{{\mathrm{\σ}}_n}{p_a}\right)(1-\frac{\mathrm{\ρ}}{{\mathrm{\ρ}}_p}){(1+\mathrm{\ρ})}^{-2}{f}_h$

式中，τ为剪应力，σ_n为法向应力，p_a是标准大气压，e为当前孔隙比，e₀为初始孔隙比，H为塑性模量，D_n0、D_s0为弹性参数，H₀、f_h为塑性模量参数，c₀为很小的常数取0.0001，D_n为法向刚度，D_s为剪切刚度，ψ为剪胀角，M_c为临界状态参数，α、γ_d、k_m为塑性流动方向参数，ρ为反弯点到当前应力状态的距离，ρ_p为峰值应力边界面上映射点到当前应力状态的距离，ρ_max为最大应力历史边界面上映射点到当前应力状态的距离，d_g为剪胀比。

c、利用求解模块计算在水合物不发生分解的条件下计算水合物储层及结构物在自重作用下的应力、应变，将该应力、应变作为模型的初始应力场、应变场，在此基础上计算结构物在荷载作用下与水合物作用的应力场、应变场。

d、将步骤c中获得的应力、应变场输入到计算模块中，将位移值置零，利用随时间变化参数的本构模型与应力释放法,计算在水合物分解条件下水合物储层与结构物的应力、变形等反应。利用可视化后处理模块处理结果数据研究水合物储层-结构物相互作用的规律。

所述应力释放法是指同时输入分解前和分解后的材料参数，运用分解前的材料参数和应力应变关系求出各单元的初始状态(应力σ_A、应变ε_A)，然后对单元进行约束，固定其应变，再根据分解后的材料参数和应力应变关系重新计算单元应力。这样，必然会产生不平衡力。把不平衡力F加到土体上，按分解后的材料参数和应力应变关系，则可以求出分解产生的附加应力和变形。开发了应力释放荷载类来控制随时间变化的应力释放计算。这样通过特殊的本构模型和荷载组合实现了应力释放法。其中，不平衡力F用下面的式子表示：

F＝∫_VB^T[σ_A-σ_B(ε_A)]dV

其中，B为应变－位移矩阵；σ_A为分解前的应力；σ_B是与分解前应变ε_A对应的按分解后应力～应变关系计算得到的应力。

案例一

为水合物储层与桩基础相互作用，一地基分上、中、下三层土，水合物层厚度30m，上部和下部粘土层厚度也均为30m。一根55m长混凝土桩竖直埋入地基正中央，入土深度为50m，桩直径0.5m，桩顶承受100t荷载。上覆和下卧土层参数均取为ρ＝1.98g/cm³,c＝30kPa。水合物分解时间48h，利用参数标定模块拟合并假定参数值按分解时间线性变化，其中弹性模量E由186MPa变为18.6MPa，粘聚力c为0，内摩擦角φ由39.4°变为15°，重度ρ由19.9KN/m³变为19.2KN/m³。本例采用2维有限元模型，利用以上资料建立有限元计算模型，在前处理模块中施加约束，水合物选定随时间变参数的Mohr-Columb模型，而粘土选择Mohr-Columb模型。首先计算模型在自重下的初始应力、应变，而后将结果输入到计算模块，施加桩顶荷载计算模型的应力，应变，最后再将应力、应变结果输入到计算模块中，计算水合物分解条件下，模型的力学响应。计算得出随水合物分解，水平位移增加不大，地基和混凝土桩主要表现为竖向沉降，竖向最大位移为0.51m。

案例二

为荔湾3-1气田工程某海洋平台的打入桩基础，该基础位于海底1200m处，所在区域土层含有水合物层，水合物埋深100m，厚25m。采用三维有限元模型进行分析，模型两个水平方向均取为300m，竖直方向取200m。打入桩长度95m，桩径2.438m，入土深度90m，在距桩顶10m处施加12000KN荷载，与水平方向呈5度角。各土层选用Mohr-Columb本构模型，水合物选用随时间变参数Mohr-Columb本构模型。基础所在区域土层分布情况及相应本构参数值由地质勘探及室内实验给出。通过室内三轴试验获得分解0h、12h、24h、36h、48h的水合物应力应变曲线，将各分解时间下应力应变关系输入到参数标定与优化模块，程序将拟合出各分解时间下的最优本构参数值，即可确定随时间变化本构参数曲线，各时间段内参数值按分解时间线性内插得到。其中水合物内摩擦角由38.4°变为18.6°，粘聚力c为0，弹性模量E由7.48×10⁸Pa变为2.84×10⁸Pa，剪胀角ψ为0°，泊松比ν为0.35，有效重度γ′由9.8KN/m³变为9.2KN/m³。打入桩为钢桩采用线弹性模型E＝2×10¹¹Pa，泊松比ν为0.3。首先根据实际工程土层资料，建立三维有限元几何模型并划分单元网格，然后在前处理模块中选定土层和水合物层的本构及相应参数值，再添加接触面，接触面使用广义塑性接触模型来模拟海洋地基与桩基础的相对变形，最后施加约束，模型侧面边界约束水平方向，底面边界约束x、y、z方向，顶面边界自由，由此建立了完整的有限元计算模型。在计算模块中首先计算模型在自重及水压力作用下的应力、应变，然后将结果作为初始应力、应变输入到计算模块计算打入桩在外荷载作用下的应力、应变，然后将结果文件输入计算模块中，在此基础上计算水合物分解时水合物层、土层以及打入桩的应力、应变，最后将计算结果输入到可视化后处理模块中，分析水合物分解时模型的响应规律。由模拟结果得出：水合物分解桩基础及水合物以上土层发生整体沉降而水平方向位移较小，随着水合物分解沉降值增大，分解48h时打入桩沉降值为0.32m。

Claims (6)

1.一种天然气水合物储层-结构物相互作用的计算预测方法，其特征在于：

该计算预测方法包含三个模块：水合物参数标定模块、可视化前处理后处理模块和计算模块；

水合物参数标定模块分为本构模型模块、输入输出模块、图形显示模块和优化与标定模块，其根据水合物三轴试验数据标定和优化非线性本构模型参数；

可视化前后处理模块：读入并显示有限元网格，修改单元属性，设置材料应力模型和材料参数，施加和修改约束条件，输出应力、变形；

求解模块：对建立好的有限元计算模型求解，得出水合物与结构物的相互作用规律；

具体步骤如下：

a、在参数标定与优化模块中选取合理的本构模型，再将水合物沉积物三轴试验获取不同分解时间下的应力应变关系分别输入到水合物参数标定模块中，自动拟合得到不同分解时间下的本构参数值；根据不同分解时间本构参数值获取随时间变化的本构参数曲线，进而建立参数随时间变化的本构模型；

b、根据实际工程中水合物沉积物层、土层以及结构物范围在可视化前处理模块中建立并划分有限元网格，添加接触面，施加边界约束，设置材料应力模型并设置材料参数，其中水合物选取的是步骤a中确定的变参数本构模型，从而建立了完整的有限元计算模型；

c、利用求解模块计算在水合物不发生分解的条件下水合物储层及结构物在自重作用下的初始应力、应变，进而计算水合物物储层与结构物在受外荷载条件下的应力、应变；

d、在水合物沉积物选用随时间变化参数的本构模型基础上，在求解模块中利用应力释放法,计算在水合物分解条件下水合物储层与结构物的应力、变形。

2.根据权利要求1所述的一种天然气储层-结构物相互作用的计算预测方法，其特征在于：步骤a中所述的水合物本构模型随时间变化参数的Mohr-Coulomb模型，该模型考虑分解过程中材料属性的变化，其参数均定义了一条随时间变化的曲线，进行计算时根据时间步进行插值。

3.根据权利要求2所述的一种天然气储层-结构物相互作用的计算预测方法，其特征在于：步骤b中所述的接触面，施加于结构物与水合物储层之间，施加接触面后模型计算的结果能合理的反映水合物储层与结构物之间的相对变形；接触面模型采用的是海洋地基与结构物广义塑性接触面模型，该模型描述接触面的剪胀与剪缩，模拟接触面的滑动、张开和闭合，模拟单调和循环荷载，从而利用一套参数完成静、动力分析全过程；海洋地基与结构物广义塑性接触面模型表示为：

式中，τ为剪应力，σ_n为法向应力，p_a是标准大气压，e为当前孔隙比，e₀为初始孔隙比，H为塑性模量，D_n0、D_s0为弹性参数，H₀、f_h为塑性模量参数，c₀为很小的常数取0.0001，D_n为法向刚度，D_s为剪切刚度，ψ为剪胀角，M_c为临界状态参数，α、γ_d、k_m为塑性流动方向参数，ρ为反弯点到当前应力状态的距离，ρ_p为峰值应力边界面上映射点到当前应力状态的距离，ρ_max为最大应力历史边界面上映射点到当前应力状态的距离，d_g为剪胀比。

4.根据权利要求3所述的一种天然气储层-结构物相互作用的计算预测方法，其特征在于：步骤c所述的应力释放法是指同时输入分解前和分解后的材料参数，运用分解前的材料参数和应力应变关系求出各单元的初始状态，然后对单元进行约束，固定其应变，再根据分解后的材料参数和应力应变关系重新计算单元应力；把产生的不平衡力F加到土体上，按分解后的材料参数和应力应变关系，求出分解产生的附加应力和变形；不平衡力F用下面的式子表示：

F＝∫_VB^T[σ_A-σ_B(ε_A)]dV

其中，B为应变－位移矩阵；σ_A为分解前的应力；σ_B是与分解前应变ε_A对应的按分解后应力～应变关系计算得到的应力。

5.根据权利要求1-4任一所述的一种天然气储层-结构物相互作用的计算预测方法，其特征在于，所述的水合物储层-结构物相互作用计算模型采用二维或三维模型。

6.根据权利要求5所述的一种天然气储层-结构物相互作用的计算预测方法，其特征在于，所述的水合物储层-结构物相互作用计算模型采用的计算原则包括，

A、边界条件设定为外边界水平方向固定，竖直方向自由，底边界X、Y、Z方向固定，上部边界水平；

B、计算模型的初始应力场、应变场由水合物储层在重力作用下的自平衡得到，将位移值置零，在此基础上计算结构物在荷载作用下与水合物作用的应力场、应变场，进而计算在水合物分解时结构物以及储层的相应应力与变形。