一种堆石坝材料参数的空间随机场模拟方法与流程

本发明属于模拟方法领域，具体涉及一种堆石坝材料参数的空间随机场模拟方法。

背景技术：

因断面小、适应性强、施工简便、导流简化及安全性高等优点，堆石坝成为目前高坝大库的首选坝型。然而，一方面，堆石坝中的堆石材料是一种非常复杂的工程材料，其材料参数、物理力学特性随机性强、变异性大。同时，因为堆石坝材料分区较多且体型较大等原因，致使堆石坝材料参数难免存在空间差异。另一方面，由于试验测量误差、数值计算理论以及本构模型的局限无法完全模拟实际工程中材料的受力情况等，以及由测量误差所带来的统计不确定性等对计算结果的精度有很大影响。因此，在常规有限元计算中假定大范围的筑坝材料参数为唯一、确定的参数值显然无法反映工程实际，造成计算值与实测值存在较大的误差。而且，目前商业有限元软件无法实现随机有限元计算，随机场模拟与通用商业有限元软件的接口框架不够完善。

技术实现要素：

本发明的目的是提供一种堆石坝材料参数的空间随机场模拟方法，能够降低现有方法计算堆石料材料参数产生的误差。

本发明所采用的技术方案是，一种堆石坝材料参数的空间随机场模拟方法，按照以下步骤具体实施：

步骤1、利用有限元软件建立堆石坝有限元模型，对堆石坝有限元模型进行网格划分；

步骤2、选取e-b本构模型中m个材料参数对步骤1中的堆石坝有限元模型进行随机场模拟，m＝1,2,…,9，得到堆石坝有限元模型空间随机场的模拟分布图。

本发明的特点还在于，

步骤1按照以下步骤具体实施：

步骤1.1、利用有限元软件建立堆石坝有限元模型，对堆石坝有限元模型从上游到下游依次分区：上游盖重1b、上游斜铺盖1a、混凝土面板、垫层区2a、过渡区3a、主堆石区1：3bi、主堆石区2：3bii、次堆石区3c；

步骤1.2、对步骤1.1中的堆石坝有限元模型的每个分区采用八节点六面体单元进行有限元网格划分，得到若干个单元。

步骤2按照以下步骤具体实施：

步骤2.1，对步骤1垫层区2a中每个单元进行编号：ai，i＝1,2,...,n，n为单元总数；

步骤2.2，计算步骤2.1中每个单元的中心点坐标值矩阵：

其中，i＝1,2,...,n，xij(j＝1,2,...,8)为第i个单元第j个节点的横坐标值，yij(j＝1,2,...,8)为第i个单元第j个节点的纵坐标值；xi为第i个单元的中心点横坐标值，yi为第i个单元的中心点纵坐标值；

步骤2.3采用lhs方法对垫层区2a每个单元的m个材料参数进行抽样，构建垫层区2a的独立标准正态随机样本矩阵：

其中i＝1,2,...n，ξi为第i次抽样得到的矩阵，slj第l个单元的第j个材料参数的抽样结果；

步骤2.4：选取步骤2.3中m个材料参数的等效互相关系矩阵rm×m：

对上述矩阵rm×m进行cholesky分解，得到下三角矩阵l1(m×m)：

将上述得到的下三角矩阵l1进行转置得到

步骤2.5：利用步骤2.2得到的各单元中心点横纵坐标之间的相对值计算高斯型自相关函数：

其中，τx，τy分别为步骤2.2计算的各单元中心点横纵坐标之间的相对值，τx＝|xi-xj|，τy＝|yi-yj|，δh为水平波动范围，δv为垂直波动范围；

通过上述高斯型自相关函数计算垫层区2a中n个单元的自相关函数矩阵ρn×n：

其中，pij为第i个单元与第j个单元间的相关系数；

对上述矩阵ρn×n进行cholesky分解得到下三角矩阵l2(n×n)：

步骤2.6：将步骤2.3中独立标准正态随机样本矩阵ξ右乘步骤2.5得到的下三角矩阵l2，再左乘步骤2.4得到的最终得到相关标准高斯随机场样本

其中，i＝1,2,...n，，为l1的转置。

步骤2.7：对步骤2.6得到的相关标准高斯随机场样本取指数，得到垫层区2a的相关对数正态随机场矩阵hi(x,y)：

其中，为每个单元随机场参数的均值，为每个单元随机场参数的变异系数；

步骤2.8，对堆石坝有限元模型的其他分区依次重复步骤2.1～步骤2.7，进而得到每个分区的相关对数正态随机场矩阵hi(x,y)，将每个分区的相关对数正态随机场hi(x,y)的模拟值赋值给分区内相应的单元，得到本发明堆石坝有限元模型的空间随机场h(x,y)n×n的模拟分布图。

