中图分类号：TV315 文献标识码：A 文章编号：0559-9350（2011）06-0737-06
Determination of imPACt factors on temperature field in dam based on global sensitivity analysis method
WU Zhi-wei1，2，SONG Han-zhou1
Abstract：Coupled model of temperature field and seepage field in dam contains many parameters under operation. Sensitivity analysis method can be used to determine the main impact factors，which is able to reduce workload of model calibration. Based on a steady heat-fluid coupled model，both local sensitivity analysis and global sensitivity analysis were carried out. A case study is taken as an example. The results show that，the temperatures at the bottom of a reservoir and on the downstream dam surface are most sigNIficant factors and others are in the sequence of apparent hydraulic gradient， specific permeability， and conductive coefficient of medium，factor of porosity and heat capacity. It also shows that the global sensitivity has advantages in analyzing interactions among parameters.
Key words：temperature field in dam；local sensitivity；global sensitivity；numerical calculation
1 研究背景
正常运行期的大坝坝基与坝体（堆石坝、碾压混凝土坝等）温度场受渗流条件的影响。实测资料表明，大坝渗流温度场各处差异与水文地质条件较吻合，监测大坝温度场的变化可以反演大坝渗流场的分布及大坝渗透缺陷［1］。基于此，许多学者建立了反映大坝温度场与渗流场相互作用的耦合计算模型［2-3］。此类耦合模型均包含诸多参数，在进行模型校正时，工作量较大。为了节省计算资源，有必要对模型做灵敏度分析，定量研究各影响因素对模型的影响，找出决定温度场变化的主要因素。
这样就能集中主要精力提高这些主要影响因素的观测精度，即使较少考虑或忽略那些次要因素，也可以得到模拟结果与实际现象接近的数值模型。
灵敏度分析方法是用来研究系统影响因子在微小摄动时系统的变异情况，可以定量分析各影响因子对系统的贡献［4］。目前灵敏度分析方法在诸多领域都有广泛的应用［5-8］，对于大坝温度场，已有文献多是将灵敏度分析方法应用于大体积混凝土浇注的温度应力控制研究［9］。而针对运行期受渗流影响的大坝温度场，同样可以采用该方法探讨其主要影响因子，但相关的研究还少见报道。本文主要采用灵敏度分析方法探讨运行期坝址温度场的主要影响因子。
灵敏度分析包含局部灵敏度分析和全局灵敏度分析。本文结合大坝热流耦合模型，引入两种不同的灵敏度分析方法来探讨运行期大坝温度场的主要影响因子，为针对大坝温度场的模型校正提供借鉴，并结合工程实例探讨两种分析方法各自的优缺点。
2 大坝热-流耦合模型
首先做以下基本假设：稀浓度溶液，液-固相间不发生质量交换；流体密度只受温度和压强的影响；水流运动符合Darcy定律；含水层及坝体是均质、多孔介质；各物理场都是稳态的。
常规耦合都假设流体是定常的，实际上流体的密度和黏度受温度的影响较显著。浅部地下水的温度多在0~40℃之间，根据该温度范围内水的物性资料（据Weast R C［10］做整理），得如下关系：
μ(T) =0.001×(1.7920-0.0620T +0.0017T 2 -3.5481×10-5T3 +5.9912×10 -7T 4 -3.0042×10-9 T 5) （1）
ρw (T) =-0.0088T 2 +0.0851T +999.7 0℃< T <4℃（2）
ρw (T) = -0.0052T 2 - 0.0236T +1000.1 4℃≤T≤40℃（3）
式中：温度（T）的单位为℃；密度（ρw）的单位为kg/m3；动力黏度（μ）的单位为Pa·s。
式（1）—（3）的相关系数R2都在0.9998以上，可见能准确反映动力黏度和密度与温度的关系。
在含水层内任取一无限小的平行六面均衡单元体，根据质量守恒定律，无源汇项的稳态渗流满足：
式中：v 为地下水流速。
如果考虑温度对流体密度的影响，流体密度在各处不均一，式（4）中密度的偏微分项不可省去。利用文献［11］的Darcy流速表示方法，可得流体密度黏度非均一的地下水稳定渗流方程：
式中：p 为压强；g 为重力加速度；k 为介质渗透率；j 为笛卡尔坐标系，分别表示x、y、z 方向；e1=e2=0，e3=1，将其代入式（3），可得多孔介质稳态对流传热方程为：
式中：λeq为地层等效导热系数， λeq =φλw + (1 - φ)λs ，λs、λw分别为固体介质和水的导热系数，φ为孔隙率。
式（5）与式（6）分别为渗流场和温度场控制方程（详细推导过程见文献［12］），在耦合计算时调用式（1）、式（2）（或式（3））以更新求解域中流体的物理性质，结合相应的边界条件可实现双场耦合。
3 灵敏度分析
灵敏度分析包括局部灵敏度分析和全局灵敏度分析。
3.1 局部灵敏度分析
局部灵敏度分析只检验单个参数的变化对模型结果的影响程度。分析时只改变一个参数的值，而其它参数均取其中心值，计算在该参数发生变化时的模型变化量来衡量模型对其的灵敏度。针对不同的问题有两种变换法：第一种是因子变化法。如将预分析的参数增加10％或减少10％，计算模型的改变量；另一种方法是偏差变化法，如将预分析的参数增加一个标准偏差或减少一个标准偏差。前者操作较为方便，后者只适用于参数取值符合一定概率分布的情况。
3.2 全局灵敏度分析
全局灵敏度分析能检验模型多个参数的变化对其输出结果产生的总的影响，并能分析每一个参数及参数之间的相互作用所产生的影响。
关于全局灵敏度分析方法的研究成果已很多，但由于全局灵敏度分析方法难以操作，将其引入到大坝温度场数值模拟中的研究成果还未见报道。全局灵敏度分析方法主要有多元回归法、Morris法、傅里叶幅度灵敏度检验法（FAST）以及基于方差分析的Sobol法等［13］。本文采用Morris法［14］做大坝温度场的全局灵敏度分析。Morris法是假设数值模型中有k 个参数，每个参数的取样点个数为p、m 个参数分别在p 个取样点上取值，就可以获得向量X=［x1、x2、…、xk］，构造m×k（m=k+1）阶矩阵B：
矩阵的第j 列代表第j 个参数，矩阵值1和0分别代表参数值改变与未改变。各行参数组合的计算值之差就是相应两个或多个参数同时改变的全局灵敏度［12］。
4 工程实例
4.1 模型设定
以文献［2］中的工程实例作为本文的计算模型。
（1）几何模型。龙滩碾压混凝土（RCC）重力坝位于红水河上，最大坝高216.5m。选取此坝的第12号坝段进行灵敏度分析，计算模型如图1所示。碾压混凝土坝不同于常态混凝土坝，它本身是多孔隙低透水性介质，因此坝体内部也会形成渗流场，这里将其与坝基部分同时分析。坝体上游布有排水孔，根据经验，坝体水头在排水孔附近骤降，因此假设坝体渗流自由面形态近似用CJ表示。这样在自由面之上的坝体部分就可以不考虑连续渗流。
图1 计算模型（引自文献［2］，A~L 为边界编号）
（2）初始计算参数。对原文中各计算参数的单位进行折算，取15℃时的地下水密度和黏度把渗透系数折算为渗透率。得到模型所采用的参数值，如表1所示。
表1 计算参数（据文献［2］）
（3）边界条件与初始条件。如图1所示，渗流场边界条件为： H | DE =150m 、H | HIJ =0m ；DC的水头按其埋深呈线性变化；EF、FG、GH以及z=0m和z=3m断面为零通量面（q=0）；BCJ和JI分别为自由面和溢出面边界。温度场边界条件为：上游坝面B 处温度为多年平均表层水温21℃，水深60m以下保持12.1℃不变，水深0~60m之间按线性变化；下游坝面按多年平均气温并考虑太阳幅射作用取为22.1℃；EF、FG、GH以及x=0m和x=3m断面为绝热面。初始渗流场水头取零，初始温度场温度取15℃。
4.2 计算方法
由于考虑的是一个耦合问题，因此大坝温度场应该同时受到区域热状态和渗流场的共同影响。在计算区域相当大的假设条件下，可以认为大坝渗流体系与深部地质体之间没有热量交换，则温度边界只有上游库底水温和下游坝面温度两个温度已知边界条件。渗流介质的基本参数对温度场有一定的影响，可分别用渗透率、比热容、孔隙率和热传导系数等加以反映。另外上下游的水头边界也是温度场的影响因素，这里用水力梯度i=（H 上-H 下）/L 来综合反映水头边界条件，L 是坝体底宽。因此，选用的大坝温度场摄动因子有：上游库底水温T 上；下游坝面温度T 下；渗透率k；介质比热cs；孔隙率φ；饱和热传导系数λ；水力梯度i。
因为各物理量的单位不尽相同，必须对评价指标归一化，这里的做法是用改变值和初始值的差与初始值的比值，以百分数来衡量。
额定状态即4.1节所用参数和边界条件所对应的大坝温度场状态。
大坝温度场具有空间特性，由于数据量较大，所以取一个横向断面来分析。这里只分析纵向断面，把因参数变动的温度场与额定温度场在每个有限元单元节点上做差值，然后取平均来表示温度的变化。
4.3 局部灵敏度分析结果评价
局部灵敏度分析的具体方法是假设各摄动因子波动±5%、±10%、±15%、±20%，应用数值模拟方法得到各参数变动后的温度场分布，与额定温度场作对照，即可揭示大坝温度场对各参数的灵敏性，如图2所示。
图2 参数波动与模型结果的关系
上游温度波动±20%时，温度场分别波动约+12.3%和-12.2%，它是对温度场影响最为显著的指标；下游温度波动±20%时，温度变幅约为±6.6%，与上游温度边界条件相比，大坝温度场对其灵敏度减半。其它参数变动相同的幅度，温度场改变都在±1%以下。在这些参数中，反映流场水力条件的上下游水力梯度对温度分布影响相对较强，为+0.6%~-0.8%，另外几个与介质特性相关的参数对温度场分布贡献甚微。由温度波动范围的大小可以确定上游水库库底温度和下游坝面温度对温度场的影响最明显，其次依次为水力梯度（i）、渗透率（k）、岩石热导系数（λ）、孔隙率（Φ）和介质比热（cs）。造成这种现象的原因可认为是温度边界相当于温度场的稳定热源，决定系统能量的大小，而流场条件和介质热学参数主要决定热量运移过程，在这个稳态问题中难以体现出来。
各参数的变动与温度场的变化基本都是线性关系，与温度场正相关的参数有T 上、T 下和λ，即参数值变大，温度变大，反之，则减小；与大坝温度场负相关的参数有k、φ和i，这都符合定性分析结果。
局部灵敏度分析的优点是易于操作，但也存在一些缺陷，如一次只能分析一个参数，不能考虑参数之间相互作用对模型输出结果的影响。另外，对某一参数进行分析时，其它参数都取中间值。实际上，其它参数的不同取值会影响该参数的灵敏度，模拟结果的变化是所有参数共同作用产生的，所以局部灵敏度分析在计算分析上存在一定的局限性。
4.4 全局灵敏度分析结果分析
由局部灵敏度分析可知各参数的灵敏度函数多为线性，采用Morris法做全局灵敏度分析时，可只假设参数增加10%，这样既可以得到可信的分析结果又能简化计算。分别将矩阵B 各行参数组合代入到温度场耦合计算模型中，并依次计算模拟结果之差，可得各参数及参数间相互作用的灵敏度。根据一种参数排列方式，由Morris法得28种参数组合的全局灵敏度，如表2所示。
为了更加直观地反映分析结果，这里用条形图来显示各种组合方式的温度场变化幅度（%），见图3。由图3与表2可知：
（1）组别1-7表示在全局灵敏度分析时只有单个参数变化所引起的模型变动，其总体规律与局部灵敏度分析相似。模型对各参数的灵敏度大小顺序一致，但两种方法的模型变幅有小幅差异，这是因为前者考虑了其它参数不同取值对模型结果的影响。
（2）组别8、14、19、23、26和28是温度场变幅较大的几种参数组合形式，而它们的参数组合中都同时包含了上游温度和下游温度，说明多个参数相互组合影响时，上游温度和下游温度仍是决定大坝温度场的主要因素。然而从具体数值来看，这6种组合方式对温度场的作用幅度又是不同的，这反映了各种参数会相互作用，共同影响温度系统的响应。
（3）所有参数同时增大10%的系统温度变动（组别28）反而小于组别8、23、26。由局部灵敏度分析可知，参数与大坝温度场存在正相关和负相关两种关系，因此，可以认为这种现象由多个参数相互制约引起的，但数值上又不是简单的代数叠加。
（4）结合局部灵敏度分析知，文中所选的7 个参数中，与系统正相关的3 个参数对系统作用明显，那么在正相关的这些参数中可以舍去部分对系统作用相对较弱的参数。而负相关的参数虽然绝对值较小，但由全局灵敏度分析知，它们所带来的系统变动能够与正相关参数的作用相中和，对系统的影响也会很显著，因此在做模型校正时不能忽略。
全局灵敏度分析结果证实，其它参数的不同取值对某个参数灵敏度计算是有影响的，采用全局灵敏度分析得到的单个参数灵敏度与局部灵敏度分析结果不尽相同，显然全局灵敏度分析结果更加可信。参数对模型结果的影响也不能用简单的数学方法加以叠加，而全局灵敏度能够计算多个参数共同作用对模型结果的影响，在做模型校正过程中有其优越性。
5 结语
（1）两种灵敏度分析方法都说明，运行期大坝温度场受上游水库库底温度和下游坝面温度的影响最明显，其次依次为水力梯度（i）、渗透率（k）、岩石热导系数（λ）、孔隙率（Φ）和介质比热（cs）。与温度场呈负相关的是水力梯度、渗透率和孔隙率，其余为正相关；（2）全局灵敏度分析表明多参数共同影响下的模型响应并不是单个参数作用效果的简单叠加，全局灵敏度分析能用来定量研究参数间共同作用对数值模型的影响，与局部灵敏度分析方法相比有其优越性，应加以推广；（3）本文研究对象是运行状态下的大坝温度场，并对物理模型做了诸多简化处理，但就寻找大坝温度场主要影响因子来说，这种简化是可行的。
表2 基于全局灵敏度分析的温度变幅（%）
图3 基于全局灵敏度分析的温度变幅
参考文献：
［ 1 ］肖才忠，潘文昌. 由温度场研究坝基渗流初探［J］. 人民长江，1999，30（5）：21-23 .
［ 2 ］柴军瑞. 混凝土坝渗流场与稳定温度场耦合分析的数学模型［J］. 水力发电学报，2000（1）：27-35 .
［ 3 ］许增光，李康宏，柴军瑞，等. 考虑渗流热学效应的大坝稳定温度场有限元数值分析［J］. 红水河，2006，25（2）：112-115 .
［ 4 ］ Crosetto M，Tarantola S . Uncertainty and sensitivity analysis：Tools for GIS-based model implementation［J］. International Journal of Geographical Information Science，2001，l5（5）：4l5-437 .
［ 5 ］李森，陈家军，叶慧海，等. 地下水流数值模拟中随机因素的灵敏度分析［J］. 水利学报，2006，37（8）：977-984 .
［ 6 ］王浩昌，杜鹏飞，赵冬泉，等. 城市降雨径流模型参数全局灵敏度分析［J］. 中国环境科学，2008，28（8）：725-729 .
［ 7 ］李国敏，Chin-Fu Tsang . 地下非均匀非饱和带中地下洞室的渗流问题数值模拟——介质参数的灵敏度分析［J］. 地球科学，2003，28（5）：497-504 .
［ 8 ］徐崇刚，胡远满，常禹，等. 生态模型的灵敏度分析［J］. 应用生态学报，2004（6）：1057-1062 .
［ 9 ］朱伯芳，等. 大体积混凝土的温度应力与温度控制［M］. 北京：中国电力出版社，1999 .
［ 10］ Weast R C . Handbook of Chemistry and Physics，63rd ed .［M］. Boca Raton . Fla . ：Chemistry Rubra Pun . Co .Press，1982：261-263 .
［ 11］张勇，薛禹群，谢春红. 高温差条件下达西定律的理论推导［J］. 水科学进展，1999，10（4）：362-367 .
［ 12］吴志伟，宋汉周. 坝址温度场与变物性渗流场全耦合分析［J］. 水利学报，2010，41（6）：703-710 .
［ 13］束龙仓，王茂枚，刘瑞国，等. 地下水数值模拟中的参数灵敏度分析［J］. 河海大学学报（自然科学版），2007，35（5）：491-495 .
［ 14］ Morris M D . Factorial sampling plans for preliminary computational experiments［J］. Technometrics，1991，33 （2）：161-174 .
作者简介：吴志伟（1985-），男，河南信阳人，博士生，主要从事地下水、应用地球物理探测研究。
