本发明属于水利工程技术领域，具体涉及一种基于comsol数学模块的河岸潜流带水热耦合建模方法及系统。背景技术：河岸潜流带是指河岸带内地表水与地下水之间相互作用的过渡区域，它是河流中物质能量循环和动植物栖息的重要场所，被誉为河流的“肝脏”。目前，河岸潜流带已成为多个学科研究的热点。河岸潜流带与河流地表水交换过程时刻伴随着能量的传递和交换，因此，温度作为能量的直观载体，是能够反映交换过程时空变化的表征因子，其易于观测、无污染，作为天然示踪剂具有一定的优势。近年来，随着温度场与渗流场耦合机理研究的深入，对河岸潜流带的研究已经从传统的水文学及水动力学方法发展到温度示踪法。虽然关于利用室内或者野外温度示踪试验来研究河流相关问题很多，但是试验成本比较高，且缺乏有效的水热耦合模型来分析河岸潜流带水热动态变化。土壤导热系数作为水热耦合模型中的关键参数，直接影响着模型模拟的精度。然而，在现有的耦合模型或者建模软件中都是假定土壤的导热系数为固定值。理论上来说这是不准确的，因为土壤的导热系数与土壤类型、粒径分布、土壤结构、孔隙度、饱和度和温度等因素有关，其中孔隙度和饱和度对它的影响最大。这对于河岸非饱和带的土体水热耦合研究而言显得尤为重要。尽管关于预测导热系数的经验模型已被广泛的提出，但是对于将导热系数经验模型考虑到水热耦合模型中的研究却鲜有报道。技术实现要素：本发明的目的在于克服现有技术中的不足，提出了一种基于comsol数学模块的河岸潜流带水热耦合建模方法，解决了现有技术中在利用comsol软件对河岸潜流带水热运移过程进行模拟时不能够将导热系数经验模型考虑到水热耦合模型中的技术问题。为解决上述技术问题，本发明提供了一种基于comsol数学模块的河岸潜流带水热耦合建模方法，其特征是，包括以下过程：在comsol软件中构建河岸潜流带水热耦合模型，所述河岸潜流带水热耦合模型包括水流运动方程、热量运移方程和导热系数经验模型；利用comsol数学模块中的系数型偏微分方程对热量运移方程进行等效；导入河岸潜流带几何模型，并定义河岸潜流带土体水力参数、热性质参数和导热系数经验模型；对河岸潜流带水热耦合模型进行边界条件设置；对河岸潜流带几何模型进行网格剖分；计算得到河岸潜流带温度场。进一步的，所述水流运动方程为：河岸潜流带水流运动方程采用richard方程进行描述：式中：ρw为水的密度，kg/m3；cm为容水度，m-1；g为重力加速度，m/s2；se为土体的相对饱和度；ss为弹性贮水率，pa-1；p为压强，pa；t为时间，s；为拉普拉斯算子；θ为体积含水率，m3/m3；ks为介质饱和渗透系数，m/s；kr(θ)为非饱和带相对渗透系数，m/s，是体积含水率θ的函数；μ为水的动力粘度，pa·s，μ(t)＝0.00002424×10[247.8/(t+133.16)]；t为温度；z为计算点高程，m；qm为水流源；土壤水力函数选用vangenuchten模型进行描述：θ＝θr+se(θs-θr)(2)式中：θr为残余含水率，m3/m3；θs为饱和含水率，m3/m3；hp为压力水头，m；α为水分特征曲线进气值得倒数，m-1；β为水分特征曲线坡度的指示参数，通过拟合土壤水分特征曲线得到；m＝1-1/β。进一步的，所述热量运移方程为：河岸潜流带热量运移方程可用热对流传热方程表示：式中：cw为水的比热容，j/kg·℃；cs为土体的比热容，j/kg·℃；t为温度，℃；keff为有效导热系数，w/m3·℃；dh为水动力弥散系数，m2/s；u为平均流速，m/s，在数值上u＝v/θ，v为darcy渗流流速。进一步的，所述导热系数经验模型为：导热系数经验模型可选用chung&horton模型或campbell模型，模型公式分别为：keff＝b1+b2θ+b3θ0.5(7)式中：b1、b2和b3分别为0.228，-2.406和4.909；keff＝a+bθ-(a-d)exp[-(cθ)e](8)式中：mclay为黏土含量，ρb为堆积密度(g/cm3)。进一步的，所述利用comsol数学模块中的系数型偏微分方程对热量运移方程进行等效包括以下过程：comsol中提供的系数型偏微分方程如下：式中：ea、da、c、α、γ、β、a、q和h均为自定义系数，ω为求解域，是ω的外边界，n是的外法线方向；f、g、r分别为求解域和边界上的源项；利用系数偏微分方程的形式来等效热量运移方程，等效过程如下：在热量运移方程中，可构成方程的时变项，可设定为对流传导弥散项，f则可作为热量源汇项；即设定ea＝0，da＝θcw+(1-n)cs，c＝keff+θcwdh，α＝0，γ＝0，β＝θcwu，a＝0，f＝0，此时pde方程即等效于热量运移方程；热量运移方程通过土壤含水率与richards方程耦合，又通过导热系数keff实现导热系数经验模型与热量运移方程之间的耦合，即可实现利用偏微分方程对热量运移方程进行等效。进一步的，所述对河岸潜流带水热耦合模型进行边界条件设置包括：对于渗流场，河岸潜流带左右边界设置为无流动边界；底部边界设置为透水层边界；土体与河水交界面设置为变水头边界，压力水头初始值设定为0m；对于温度场，河岸潜流带左、右和底部边界设置为绝热边界；与大气接触边界设置为大气边界，与河水接触边界设置为水温边界，初始温度场根据实际工程中所测的初始均值进行设置。与现有技术相比，本发明所达到的有益效果是：实现了将不同导热系数经验模型耦合到水热耦合模型当中，并能够对不同导热系数经验模型下的水热耦合模型模拟效果进行比较和评价，解决了现有技术中存在的不能够将导热系数经验模型耦合到水热耦合模型中，且无法对多种导热系数经验模型在实际工程中的应用性能表现进行比较分析的问题。本发明基于comsol数学模块的河岸潜流带水热耦合建模方法可以将多种导热系数经验模型耦合到水热耦合模型当中，从而得到最优的方案，在一定程度上可以提高河岸潜流带水热耦合模型的模拟精度。附图说明图1是comsol求解河岸潜流带水热耦合模型流程图；图2是实测河道水温在研究时段内的变化图；图3是实测环境温度在研究时段内的变化图；图4是实测水头在研究时段内的变化图；图5是chung&horton导热系数经验模型表达式写入图；图6是河岸潜流带温度测点布置示意图；图7是chung&horton模型和campbell模型下的河岸潜流带水热耦合模型模拟温度与观测温度对比时序曲线图：(a)为0.95m温度测点处温度实测值与chung&horton模型和campbell模型结果对比；(b)为1.40m温度测点处温度实测值与chung&horton模型和campbell模型结果对比；(c)为2.00m温度测点处温度实测值与chung&horton模型和campbell模型结果对比具体实施方式下面结合附图对本发明作进一步描述。以下实施例仅用于更加清楚地说明本发明的技术方案，而不能以此来限制本发明的保护范围。本发明的一种基于comsol数学模块的河岸潜流带水热耦合建模方法，参见图1所示，具体包括以下过程：步骤1，在comsol软件中构建河岸潜流带水热耦合模型，所述河岸潜流带水热耦合模型包括水流运动方程、热量运移方程和导热系数经验模型。comsol软件是现有技术中已知的多物理场仿真软件，在comsol软件中构建河岸潜流带水热耦合模型，所述河岸潜流带水热耦合模型包括水流运动方程、热量运移方程和导热系数经验模型。水流运动方程计算的是河岸潜流带的渗流场，热量运移方程是为了求解河岸潜流带温度场。模型的构建过程如下：1)水流运动方程河岸潜流带水流运动方程采用richards方程进行描述：式中：ρw为水的密度，kg/m3；cm为容水度，m-1；g为重力加速度，m/s2；se为土体的相对饱和度；ss为弹性贮水率，pa-1；p为压强，pa；t为时间，s；为拉普拉斯算子；θ为体积含水率(水的体积/总体积)，m3/m3；ks为介质饱和渗透系数，m/s；kr(θ)为非饱和带相对渗透系数，m/s，是体积含水率θ的函数；μ为水的动力粘度，pa·s，μ(t)＝0.00002424×10[247.8/(t+133.16)]；t为温度；z为计算点高程，m；qm为水流源汇项。土壤水力函数选用vangenuchten模型进行描述：θ＝θr+se(θs-θr)(2)式中：θr为残余含水率，m3/m3；θs为饱和含水率，m3/m3；hp为压力水头，m；α为水分特征曲线进气值得倒数，m-1；β为水分特征曲线坡度的指示参数，通过拟合土壤水分特征曲线得到；m＝1-1/β。2)热量运移方程河岸潜流带热量运移方程可用热对流传热方程表示：式中：cw为水的比热容，j/kg·℃；cs为土体的比热容，j/kg·℃；t为温度，℃；keff为有效导热系数，w/m3·℃；dh为水动力弥散系数，m2/s；u为平均流速，m/s，在数值上u＝v/θ，v为darcy渗流流速。公式6中，等式左边表示温度在变饱和条件下随时间的变化，即为时变项；等式右边第一项表示热传导项，第二项表示热弥散项，第三项表示热对流项，第四项(qs)表示热量源汇项。3)导热系数经验模型导热系数经验模型可选用chung&horton模型或campbell模型，模型公式分别为：keff＝b1+b2θ+b3θ0.5(7)式中：b1、b2和b3分别为0.228，-2.406和4.909。keff＝a+bθ-(a-d)exp[-(cθ)e](8)式中：mclay为黏土含量，ρb为堆积密度(g/cm3)。该模型构建方法可以灵活的更换导热系数经验模型，并且可以形成不同导热系数经验模型下水热耦合模型计算结果的对比，这一点是传统模型构建方法所不能实现的。步骤2，利用comsol数学模块中的系数型偏微分方程对热量运移方程进行等效；comsol中自带计算传热模块是用热量运移方程进行计算的，但是不能修改导热系数经验模型，这里用pde模块来代替comsol自带的传热模型进行计算，pde是偏微分方程，相当于人为编辑公式计算，可以灵活修改导热系数经验模型，因此必须要进行等效。comsol中提供的系数型pde如下：式中：ea、da、c、α、γ、β、a、q和h均为自定义系数，ω为求解域，是ω的外边界，n是的外法线方向；f、g、r分别为求解域和边界上的源项。利用系数偏微分方程(pde)的形式来等效热量运移方程(6)，等效过程如下：在热量运移方程中，可构成方程的时变项，可设定为对流传导弥散项，f则可作为热量源汇项。即设定ea＝0，da＝θcw+(1-n)cs，c＝keff+θcwdh，α＝0，γ＝0，β＝θcwu，a＝0，f＝0，此时pde方程即等效于热量运移方程。热量运移方程通过土壤含水率(θ)与richards方程耦合，又通过导热系数(keff)实现导热系数经验模型与热量运移方程之间的耦合，即可实现利用偏微分方程对热量运移方程进行等效。(热量运移方程和richards方程中均包含了土壤含水率(θ)，计算求解时相当于利用土壤含水率(θ)将两个式子相关联，即构成耦合，同理导热系数经验模型和热量运移方程也是如此。)步骤3，导入河岸潜流带几何模型，并定义河岸潜流带土体水力参数、热性质参数和导热系数经验模型；前面河岸潜流带模型指的是数学模型，即式1～9，几何模型模型根据实际工程断面数据绘制，仅仅是字面上的几何尺寸模型。在水利工程中，水力参数指的是求解水流方程的参数(即式1-5中未知参数，具体包括ks、θs、θr、ρs、α、β)，热性质参数指的是求解热流运移方程所需参数(即式6中未知参数n、cs、cw；另外，导热系数实际是式6中的参数，只不过本专利采用经验公式计算，按照归类，求解导热系数的参数也应该属于热性质参数，因此mclay、ρb和ρs也作为热性质参数。这些参数均列于表1中)。取实际工程的河岸潜流带土体水力参数和热性质参数，如表1所示。图2、图3和图4分别给出了水温、水头和环境温度在研究时段的变化。利用comsol对河岸潜流带水热运移特性进行模拟，如果模型能够于实际吻合，后期可免除实地测量，而直接通过模型进行预测，避免费时费力、耗费成本和资源。从而保证数值计算结果能够对实际工程具有参考价值。表1河岸潜流带土体水力参数和热性质参数步骤3.1，从材料库中添加流体材料(流体和气体>>liquids>>water)和土质材料(材料库>>minernals,rocksandsoils>>soil)，并在“材料属性明细”中设定不同河岸潜流带土体的水力参数和热性质参数，如密度ρs和渗透系数ks；步骤3.2，对richards方程中计算所需参数进行设定，如饱和含水率θs、残余含水率θr和vangenuchten模型参数α和β；步骤3.3，在“组件”的“定义”下面以三次样条插值的形式(对于图2、3、4中这种离散的数据用此方法导入精度较高)将图2水温、图3环境温度及图4水头变化数据导入，以备后期定义边界条件所需；步骤3.4，在“组件”的“定义”下面添加“变量”，通过表达式的形式将导热系数经验模型写入，并定义模型公式的其余未知参数(n、mclay、ρb、cw、cs、dh)；图5以chung&horton模型为例，展示了具体的导热系数经验模型的表达式写入形式，图中①表示在河岸潜流带中设置的三个温度监测点；②为chung&horton导热系数经验模型表达式，其中c_w对应于水的比热容cw，c_s对应于土体的比热容cs，k_eff对应于有效导热系数keff，d_hx对应于水动力弥散系数dhx方向的分量，d_hy对应于水动力弥散系数dhy方向的分量，n对应于孔隙率n；③为图2、图3、图4中数据导入软件中的结果；④为材料组件；⑤为richards方程；⑥为等效热量运移方程后的系数偏微分方程；⑦为网格剖分组件；步骤3.5，在“组件”的“定义”下面添加“域点探针”，并设定域点探针监测的位置和监测的对象。因实际试验中在0.95m1.40m和2.00m处布置温度传感器，而在仿真过程中，设置域点探针相当于布设传感器，监测对象为温度。步骤4，对河岸潜流带水热耦合模型进行边界条件设置；步骤4.1，对于渗流场，河岸潜流带左右边界(如图6所示，af和bc)设置为无流动边界；底部边界(ab)设置为透水层边界；土体与河水交界面设置为变水头边界(因为水位是变化的，所以这个边没有办法具有定下来，只能描述土体与河水交接面，如图6中de边，de是一个大致的范围，它为变水头边界，位置是不断的在变，但是大致就在这个范围，变水头设置不能直接设定为一个常数，而是变化的离散数据，因此需要调用步骤3.3中的水头数据)，压力水头初始值设定为0m。图6是根据河岸带断面数据所画出；此图是整个河道剖面图，传感器布设在河岸带处，只是个示意图；边界条件是根据河岸潜流带所处的实际情况进行设置的；区域是根据土质渗透系数不同进行划分的，因为不同位置及深度，土壤的渗透系数是不一样的；图中温度测点是实际野外试验中，在此三个位置处布置的温度传感器，用于测量河岸潜流带不同位置处的温度变化。图6中⑧为0.95m、1.40m和2.00温度监测点；⑨为区域1；⑩为区域2；为区域3，x坐标轴代表河岸潜流带跨度方向，x＝0取在距离河床中心7m处，y坐标轴代表河岸潜流带深度方向，并取河床底部高程为y＝0，整个河岸潜流带几何模型是依据实测断面数据通过cad绘制出来；a～f为边界编号；步骤4.2，对于温度场，河岸潜流带左、右和底部边界(如图6所示，af、bc和ab)设置为绝热边界；与大气接触边界设置为大气边界(cd和ef，通过调用步骤3.3中的环境温度)，与河水接触边界设置为水温边界(de，通过调用步骤3.3中的水温)。初始温度场根据实际工程中温度传感器所测的初始均值进行设置。步骤5，在“网格”中设定网格的剖分形式以及网格单元尺寸；利用comsol自带的自动剖分三角形网格对河岸潜流带几何模型(即图5中a～f的闭合区域)进行网格剖分；在网格单元尺寸设定中，将网格单元尺寸校准为流体动力学，并采用预定义较细化尺寸。步骤6，设定求解器为瞬态过程(求解分为瞬态和稳态两种，根据需要设置为瞬态)，并设定求解步长为1h和求解时间为1000h，运用隐式欧拉向后差分法对水流运动方程和热量运移方程中的时间项进行离散，采用非线性迭代修正阻尼牛顿法求解方程。求解上述构建的方程的最终目的是要求温度(t)，即pde方程中的u，求解的目的是看模拟结果与实测结果之间是否一致。模型验证及对比：本发明采用均方根误差(rmse)、决定系数(r2)、平均绝对误差(mae)和平均相对误差(mre)对不同导热系数经验模型下的河岸潜流带水热耦合模型的模拟精度进行评估：式中：n为样本总量，oi和si分别为观测和模拟温度值，为平均实观测值，rmse值越小，模拟值与观测值之间的偏差越小，模拟值与观测值的一致性越好，模型的模拟结果越准确可靠；r2在-∞～1之间变化，其值越大表明模拟结果越好，一般r2>0.6认为模拟值和观测值之间一致性比较好；mae和mre的值越小表明结果越好。根据实际工程测得的河岸潜流带温度变化数据，结合上述模型构建方法，图7给出了chung&horton和campbell导热系数经验模型下河岸潜流带水热耦合模型的模拟与实测温度随时间的变化曲线，图7(a)为0.95m温度测点处温度实测值与chung&horton模型和campbell模型结果对比；图7(b)为1.40m温度测点处温度实测值与chung&horton模型和campbell模型结果对比；图7(c)为2.00m温度测点处温度实测值与chung&horton模型和campbell模型结果对比。各模型的性能表现通过表2给出。表2模型性能评价指标结果从图7和表2可以看出，无论从哪个评价指标来看，基于chung&horton模型的河岸潜流带水热耦合模型具有较高的模拟精度，各项指标均优于campbell模型。因此，本发明方法能够实现将不同导热系数经验模型耦合到水热耦合模型当中，并能够实现对不同导热系数经验模型下的水热耦合模型模拟效果进行比较和评价。以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明技术原理的前提下，还可以做出若干改进和变型，这些改进和变型也应视为本发明的保护范围。}

【摘要】： 五十年代以来,应用能量观点研究非饱和土壤(包气带)水分运动,有了很大的发展。在应用这个理论对土壤水分运动进行定量分析时,求其基本方程的解或数值解,都离不开方程中出现的几个基本参数——导水率K(θ)(又叫水力传导度或非饱和土渗透系数)、扩散率D(θ)(又叫扩散度或扩散系数)和容水度C(θ)(又叫比水容量)。而且,对于某种密度的同一种土壤,它们都不是常数,是随土壤含水率θ(体积比)而变化的。另外,研究土壤水运动,对土壤水的相对能量值(土水势)的测定也是不可缺少的资料,}