步骤2.5中高斯型自相关函数还可以是指数型自相关函数、二阶自回归型自相关函数、指数余弦型自相关函数、三角型自相关函数中的任一种。

将步骤2.8得到的随机场的模拟值采用有限元软件子程序的二次开发软件进行计算，得到堆石坝有限元模型坝体沉降、水平位移以及各主应力的分布云图，对得到的堆石坝有限元模型随机场的坝体沉降、水平位移以及各主应力的分布云图进行验证，验证的具体步骤为：

步骤a、采用有限元软件子程序的二次开发软件将步骤2得到的各单元的材料参数随机场模拟值替换步骤2中各单元的原始材料参数；

步骤b、采用有限元软件子程序的二次开发软件调用有限元软件，将步骤a中的随机场模拟值代入有限元计算，得到堆石坝有限元模型坝体沉降、水平位移以及各主应力的分布云图；

步骤c、将步骤b得到的堆石坝有限元模型坝体沉降、水平位移以及各主应力的分布云图及其数据与实测的坝体数据进行对比，得到验证结果。

步骤3中有限元软件子程序的二次开发包括以下具体步骤实施：

步骤ⅰ，设置堆石坝e-b本构模型的9个材料参数的初始取值；

步骤ⅱ，按照单元序号，将步骤2中得到的各单元材料参数随机场模拟值替换步骤1中9个原材料参数；

步骤ⅲ，计算弹性模量et和体积模量bt，判断加载或者卸载，采用加载或者卸载的弹性模量，形成切线模量矩阵，最终形成应力-应变关系矩阵；

步骤ⅳ，定义荷载施加步：前30步为分期填筑步，第30～34步为分期施加水荷载步；

步骤ⅴ，采用堆石坝的博格斯模型计算坝体后期流变；

步骤ⅵ，设置自定义输出项：水平位移、竖向沉降、最大主应力、最小主应力、应力水平，以及面板轴向与顺坡向应力；

步骤ⅶ，在dos环境下加载二次开发程序和模型数据文件直接调用有限元软件marc进行计算，重复步骤ⅰ～步骤ⅴ进行迭代计算；

步骤ⅷ，对步骤ⅶ计算的结果文件进行处理，得到应力应变等响应量的云图。

本发明的有益效果是：现有方法无法考虑堆石料材料参数在空间上的变异性，使得有限元计算应变值往往小于实测值，本发明开发了面板堆石坝非侵入式随机有限元方法，不仅使随机过程和有限元计算独立进行，可应用于不同的工程实例中，在不同的有限元软件中只需修改接口程序即可，更有利于随机有限元方法的推广和使用，而且对准确计算和把握混凝土面板堆石坝的应力应变问题具有重要的意义，有效降低了有限元计算的应变值与实测值的误差。

附图说明

图1是本发明一种堆石坝材料参数的空间随机场模拟方法的流程图；

图2是本发明有限元软件二次开发的流程图；

图3是本发明面板堆石坝有限元模型的剖面尺寸图；

图4是本发明面板堆石坝有限元模型有限元网格划分图；

图5是本发明面板堆石坝e-b本构模型参数k的随机场模拟图；

图6是本发明面板堆石坝e-b本构模型参数的随机场模拟图；

图7是本发明面板堆石坝e-b本构模型参数kb的随机场模拟图；

图8是发明所提方法计算的面板堆石坝沉降等值线图；

图9是本发明一种堆石坝材料参数的空间随机场模拟方法与传统方法的对比图。

具体实施方式

下面结合附图和具体实施方式对本发明进行详细说明。

如图1所示，是本发明一种堆石坝材料参数的空间随机场模拟方法的流程图，按照以下步骤具体实施：

步骤1、利用有限元软件建立堆石坝有限元模型，对堆石坝有限元模型进行网格划分：

步骤1.1：利用有限元软件建立堆石坝有限元模型，对堆石坝有限元模型从上游到下游依次分区：上游盖重(1b)、上游斜铺盖(1a)、混凝土面板、垫层区(2a)、过渡区(3a)、主堆石区1(3bi)、主堆石区2(3bii)、次堆石区(3c)。

步骤1.2：对步骤1.1中的堆石坝有限元模型的每个分区采用八节点六面体单元进行有限元网格划分，得到若干个单元。

步骤2，选取e-b本构模型中m个材料参数对步骤1中的堆石坝有限元模型进行随机场模拟，m＝1,2,…,9，得到堆石坝有限元模型空间随机场的模拟分布图，具体随机场模拟的过程为：

步骤2.1，对步骤1中垫层区2a分区中每个单元进行编号：ai(i＝1,2,...,n，n为单元总数)；

步骤2.2，计算步骤2.1中每个单元的中心点坐标值矩阵：

