一种土石坝流热耦合模型参数敏感性分析构建方法与流程

本发明属于水利工程
技术领域：
，涉及一种土石坝流热耦合模型参数敏感性分析构建方法。
背景技术：
：正常运行期的土石坝坝基与坝体温度场主要受渗流条件的影响，岩土体介质中的水分运移控制着温度场的分布，温度场的时空分布特征能反演渗流场特征。因此，研究土石坝温度场与渗流场的耦合作用，可以为通过温度场分布特征反馈坝体渗流安全性态提供理论基础。随着土石坝渗流热监测技术的发展，通过温度场来反馈渗流场及渗透破坏点已成为当前热门的工程应用技术之一。如若只分析土石坝渗流场或温度场的单场分布，而忽略两场的耦合作用与实际工程并不相符。基于此，很多学者建立了反映土石坝温度场和渗流场耦合的计算模型。然而，这些模型往往参数众多，并涉及复杂非线性动态过程。模型参数尺度多样化和制造误差等多种不确定性因素均会导致模型预测不确定性。因此，分析模型参数不确定性和敏感性，可有效研究参数对模型输出变量作用、参数与参数间相互作用对模型结果影响。通过定量研究各影响因素对模型的影响，找出决定温度场变化的主要因素，这样就能集中主要精力提高这些主要影响因素的观测精度，即使较少考虑或忽略那些次要因素，也可以得到模拟结果与实际现象接近的数值模型。目前，敏感性方法在诸多的领域都有着广泛的应用。然而，对于大坝温度场，已有的研究多是将敏感性分析方法应用于大体积混凝土浇筑的温度应力控制研究，而针对运行期受渗流影响的土石坝温度场，同样可以采用该方法探讨其主要影响因子，但相关的敏感性分析构建方法还未见报道。技术实现要素：本发明的目的是提供一种土石坝流热耦合模型参数敏感性分析构建方法，通过构建土石坝饱和-非饱和流热耦合模型，对土石坝内部渗流场与温度场进行表达，得到影响模型温度输出的主要和次要影响因子。本发明所采用的技术方案是，一种土石坝流热耦合模型参数敏感性分析构建方法，具体按照以下步骤实施：步骤1：构建土石坝饱和-非饱和流热耦合模型，包括建立流动方程以及建立土石坝热量运移模型，用于描述土石坝内部渗流场与温度场的关系；步骤2：对土石坝饱和-非饱和流热耦合模型设置边界条件，包括对渗流场设置边界条件以及温度场设置边界条件；步骤3：确定土石坝饱和-非饱和流热耦合模型的参数，选定其参数的范围；步骤4：对选定的流热耦合模型参数范围进行样本抽样；步骤5：将抽样得到的样本参数组分别代入构建的土石坝饱和-非饱和流热耦合模型中进行计算；步骤6：根据土石坝饱和-非饱和流热耦合模型的计算结果，求出土石坝模型中各监测点温度输出值的均值与方差。本发明的特征还在于，步骤1中建立的流动方程具体为土石坝饱和-非饱和瞬态渗流场，采用richards方程进行描述：式中：ρw为水的密度，cm为土体容水度，g为重力加速度，se为土体相对饱和度，ss为弹性贮水率，为拉普拉斯算子，θ为土体含水率，ks为土体饱和水力传导率，kr(θ)为坝体非饱和带土体相对水力传导率，μ(t)为水的动力粘度，为温度的函数，μ(t)＝0.00002414×10(247.8/(t+133.16))，p为压强，z为计算点位置高程，qm为渗流场源汇项，t为温度。步骤1中的土体含水率θ具体为：θ＝θr+se(θs-θr)(2)；其中，θr为残余含水率，θs为饱和含水率，α为水分特征曲线进气值的倒数，hp为压力水头，nv为水分特征曲线坡度的指示参数，通过拟合土体水分特征曲线得到，m为模型参数，m＝1-1/nv。步骤1中的土体容水度cm具体为：步骤1中的坝体非饱和带土体相对水力传导率kr(θ)具体为：步骤1中的土石坝热量运移模型采用热对流传热方程描述，具体为：式中：ρeq为等效密度，ceq为等效比热容，t为温度，t为时间，λeq为等效导热系数，cw为水的比热容，dh为水动力弥散系数，u为水的平均流速，qs为温度场源汇项。公式(6)中的水动力弥散系数dh具体为：dh＝αt|v|δij+(αl-αt)vivj|v|(7)式中，αt为横向弥散度，αl为纵向弥散度，vi为i方向流速矢量，δij为克里格常量，当i＝j时为1，否则为0；|v|为流速矢量的大小；等效密度ρeq具体为：ρeq＝(1-n)ρs+θρw+(n-θ)ρg(8)等效比热容ceq具体为：ceq＝(1-n)cs+θcw+(n-θ)cg(9)其中，cw为水的比热容，ρw为水的密度，cs为坝体土体的比热容，ρs为坝体土体的密度，cg为空气的比热容，ρg为空气的密度，n为坝体土体的孔隙率，当处于饱和状态时，θ＝n，方程(6)就是饱和含水层对流传热方程。步骤2具体为：对渗流场设置边界条件：上游坝体为定水头边界，水头为13.2m，坝面水头按其埋深呈线性变化；下游无水，下游边界为定水头边界，水头为0m；其他边界为零通量边界；对温度场设置边界条件：库水面以下边界按库水温年变化选取，为：式中，y为库水位高程，单位为m，为月份；坝体底部边界为绝热边界；库水面以上边界与大气接触，具有季节波动性，按日平均气温取值，这里采用当地气温散点资料，由下式进行拟合：式中，t0为年平均气温，l为傅里叶级数的阶数，l＝1,2…,γ，al和bl为傅里叶级数的系数，l为周期；渗流场初始条件假设为压力水头为0m，温度场假设初始温度为年平均气温。步骤3中确定土石坝流热耦合模型参数为：水力传导率ks、饱和含水率θs、残余含水率θr、坝体土体的密度ρs、坝体土体的比热容cs、坝体土体的孔隙率n和等效导热系数λeq，然后根据土石坝相关设计规范及所选取的实际工程设计资料，确定上述7个参数的范围。步骤4具体为：步骤4.1，根据步骤3选定的7个参数，设模型输出函数y＝f(x1,x2…,xk)，含k个参数，k＝7，结合各参数所服从的概率分布，将每个参数的变化范围映射到区间[0,1]，并有预设的抽样水平p＝k+1将其离散化，使每个参数只能从中取值，让k个参数都在p个抽样点上随机抽样一次，获得的向量表示为x＝(x1,x2,…,xk)；步骤4.2，从多个向量x＝(x1,x2,…,xk)中随机抽取一个向量作为输入参数“基值”向量，记为：然后从向量出发进行8次抽样，每次抽样都对中的一个或几个参数施加δ的变化量，使抽样得到的向量中每相邻2个向量都只有1个参数发生δ的变化量，并让这k个参数轮流变化一次，抽样得到的这个k+1个相邻向量被称为一条轨道。步骤5具体为：步骤5.1，在土石坝模型中设定监测温度输出值的监测点；步骤5.2，将抽样得到的参数组依次代入土石坝流热耦合模型中进行计算；步骤5.3，整理各参数组计算得到各监测点温度输出值，并列表。步骤6具体为：步骤6.1，根据整理计算得到各监测点温度输出值，计算各监测点模型计算输出的温度均值和方差，并列表；步骤6.2，将得到的均值和方差取绝对值，绘制散点图，得到影响模型温度输出的主要和次要影响因子。本发明的有益效果是：(1)本发明构建的土石坝流热耦合模型能够较真实地反映土石坝所处的流热状态，首先是考虑了坝体内部温度场与渗流场相互作用的过程；其次是考虑周期性环境温度变化的影响，属于强耦合过程。(2)土石坝内部温度场输出影响与监测点位置之间有一定的变化规律：自坝体上游至坝体下游水平方向与防渗墙前后坝顶至坝基垂直方向，参数对温度值的影响呈现增大的趋势，各参数影响模型温度值输出时与其他参数的相互作用大小也呈现增大的趋势。具体实施方式下面结合附图和具体实施方式对本发明进行详细说明。本发明一种土石坝流热耦合模型参数敏感性分析构建方法，具体按照以下步骤实施：步骤1：构建土石坝饱和-非饱和流热耦合模型，包括建立流动方程以及建立土石坝热量运移模型，用于描述土石坝内部渗流场与温度场的关系；步骤1中建立的流动方程具体为土石坝饱和—非饱和瞬态渗流场，采用richards方程进行描述：式中：ρw为水的密度，cm为土体容水度，g为重力加速度，se为土体相对饱和度，ss为弹性贮水率，为拉普拉斯算子，θ为土体含水率，ks为土体饱和水力传导率，kr(θ)为坝体非饱和带土体相对水力传导率，μ(t)为水的动力粘度，为温度的函数，μ(t)＝0.00002414×10(247.8/(t+133.16))，p为压强，z为计算点位置高程，qm为渗流场源汇项，t为温度；步骤1中的土体含水率θ具体为：θ＝θr+se(θs-θr)(2)；其中，θr为残余含水率，θs为饱和含水率，α为水分特征曲线进气值的倒数，hp为压力水头，nv为水分特征曲线坡度的指示参数，通过拟合土体水分特征曲线得到，m为模型参数，m＝1-1/nv；步骤1中的土体容水度cm具体为：步骤1中的坝体非饱和带土体相对水力传导率kr(θ)具体为：土石坝热量运移模型采用热对流传热方程描述，具体为：式中：ρeq为等效密度，ceq为等效比热容，t为温度，t为时间，λeq为等效导热系数，cw为水的比热容，dh为水动力弥散系数，u为水的平均流速，qs为温度场源汇项。公式(6)中的水动力弥散系数dh具体为：dh＝αt|v|δij+(αl-αt)vivj|v|(7)式中，αt为横向弥散度，αl为纵向弥散度，vi为i方向流速矢量，δij为克里格常量，当i＝j时为1，否则为0；|v|为流速矢量的大小；等效密度ρeq具体为：ρeq＝(1-n)ρs+θρw+(n-θ)ρg(8)等效比热容ceq具体为：ceq＝(1-n)cs+θcw+(n-θ)cg(9)其中，cw为水的比热容，ρw为水的密度，cs为坝体土体的比热容，ρs为坝体土体的密度，cg为空气的比热容，ρg为空气的密度，n为坝体土体的孔隙率，当处于饱和状态时，θ＝n，方程(6)就是饱和含水层对流传热方程；步骤2：对土石坝饱和-非饱和流热耦合模型设置边界条件，包括对渗流场设置边界条件以及温度场设置边界条件；对渗流场设置边界条件：上游坝体为定水头边界，水头为13.2m，坝面水头按其埋深呈线性变化；下游无水，下游边界为定水头边界，水头为0m；其他边界为零通量边界；对温度场设置边界条件：库水面以下边界按库水温年变化选取，为：式中，y为库水位高程，单位为m，为月份；坝体底部边界为绝热边界；库水面以上边界与大气接触，具有季节波动性，按日平均气温取值，这里采用当地气温散点资料，由下式进行拟合：式中，t0为年平均气温，l为傅里叶级数的阶数，l＝1,2…,γ，al和bl为傅里叶级数的系数，l为周期；渗流场初始条件假设为压力水头为0m，温度场假设初始温度为年平均气温；步骤3：确定土石坝饱和-非饱和流热耦合模型的参数，选定其参数的范围；由于考虑的是一个耦合问题，因此大坝温度场应该同时受到区域热状态和渗流场的共同影响。在计算区域相当大的假设条件下，假定大坝渗流体系与深部地质体之间没有热量交换，则温度边界只有上游库底水温和下游坝面温度两个温度已知边界条件，渗流介质的基本参数对温度场有一定的影响，确定土石坝流热耦合模型参数为：水力传导率ks、饱和含水率θs、残余含水率θr、坝体土体的密度ρs、坝体土体的比热容cs、坝体土体的孔隙率n和等效导热系数λeq，然后根据土石坝相关设计规范及所选取的实际工程设计资料，确定上述7个参数的范围；步骤4：对选定的流热耦合模型参数范围进行样本抽样；具体为：步骤4.1，根据步骤3选定的7个参数，设模型输出函数y＝f(x1,x2…,xk)，含k个参数，k＝7，结合各参数所服从的概率分布，将每个参数的变化范围映射到区间[0,1]，并有预设的抽样水平p＝k+1将其离散化，使每个参数只能从中取值，让k个参数都在p个抽样点上随机抽样一次，获得的向量表示为x＝(x1,x2,…,xk)；步骤4.2，从多个向量x＝(x1,x2,…,xk)中随机抽取一个向量作为输入参数“基值”向量，记为：然后从向量出发进行8次抽样，每次抽样都对中的一个或几个参数施加δ的变化量，使抽样得到的向量中每相邻2个向量都只有1个参数发生δ的变化量，并让这k个参数轮流变化一次，抽样得到的这个k+1个相邻向量被称为一条轨道；步骤5：将抽样得到的样本参数组分别代入构建的土石坝饱和-非饱和流热耦合模型中进行计算；具体为：步骤5.1，在土石坝模型中设定监测温度输出值的监测点；步骤5.2，将抽样得到的参数组依次代入土石坝流热耦合模型中进行计算；步骤5.3，整理各参数组计算得到各监测点温度输出值，并列表；步骤6：根据土石坝饱和-非饱和流热耦合模型的计算结果，求出土石坝模型中各监测点温度输出值的均值与方差；具体为：步骤6.1，根据整理计算得到各监测点温度输出值，计算各监测点模型计算输出的温度均值和方差，并列表；步骤6.2，将得到的均值和方差取绝对值，绘制散点图，得到影响模型温度输出的主要和次要影响因子。实施例一种土石坝流热耦合模型参数敏感性分析构建方法，具体按照以下步骤实施：步骤1：构建土石坝饱和-非饱和流热耦合模型，包括建立流动方程以及建立土石坝热量运移模型，用于描述土石坝内部渗流场与温度场的关系；建立的流动方程具体为土石坝饱和—非饱和瞬态渗流场，采用richards方程进行描述：式中：ρw为水的密度(kg/m3)，cm为土体容水度(1/m)，g为重力加速度(m/s2)，se为土体相对饱和度(％)，ss为弹性贮水率(1/pa)，为拉普拉斯算子，θ为土体含水率(m3·m-3)，ks为土体饱和水力传导率(m/s)，kr(θ)为坝体非饱和带土体相对水力传导率(m/s)，μ(t)为水的动力粘度(pa·s)，为温度的函数，μ(t)＝0.00002414×10(247.8/(t+133.16))，p为压强(pa)，z为计算点位置高程(m)，qm为渗流场源汇项，t为温度(℃)；土体含水率θ具体为：θ＝θr+se(θs-θr)(2)；其中，θr为残余含水率(m3·m-3)，θs为饱和含水率(m3·m-3)，α为水分特征曲线进气值的倒数(1/m)，hp为压力水头(m)，nv为水分特征曲线坡度的指示参数，通过拟合土体水分特征曲线得到，m为模型参数，m＝1-1/nv；步骤1中的土体容水度cm具体为：步骤1中的坝体非饱和带土体相对水力传导率kr(θ)具体为：土石坝热量运移模型采用热对流传热方程描述，具体为：式中：ρeq为等效密度(kg/m3)，ceq为等效比热容(j/(kg·℃))，t为温度(℃)，t为时间(s)，λeq为等效导热系数(w/(m·℃))，cw为水的比热容(j/(kg·℃))，dh为水动力弥散系数，u为水的平均流速(m/s)，qs为温度场源汇项。公式(6)中的水动力弥散系数dh具体为：dh＝αt|v|δij+(αl-αt)vivj|v|(7)式中，αt为横向弥散度(m)，αl为纵向弥散度(m)，vi为i方向流速矢量，δij为克里格常量，当i＝j时为1，否则为0；|v|为流速矢量的大小；等效密度ρeq具体为：ρeq＝(1-n)ρs+θρw+(n-θ)ρg(8)等效比热容ceq具体为：ceq＝(1-n)cs+θcw+(n-θ)cg(9)其中，cw为水的比热容(j/(kg·℃))，ρw为水的密度(kg/m3)，cs为坝体土体的比热容(j/(kg·℃))，ρs为坝体土体的密度(kg/m3)，cg为空气的比热容(j/(kg·℃))，ρg为空气的密度(kg/m3)，n为坝体土体的孔隙率(％)，当处于饱和状态时，θ＝n，方程(6)就是饱和含水层对流传热方程；步骤2：对土石坝饱和-非饱和流热耦合模型设置边界条件，包括对渗流场设置边界条件以及温度场设置边界条件；如图1所示，图1给出了土石坝坝体结构图，对渗流场设置边界条件：上游坝体h-a-b-c为定水头边界，水头为13.2m，坝面水头按其埋深呈线性变化；下游无水，下游边界d-e边界为定水头边界，水头为0m；其他边界为零通量边界；对温度场设置边界条件：库水面以下边界按库水温年变化选取，为：式中，y为库水位高程(m)，为月份，坝体底部边界e-f-g-h为绝热边界；库水面以上边界b-c-d-e与大气接触，具有季节波动性，按日平均气温取值，这里采用当地气温散点资料(资料来源于中国气象数据网)，由下式进行拟合：式中，t0为年平均气温(℃)，l为傅里叶级数的阶数，l＝1,2…,γ，al和bl为傅里叶级数的系数，l为周期，本文中周期为365天；由式(11)拟合得到该地年气温波动模型；取5阶傅立叶级数就可以满足精度要求，该地年气温回归模型系数分别为：t0＝8.8141，a1＝-14.241，a2＝-1.0891，a3＝-0.3545，a4＝-0.0169，a5＝0.0563，b1＝-2.5332，b2＝-0.3466，b3＝-0.2486，b4＝0.1595，b5＝-0.0471，采用此模型作为周期性温度边界，因此得到年平均气温为8.8℃，采用此模型可以确定各时刻的气温波动；渗流场初始条件假设为压力水头为0，温度场假设初始温度为年平均气温；步骤3：确定土石坝饱和-非饱和流热耦合模型的参数，选定其参数的范围；确定土石坝流热耦合模型参数为：水力传导率ks、饱和含水率θs、残余含水率θr、坝体土体的密度ρs、坝体土体的比热容cs、坝体土体的孔隙率n和等效导热系数λeq，然后根据土石坝相关设计规范及所选取的实际工程设计资料，确定上述7个参数的范围，如表1所示；表1土石坝流热耦合模型参数服从的概率分布参数符号单位参数服从的概率分布水力传导率ksm/sks～(10-9，10-7)饱和含水率θsm3·m-3θs～(0.36，0.49)残余含水率θrm3·m-3θr～(0.013，0.190)坝体土体的密度ρskg/m3ρs～(1390，2690)坝体土体的比热容csj/(kg·℃)cs～(400，1500)坝体土体的孔隙率n％n～(21，68)等效导热系数λeqw/(m·℃)λ～(0.25，5.03)步骤4：对选定的流热耦合模型参数范围进行样本抽样；具体为：步骤4.1，根据步骤3选定的7个参数，设模型输出函数y＝f(x1，x2…，xk)，含7个参数，结合各参数所服从的概率分布，将每个参数的变化范围映射到区间[0，1]，并有预设的抽样水平p＝k+1将其离散化，使每个参数只能从中取值，让k个参数都在p个抽样点上随机抽样一次，获得多个向量x＝(x1，x2，…，xk)；步骤4.2，从多个向量x＝(x1，x2，…，xk)中随机抽取一个向量作为输入参数“基值”向量，记为：然后从向量出发进行8次抽样，每次抽样都对中的一个或几个参数施加δ的变化量，使抽样得到的向量中每相邻2个向量都只有1个参数发生δ的变化量，并让这k个参数轮流变化一次，抽样得到的这个k+1个相邻向量被称为一条轨道；现将第1条轨道的抽样结果列表如表2，其他9条轨道不列举。表2土石坝流热耦合模型参数抽样样本样本轨道11(7.17×10-8,0.434,0.038,2504,400,34.4,2.299)2(7.17×10-8,0.434,0.038,2504,400,27.7,2.299)3(7.17×10-8,0.434,0.038,2690,400,27.7,2.299)4(5.76×10-8,0.434,0.038,2690,400,27.7,2.299)5(5.76×10-8,0.434,0.038,2690,400,27.7,2.981)6(5.76×10-8,0.434,0.038,2690,557,27.7,2.981)7(5.76×10-8,0.434,0.013,2690,557,27.7,2.981)8(5.76×10-8,0.416,0.013,2690,557,27.7,2.981)步骤5：将抽样得到的样本参数组分别代入构建的土石坝饱和-非饱和流热耦合模型中进行计算；具体为：步骤5.1，在土石坝模型中设定监测温度输出值的监测点，沿着水平方向设4个监测点，其坐标分别为1#(20,5),2#(40,5),3#(60,5),4#(80,5)；垂直方向上，在防渗墙前后分别设置3个监测点，坐标为5#(55,10),6#(55,-5),7#(55,-20),8#(60,10),9#(60,-5),10#(60,-20)，如图2所示；步骤5.2，将抽样得到的参数组依次代入土石坝流热耦合模型中进行计算；步骤5.3，整理各参数组计算得到各监测点温度输出值，并列表；如表3所列；表3土石坝流热耦合模型各监测点温度值步骤6：根据土石坝饱和-非饱和流热耦合模型的计算结果，求出土石坝模型中各监测点温度输出值的均值与方差；具体为：步骤6.1，根据整理计算得到各监测点温度输出值，计算各监测点模型计算输出的温度均值和方差，并列表，如表4所示，；表4土石坝流热耦合模型监测点的morris敏感性分析结果步骤6.2，将得到的均值和方差取绝对值，绘制散点图如图3所示，得到影响模型温度输出的主要和次要影响因子。本发明一种土石坝流热耦合模型参数敏感性分析构建方法：(1)构建的土石坝流热耦合模型能够较真实地反映土石坝所处的流热状态，首先是考虑了坝体内部温度场与渗流场相互作用的过程；其次是考虑周期性环境温度变化的影响，属于强耦合过程。(2)土石坝内部温度场输出影响与监测点位置之间有一定的变化规律：自坝体上游至坝体下游水平方向与防渗墙前后坝顶至坝基垂直方向，参数对温度值的影响呈现增大的趋势，各参数影响模型温度值输出时与其他参数的相互作用大小也呈现增大的趋势。(3)总得来看，水力传导率(ks)和坝体土体的比热容(cs)的温度均值较大，它们是对模型温度值的输出影响较大的参数；而导热系数(λeq)、饱和含水率(θs)、残余含水率(θr)、孔隙率(n)和土体密度(ρs)的μ值较小，它们是敏感性小的参数，对模型输出影响很小。其中，水力传导率(ks)对模型输出影响最大，饱和含水率(θs)、残余含水率(θr)和孔隙率(n)的温度均值为0，对模型输出影响最小，所以在热传导过程中这三个参数对坝体温度影响可以忽略。(4)水力传导率(ks)和坝体土体的比热容(cs)等参数的方差值较大，说明其在影响模型温度值输出时，与其他参数产生较大的相互作用，或者说它们对温度值输出影响是非线性的，而导热系数(λeq)、饱和含水率(θs)、残余含水率(θr)、孔隙率(n)和土体密度(ρs)的方差值相对较小，说明影响模型输出时，与其他参数相互作用较小。均值较大，相应的方差值一般也较大，说明一个参数对模型的输出值影响越大时，往往与其他参数相互作用也越大。因此，在实际研究当中，提高对水力传导率(ks)和比热容(cs)的观测精度，这样，即使较少考虑或忽略那些次要因素，也可以得到模拟结果与实际现象接近的数值模型。(5)各参数对模型温度值的影响既有正作用关系，又有负作用关系。这种正负作用关系能够为模型校正提供一定的参考。正作用关系表明在调参过程中参数的增大与减小与温度值的变化一致，而负作用关系表明参数的增大与减小与温度值的变化相反。就本研究而言，在调参时，当坝体模拟温度值小于实测值时，我们就可以通过增大正作用参数值或者减小负作用参数值来调整整个模型，从而使得模型模拟更加接近于实际情况。当前第1页12