用于弹塑性模型的粗化方法

摘要

本发明涉及一种用于弹塑性模型的粗化方法。该方法包括：获取该弹塑性模型粗化过程中所需要确定的参数及其确定方法；计算过程中首先使粗化区域达到应力平衡状态，根据所需确定的粗化模型参数数目，对所粗化区域施加应力刺激，获取对应组数的平均化的应力‑应变关系；其次根据粗化区域网格材料属性，确定粗化模型孔隙比计算公式和e‑ln(p')曲线斜率λ'、k'和孔隙比数值；根据平均化的应力‑应变关系反推获取弹塑性模型的矩阵参数m；之后可以根据所获取的参数获取粗化后的弹塑性模型，本发明可以快速准确的对非均质区域进行网格粗化，提高计算效率。

权利要求

1.一种用于弹塑性模型的粗化方法，其特征在于，包括：

根据修正剑桥弹塑性模型理论推导获取该弹塑性模型粗化过程中所需要确定的参数，该参数包括压缩过程中e-ln(p')曲线斜率λ'、反弹过程中e-ln(p')曲线斜率k'、矩阵参数m和平均孔隙比 其中，e是孔隙比，p'是有效平均主应力；

施加垂直应力和水平应力使粗化区域达到应力平衡状态；以此平衡状态作为计算基础；之后根据所需确定的粗化模型参数数目，对所粗化区域施加应力刺激，获取对应组数的应力应变关系；

对粗化区域进行应力应变平均，以获取平均化的应力和体应变εv和剪切应变εs；

根据粗化区域每个网格材料属性，确定粗化模型孔隙比计算公式；根据所述孔隙比计算公式获取e-ln(p')曲线斜率λ'和k'；

假定所粗化区域网格有相同的应力，获取粗化区域的平均化的孔隙比数值；

根据所施加的垂直和水平应力刺激结合平均化的体应变εv和剪切应变εs，结合其它所得参数，通过反推获取弹塑性模型的矩阵参数m；

根据所获取的参数获取粗化后的弹塑性模型，并利用所获取的粗化后的弹塑性模型计算粗化区域的应力应变。

2.根据权利要求1所述的粗化方法，其特征在于，所述修正剑桥弹塑性模型采用以下公式表示：式中，εv和εs分别为由dp'和dq导致的变形,e是孔隙比,η等于q/p',p'是有效平均主应力,q是偏应力，m是临界状态线的斜率，εv是体应变，εs是剪切应变，dp'是有效平均主应力变化量，dq是偏应力变化量，λ'是压缩过程中e-ln(p')曲线斜率，k'是反弹过程中e-ln(p')曲线斜率。

3.根据权利要求1所述的粗化方法，其特征在于，所述根据粗化区域每个网格材料属性，确定粗化模型孔隙比计算公式；根据所述孔隙比计算公式获取e-ln(p')曲线斜率λ'和k'的步骤中，当只有一个网格时，采用如下孔隙比计算公式：e＝a+bln(p')

式中，e为孔隙比，a和b是常数。

4.根据权利要求1所述的粗化方法，其特征在于，所述根据粗化区域每个网格材料属性，确定粗化模型孔隙比计算公式；根据所述孔隙比计算公式获取e-ln(p')曲线斜率λ'和k'的步骤中，当有多个网格时，采用如下孔隙比计算公式：式中，为所有网格的平均孔隙比，n为多个网格的数量，ai和bi为对应每种材料表征孔隙比和平均有效应力方程中的系数。

5.根据权利要求4所述的粗化方法，其特征在于，假定所粗化区域网格具有相同的应力时，采用如下孔隙比计算公式：式中，ei为第i个网格的孔隙比。

6.根据权利要求1所述的粗化方法，其特征在于，所述根据所施加的垂直和水平应力刺激结合平均化的体应变εv和剪切应变εs，获取弹塑性模型的矩阵参数m的步骤中，采用最小二乘法获取矩阵参数m。

说明书

用于弹塑性模型的粗化方法

技术领域

[0001] 本发明涉及工程力学技术领域，尤其涉及一种用于弹塑性模型的粗化方法。

背景技术

[0002] 土力学是研究土体在力的作用下的应力-应变或应力-应变-时间关系和强度的应用学科，是工程力学的一个分支。

[0003] 现有技术中，剑桥模型是研究地质作用的常用模型，如图1所示，根据正常固结粘土和弱超固结粘土的三轴试验，采用状态边界面的概念，由塑性理论的流动法则和塑性势理论，采用简单曲线配合法，建立塑性与硬化定律的函数。该模型考虑了静水压力屈服特性、压硬性、剪缩性，但破坏面有尖角，该点的塑性应变方向不易确定，其假定的弹性墙内加载仍会产生塑性变形。试验证明，如图2所示，对于正常固结粘土和弱固结的饱和重塑粘土，孔隙比e与外力p、q之间存在有唯一的。当应力状态在破坏面内时，材料表现出弹性特征。当应力路径穿过屈服面材料将展现出弹塑性特征。当材料处于临界状态，材料孔隙比不会发生变化，但是会发生偏应变。图2中，从状态A到状态D为弹塑性变形过程，代表正常固结过程；从状态C到状态B为弹性变形过程，受力体所受压力降低体积变大的反弹过程符合此规律。

[0004] 在岩土工程、盆地模拟等工程学科，广泛的存在非均质性现象，例如多层非均质、平面非均质。在数值模拟工作中岩土体被划分为很多网格，由于计算机计算水平和实际数据处理的限制，对网格数量尺寸有一定要求。根据网格划分，网格内不是均质的，而且网格数量不能过多，因此采用粗化方法即用较少的网格代替较多的网格。例如三维地质建模一般产生几十万至几百万个网格数据体，要将其用于数值模拟，需对网格进行粗化(upscaling)。目前关于弹塑性力学模型的粗化算法均采用简单平均的方法，这样会导致应变变化与原始模型不同。

发明内容

[0005] 本发明的其中一个目的在于提供一种用于弹塑性模型的粗化方法，以解决现有技术中粗化方法采用简单平均方法所引起的应变变化与原始模型不同的技术问题。

[0006] 为实现上述发明目的，本发明实施例提供了一种用于弹塑性模型的粗化方法，包括：

[0007] 根据修正剑桥弹塑性模型理论推导获取该弹塑性模型粗化过程中所需要确定的参数，该参数包括压缩过程中e-ln(p')曲线斜率λ'、反弹过程中e-ln(p')曲线斜率k'、矩阵参数m和平均孔隙比

[0008] 施加垂直应力和水平应力使粗化区域达到应力平衡状态；以此平衡状态作为计算基础；之后根据所需确定的粗化模型参数数目，对所粗化区域施加应力刺激，获取对应组数的应力应变关系；

[0009] 对粗化区域进行应力应变平均，以获取平均化的应力和体应变εv和剪切应变εs；

[0010] 根据粗化区域每个网格材料属性，确定粗化模型孔隙比计算公式；根据所述孔隙比计算公式获取e-ln(p')曲线斜率λ'和k'；

[0011] 假定所粗化区域网格有相同的应力，获取粗化区域的平均化的孔隙比数值；

[0012] 根据所施加的垂直和水平应力刺激结合平均化的体应变εv和剪切应变εs，结合其它所得参数，通过反推获取弹塑性模型的矩阵参数m；

[0013] 根据所获取的参数获取粗化后的弹塑性模型，并利用所获取的粗化后的弹塑性模型计算粗化区域的应力应变。

[0014] 可选地，所述修正剑桥模型采用以下公式表示：

[0015]

[0016] 式中，εv和εs分别为由dp'和dq导致的变形,e是孔隙比,η等于q/p',p'是有效平均主应力,q是偏应力，m是临界状态线的斜率，εv是体应变，εs是剪切应变，dp'是有效平均主应力变化量，dq是偏应力变化量，λ'是压缩过程中e-ln(p')曲线斜率，k'是反弹过程中e-ln(p')曲线斜率。

[0017] 可选地，所述根据粗化区域每个网格材料属性，确定粗化模型孔隙比计算公式；根据所述孔隙比计算公式获取e-ln(p')曲线斜率λ'和k'的步骤中，当只有一个网格时，采用如下孔隙比计算公式：

[0018] e＝a+bln(p')

[0019] 式中，e为孔隙比，a和b是常数。

[0020] 可选地，所述根据粗化区域每个网格材料属性，确定粗化模型孔隙比计算公式；根据所述孔隙比计算公式获取e-ln(p')曲线斜率λ'和k'的步骤中，当有多个网格时，采用如下孔隙比计算公式：

[0021]

[0022] 式中，为所有网格的平均孔隙比，n为多个网格的数量，ai和bi为对应每种材料表征孔隙比和平均有效应力方程的中的系数。

[0023] 可选地，假定所粗化区域网格具有相同的应力时，采用如下孔隙比计算公式：

[0024]

[0025] 式中，ei为第i个网格的孔隙比。

[0026] 可选地，所述根据所施加的垂直应力和水平应力结合体应变εv和剪切应变εs，获取弹塑性模型的矩阵参数m的步骤中，采用最小二乘法获取矩阵参数m。

[0027] 本发明可以快速准确的对非均质区域进行网格粗化，提高计算效率和满足某些情况下网格尺寸的要求。

附图说明

[0028] 通过参考附图会更加清楚的理解本发明的特征和优点，附图是示意性的而不应理解为对本发明进行任何限制，在附图中：

[0029] 图1是现有技术中剑桥模型与修正剑桥模型的示意图；

[0030] 图2是现有技术中孔隙比与有效应力关系示意图；

[0031] 图3是本发明实施例提供的弹塑性模型的粗化方法示意图；

[0032] 图4是本发明实施例提供的四种不同黏土的加压泄压演化图；

[0033] 图5是本发明实施例提供的弹塑性变形网格材料布置示意图；

[0034] 图6-1～图6-2是本发明实施例提供的粗化后dεv(左)和dεs(右)的相对误差随附加压力变化图；

[0035] 图7是本发明实施例提供的非均质通道模型示意图；

[0036] 图8是本发明实施例提供的9个非均质1×1m小网格粗化为一个3×3m的网格示意图；

[0037] 图9是本发明实施例提供的粗化后的整体区域网格模型示意图；

[0038] 图10-1～图10-2是粗化前后dεv数值变化示意图；

[0039] 图11-1～图11-2是粗化前后dεs数值变化示意图。

具体实施方式

[0040] 下面结合附图和实施例，对本发明的具体实施方式作进一步详细描述。以下实施例用于说明本发明，但不用来限制本发明的范围。

[0041] 本发明实施例提供了一种用于弹塑性模型的粗化方法，如图3所示，包括：

[0042] 根据修正剑桥弹塑性模型理论推导获取该弹塑性模型粗化过程中所需要确定的参数，参数包括压缩过程中e-ln(p')曲线斜率λ'、反弹过程中e-ln(p')曲线斜率k'、矩阵参数m和平均孔隙比

[0043] 施加垂直应力和水平应力使粗化区域达到应力平衡状态；以此平衡状态作为计算基础；之后根据所需确定的粗化模型参数数目，对所粗化区域施加应力刺激，获取对应组数的应力应变关系；

[0044] 对粗化区域进行应力应变平均，以获取平均化的应力和体应变εv和剪切应变εs；

[0045] 根据粗化区域每个网格材料属性，确定粗化模型孔隙比计算公式；根据所述孔隙比计算公式获取e-ln(p')曲线斜率λ'和k'；

[0046] 假定所粗化区域网格有相同的应力，获取粗化区域的平均化的孔隙比数值；

[0047] 根据所施加的垂直和水平应力刺激结合平均化的体应变εv和剪切应变εs，结合其它所得参数，通过反推获取弹塑性模型的矩阵参数m；

[0048] 根据所获取的参数获取粗化后的弹塑性模型，并利用所获取的粗化后的弹塑性模型计算粗化区域的应力应变。

[0049] 下面本发明实施例中的弹塑性模型采用修正剑桥模型进行说明。

[0050] 本发明一实施例中采用表1所示的4种黏土材料(St Austell kaolinite、Salins 14illite、La Bouzule clay、Marais Poitevin mud)进行网格计算。该4种黏土材料的加压泄压演化规律参见图4，具体的说明请参见文献：：Pouya,A.,Djéran-Maigre,I.,Lamoureux-Var,V.,&Grunberger,D.(1998).Mechanical behaviour of fine grained sediments:experimental compaction and three-dimensional constitutive model.Marine and Petroleum Geology,15(2),129-143。

[0051] 表1 4种黏土材料参数表

[0052]粘土类型 λ κ ν M φ0
St Austell kaolinite 0.21 0.087 0.15 1.08 27
Salins 14 illite 0.12 0.045 0.14 1.17 29
La Bouzule clay 0.21 0.084 0.15 1.19 30
Marais Poitevin mud 0.30 0.073 0.17 1.33 33

[0053] 图5是本发明实施例提供的弹塑性变形网格材料布置示意图，利用上述4种黏土材料进行布置获取修正剑桥模型。图5中，位置1布置黏土材料Salins 14 illite，位置2布置黏土材料Marais Poitevin mud，位置3布置黏土材料St Austell kaolinite，位置4布置黏土材料La Bouzule clay。为简化计算，本发明实施例中，将上述四种材料进行规则网格布置。实际应用中也可以不规则布置，只是计算稍微复杂，本发明不再介绍。

[0054] 图5所示的修正剑桥模型中，初始平均有效应力-10MPa，考虑到破坏面形状,初始垂向有效应力σ'yy＝-13.636MPa，水平有效应力σ'xx(-σ'zz)＝-8.182MPa，从而确保修正剑桥模型发生弹塑性变形。施加在修正剑桥模型的顶部边界的附加应力为dσ'yy＝0.1MPa，固定底部边界，维持左右边界原有水平有效应力并且可以自由变形。

[0055] 对上述修正剑桥模型进行粗化，该修正剑桥模型在粗化前后的参数结果参见表2与表3。

[0056] 表2 10MPa初始平均有效应力粗化参数结果

[0057]

[0058] 表3 10MPa初始平均有效应力粗化效果对比

[0059]

[0060]

[0061] 由表3可见，粗化前后的修正剑桥模型在相同应力条件下的应力应变误差变化，并且较小的相对误差可以反映本发明提供的弹塑性模型粗化方法的有效性。

[0062] 本发明另一实施例中，修正剑桥模型的初始平均有效应力-11MPa。考虑到破坏面形状，初始垂直有效应力为σ'yy＝-14.636MPa，水平有效应力σ'xx(-σ'zz)＝-9.182MPa，从而确保修正剑桥模型发生弹塑性变形。施加在修正剑桥模型的顶部边界的附加应力为dσ'yy＝0.1MPa，固定底部边界，维持左右边界原有水平有效应力并且可以自由变形。

[0063] 对上述修正剑桥模型进行粗化，该修正剑桥模型在粗化前后的参数结果参见表4与表5。

[0064] 表4 11MPa初始平均有效应力粗化参数结果

[0065]参数 原始精细模型 粗化模型 相对误差
dεv 8.0533×10-4 8.1353×10-4 1.008％
dεs 4.3581×10-4 4.3121×10-4 1.067％
垂向平均边界位移/m 2.0684×10-4 2.1397×10-4 3.332％
-5 -5
水平向平均边界位移/m 3.7999×10 4.3767×10 1.318％

[0066] 表5 11MPa初始平均有效应力粗化结果对比

[0067]平均边界应变 原始精细模型(m) 粗化模型(m) 粗化模型的相对误差
水平方向 4.8011×10-5 4.8186×10-5 0.363176％
垂直方向 2.3769×10-4 2.3579×10-4 0.805802％

[0068] 由表5可见，粗化前后的修正剑桥模型在相同应力条件下的应力应变误差变化，并且较小的相对误差可以反映本发明提供的弹塑性模型粗化方法的有效性。

[0069] 由上可知，本发明实施例提供的弹塑性模型的粗化方法所获取的粗化参数随着初始应力条件的变化而变化。为评价本发明提供的粗化方法误差随着施加应力的变化情况，在初始平均有效应力-10MPa的情况下，施加压力从力0.1MPa增加到力0.5MPa(与上面实验条件相同)，从而可得相对误差随施加压力变化情况。如图6-1～图6-2所示，修正剑桥模型粗化后的dεv随着施加应力增加而增加，并且dεs随着施加应力增加而增加。

[0070] 当粗化区域内应力条件变化较大时，参数λ'、κ'、m和e均可以设为未知，此时需要根据所设置的未知参数设置不同的应力变化以及应变结果以获取与未知参数数量相同的方程，利用上述方程求解或者拟合未知参数，因此，本发明实施例提供的粗化方法可以在应力变化范围较大时保持有效性。

[0071] 本发明又一实施例中，在修正剑桥模型中添加两条不同材料的条状带，用来模拟非均质通道系统的情况。例如，在所有网格相同的模型中插入两条不同材料的条状带，用来模拟盆地沉积中的河道和沙道的情况，如图7所示，条状带2采用Salins 14 illite黏土材料，该条状带2周围采用Marais Poitevin mud黏土材料。区域1中的每个网格大小为1m×1m，位置靠上的条状带2大小为15m×1m，位置靠下的条状带2大小为18m×1m。对该非均质通道系统的弹塑性模型进行粗化，将图7中9个1m×1m网格粗化为一个3m×3m的网格如图8所示，粗化后的非均质通道系统参见图9。

[0072] 考虑到修正剑桥模型破坏面开关，初始垂直有效应力σ'yy＝-14.636MPa，水平有效应力σ'xx(-σ'zz)＝-9.182MPa，初始平均有效应力为11MPa，从而确保修正剑桥模型发生弹塑性变形。施加在修正剑桥模型的顶部边界的附加应力为dσ'yy＝0.1MPa，固定底部边界，维持左右边界原有水平有效应力并且可以自由变形。

[0073] 对上述修正剑桥模型进行粗化，该修正剑桥模型在粗化前后的参数结果参见表4与表5。

[0074] 通过在修正剑桥模型的右侧边界施加0.01MPa的测试应力，应变数值根据本发明提供的粗化方法求解。

[0075] 之后在修正剑桥模型的右侧边界施加0.1MPa的压力增量，根据粗化前后的数据验证本发明提供的方法的有效性，

[0076] 左侧边界固定，上下边界保持初始应力σ'yy＝-14.636MPa，并且可以自由变形，粗化前后如表6所示。

[0077] 表6条状不均质模型计算结果

[0078]  剪切模量(Pa) 杨氏模量(Pa)
8 8
网格1 1.61×10 3.77×10
网格2 1.26×108 2.94×108
网格3 3.15×108 7.37×108
网格4 1.52×108 3.56×108

[0079] dεv的变化情况如图10所示，图10-1是dεv粗化前的情况，图10-2是dεv粗化后的情况，图中左侧和下部数字为网格编号，右侧为图中灰度对应数值。dεs的变化情况如图11所示。图11-1是dεs粗化前的情况，图11-2是dεs粗化后的情况，图中左侧和下部数字为网格编号，右侧为图中灰度对应数值。

[0080] 从结果可见，本发明提供的粗化方法把握了型中两条非均质通道的特征，与原始的弹塑性模型相比误差较小；而且粗化后的弹塑料性模型反映了两条通道相互干扰造成的应力应变的突变区域特点。并且粗化后的弹塑料性模型的误差相对较小，从而说明本发明提供的粗化方法有效性。

[0081] 在本发明中，术语“第一”、“第二”、“第三”仅用于描述目的，而不能理解为指示或暗示相对重要性。术语“多个”指两个或两个以上，除非另有明确的限定。

[0082] 虽然结合附图描述了本发明的实施方式，但是本领域技术人员可以在不脱离本发明的精神和范围的情况下做出各种修改和变型，这样的修改和变型均落入由所附权利要求所限定的范围之内。