其中，i＝1,2,...,n，xij(j＝1,2,...,8)为第i个单元第j个节点的横坐标值，yij(j＝1,2,...,8)为第i个单元第j个节点的纵坐标值，xi为第i个单元的中心点横坐标值，yi为第i个单元的中心点纵坐标值。

步骤2.3采用lhs方法对垫层区2a每个单元的m个材料参数进行抽样，构建垫层区2a的独立标准正态随机样本矩阵：

其中i＝1,2,...n，ξi为第i次抽样得到的矩阵，slj为第l个单元的第j个材料参数的lhs抽样结果。

步骤2.4：选取步骤2.3中m个材料参数的等效互相关系数矩阵rm×m：

对上述矩阵rm×m进行cholesky分解，得到下三角矩阵l1(m×m)：

将上述得到的下三角矩阵l1进行转置得到

步骤2.5：利用步骤2.2得到的各单元中心点横纵坐标之间的相对值计算高斯型自相关函数：

其中，τx，τy分别为步骤2.2计算的各单元中心点横纵坐标之间的相对值，τx＝|xi-xj|，τy＝|yi-yj|，δh为水平波动范围，需根据地质统计结果确定；δv为垂直波动范围，需根据地质统计结果确定。

通过上述高斯型自相关函数计算垫层区2a中n个单元的自相关函数矩阵ρn×n：

其中，pij为第i个单元与第j个单元间的相关系数；

对上述矩阵ρn×n进行cholesky分解得到下三角矩阵l2(n×n)：

上述高斯型自相关函数还可以是指数型自相关函数、二阶自回归型自相关函数、指数余弦型自相关函数、三角型自相关函数中的任一种。

步骤2.6：将步骤2.3中独立标准正态随机样本矩阵ξ右乘步骤2.5得到的下三角矩阵l2，再左乘步骤2.4得到的最终得到相关标准高斯随机场样本

其中，i＝1,2,...n，为l1的转置。

步骤2.7：对步骤2.6得到的相关标准高斯随机场样本取指数，得到垫层区2a的相关对数正态随机场矩阵hi(x,y)：

其中，为每个单元随机场参数的均值，为每个单元随机场参数的变异系数。

步骤2.8，对堆石坝有限元模型的其他分区依次重复步骤2.1～步骤2.7，进而得到每个分区的相关对数正态随机场矩阵hi(x,y)，将每个分区的相关对数正态随机场矩阵hi(x,y)的模拟值赋值给分区内相应的单元，得到本发明堆石坝有限元模型的随机场h(x,y)n×n分布图。

将步骤2.8得到的随机场的模拟值带入有限元软件子程序的二次开发软件进行计算，得到本发明堆石坝有限元模型的坝体沉降、水平位移以及各主应力的分布云图，对得到的堆石坝有限元模型随机场的坝体沉降、水平位移以及各主应力的分布云图进行验证，验证的具体步骤为：

步骤a：采用有限元软件子程序的二次开发程序将步骤2中得到的各单元的材料参数随机场模拟值替换步骤2中各单元的原始材料参数。

步骤b：采用有限元软件子程序的二次开发程序直接调用有限元软件，将步骤a中的随机场模拟值代入有限元计算，得到堆石坝有限元模型坝体沉降、水平位移以及各主应力的分布云图。

步骤c、将步骤b得到的堆石坝有限元模型坝体沉降、水平位移以及各主应力的分布云图及其数据与实测的坝体数据进行对比，得到验证结果。

如图2所示，是本发明有限元软件子程序的二次开发程序的流程图，有限元软件子程序的二次开发由设置初始参数模块、参数赋值模块、计算步定义模块、e-b本构模型计算模块、流变模型计算模块、自定义输出模块和有限元迭代计算模块组成，计算步定义过程中考虑了分期蓄水、分期填筑以及坝体后期流变的影响。

基于有限元软件的二次开发程序，包括以下具体步骤进行开发：

步骤ⅰ、设置堆石坝e-b本构模型的9个材料参数的初始取值；

步骤ⅱ、按照单元序号，将步骤2中得到的各单元材料参数随机场模拟值替换步骤2中9个原始材料参数；

步骤ⅲ、计算弹性模量et和体积模量bt，判断加载或者卸载，采用加载或者卸载的弹性模量，形成切线模量矩阵，最终形成应力-应变关系矩阵；

步骤ⅳ、定义荷载施加步：前30步为分期填筑步，第30～34步为分期施加水荷载步；

步骤ⅴ、采用堆石坝的博格斯模型计算坝体后期流变；

步骤ⅵ、设置自定义输出项：水平位移、竖向沉降、最大主应力、最小主应力、应力水平，以及面板轴向与顺坡向应力；

步骤ⅶ、在dos环境下加载二次开发程序(.f文件)和模型数据文件(.dat文件)直接调用有限元软件marc进行计算，无需打开或修改有限元软件，以实现“非侵入式”计算，重复步骤ⅰ～步骤ⅴ进行迭代计算；

步骤ⅷ、对步骤ⅶ计算的结果文件(.t16文件)进行处理，得到应力应变等响应量的云图。

本发明一种基于堆石坝材料参数的空间随机场的模拟方法公开了堆石坝堆石材料e-b本构模型参数的空间变异性模拟方法，堆石坝有限元模型的随机模拟过程与有限元计算独立进行，可以采用不同的有限元软件去进行计算，可以应用于不同实际的工程中，增大了本发明在实际应用中的实用性。而且，本发明对有限元软件进行二次开发，可以精准计算堆石坝有限元模型的应力、应变等具体数值，减小了在堆石坝有限元模型的计算值与实测值之间的误差。本发明结合现有商业有限元软件进行了非侵入式的开发，使本发明提出的方法与有限元软件很好地结合起来，便于本发明提出的随机有限元方法的推广和应用。

实施例

步骤1：某面板堆石坝正常蓄水位2005.00m，校核洪水位2008.00m。坝顶高程2010.00m，最大坝高139.00m，坝中部布设有电磁式沉降仪测线es2；从上游到下游分区分别为：上游盖重(1b)、上游斜铺盖(1a)、混凝土面板、垫层区(2a)、过渡区(3a)、主堆石区1(3bi)、主堆石区2(3bii)、次堆石区(3c)，如图3所示。对模型进行有限元网格划分，如图4所示，共划分为1430个单元、2946个节点，其中2a区28个单元，2b区32个单元，3bi区378个单元，3bii区434个单元，3c区310个单元，其余为地基单元。

步骤2：对步骤1中的每个分区的材料参数进行随机场模拟。

步骤2.1：垫层区(2a)、过渡区(3a)、主堆石区1(3bi)、主堆石区2(3bii)、次堆石区(3c)的材料参数特征值取值如下表。

步骤2.2：以垫层区(2a)为例，随机场模拟采用与有限元计算同样的网格，计算垫层区(2a)单元的中心点坐标值矩阵为：

其中，第一列为单元序号，第二列为中心点横坐标值，第三列为中心点纵坐标值；

步骤2.3：采用lhs抽样方法对每个单元中k、kb三个材料参数进行抽样，构建每个单元的独立标准正态随机样本矩阵，每个单元的构建次数为1000次，产生的独立标准正态随机样本矩阵为：

步骤2.4：由于e-b本构模型参数之间是相互独立的，取步骤2.3中三个参数间的等效互相关系数矩阵为r：

对上述矩阵r进行cholesky分解，得到下三角矩阵l1：

将上述得到的下三角矩阵l1进行转置得到

步骤2.5：选择高斯型自相关函数得到垫层区(2a区)28个单元自相关系数矩阵ρ28×28：

对上述自相关系数矩阵ρ28×28进行cholesky分解得到另一个下三角矩阵l2：

步骤2.6：将步骤2.3得到的独立标准正态随机样本矩阵ξ左乘步骤2.5得到的下三角矩阵步骤l2，并左乘步骤2.4中得到的最终得到相关标准高斯随机场样本即：

其中，为l1的转置。

步骤2.7：对步骤2.6得到的随机场样本取指数，得到每个分区的相关对数正态随机场矩阵：

同理，对过渡区、3bi区、3bii区和3c区分别进行参数随机场的离散，得到各个分区的参数随机场矩阵的离散值，然后将各离散值赋值给分区内对应的单元，得到本发明随机场的模拟分布图，图5、图6和图7分别为坝体内三个材料参数k、kb的随机场分布图，其中图5、图6和图7中的横坐标为坝宽，纵坐标是坝高，图中每个点代表的是该单元的随机场模拟值。由于随机场单元与有限元计算采用同一套网格，因此一个单元对应一个随机场模拟值。

步骤3：采用有限元软件的二次开发程序，将各个单元的原材料参数替换为随机场模拟的参数值，图8为采用本发明方法计算的面板堆石坝沉降等值线图，曲线光滑性较差，波动性较强，出现局部闭合曲线，更加符合实际情况，图9为本发明方法与图2所示的坝体测线的实测值，以及传统方法计算的坝体沉降结果对比图，从图9可以看出：目前的常规有限元方法计算值与实测沉降值有着明显差距，尤其在坝体中部以及上部，实测沉降值明显大于计算值，这符合实际中面板堆石坝沉降值远超有限元计算值的事实；本发明所提随机有限元方法计算值与实测沉降值在坝体中部、中上部吻合度较高，在坝体下部部分测点略有误差，但整体拟合良好；且计算值与实测值比较接近，证明nsfem方法通过考虑筑坝材料的随机性和不均匀性，其计算值与实测值比较吻合。