自实施西部大开发、水资源高效利用等国家战略以来，我国已修建了一批特高拱坝[1-2]．我国特高拱坝库区地形地质条件极为复杂，工程结构性态易受到各种特殊因素影响[3]．变形能够直观反映大坝的安全状态．引起坝体变形的温度荷载，主要由大坝温度场分布变化情况决定．目前，当计算大坝温度场分布变化时，多考虑水温和气温的影响[4]，而常忽略太阳辐射的影响．国内外众多学者围绕太阳辐射下大坝温度场分布及变形性态，开展了一系列研究．Mirzabozorg等[5]探究了太阳辐射对薄高拱坝温度分布的影响；Jin等[6]采用ASHRAE(美国采暖、制冷与空调工程师学会)晴空模型和光线追踪法，准确计算了坝体所受太阳辐射强度与坝体温度场；Santillán等[7]提出了一种计算拱坝温度荷载的方法，该方法考虑了太阳辐射、夜间冷却和蒸发冷却，量化分析了不同温度场计算方法对应力和变形的影响；Žvanut等[8]通过建立相对简单的混凝土坝传热过程模型，考虑了水库水位、溢出、日照和阴影等对大坝温度场的影响；Pan等[9]当计算太阳辐射对拱坝温度场的影响时，另外考虑了气候条件对坝体受到的太阳辐射的影响．上述研究虽然从不同角度考虑了影响坝体温度场的因素，但较少探究太阳辐射对库水温度的影响，降低了大坝温度场边界条件分析的准确性．太阳辐射会对特高拱坝温度场分布产生影响，从而导致特高拱坝变形性态发生变化．为准确掌握特高拱坝运行性态，有必要针对太阳辐射对坝体变形性态的影响进行重点探究．本研究考虑特高拱坝所处地理位置、地形特征及气候条件，联合ASHRAE模型和光线追踪法计算坝面和库水实际接收的太阳辐射强度．根据特高拱坝狭长河道径流特征，运用水库水温三维计算方法模拟了太阳辐射下库水温变化情况．1 特高拱坝坝体及库水所受太阳辐射强度计算方法1.1　ASHRAE晴空模型ASHRAE晴空模型可用于求解特高拱坝及库水所受的瞬时太阳辐射强度．该模型主要从直射辐射、天空散射辐射和反射辐射三个方面求解太阳辐射强度[10]．入射到非垂直表面的总太阳辐射强度      Gt1=[max(cosθ,0)+CrFws+ρgFwg(sinβ+Cr)]GND,式中：θ为太阳光线入射角，是太阳光线与平面倾斜法线之间的夹角；Cr代表水平面接收的太阳散射辐射强度与垂直入射的直射辐射强度两者的比值；Fws为表面与天空间的角系数；ρg为周围环境的反射率；Fwg为表面对地面的角系数；β为太阳高度角；GND为垂直入射情况下的太阳直射辐射强度．入射到垂直表面的总太阳辐射强度[8]为      Gt=[max(cosθ,0)+GdVCr/GdH+ρgFwg(sinβ+Cr)]GND,式中：GdV为垂直表面天空散射辐射强度；GdH为水平表面天空散射辐射强度；若cosθ-0.2，则GdV/GdH =0.55+0.437cosθ+0.313cos2θ，否则GdV/GdH =0.45．1.2　光线追踪法由于两岸山体和坝体自身均会遮蔽一部分太阳辐射，因此将导致特高拱坝坝面和库水接收到的太阳辐射分布不均．为更准确地模拟出山体及坝体自身对太阳辐射实际的遮蔽情况，本文引入光线追踪算法[11]，对不同时刻坝体及库水实际接收到的太阳辐射进行计算．计算特高拱坝山体及坝体自身对太阳辐射遮蔽情况的过程[12]如下．a．依据特高拱坝三维有限元模型坝面单元，确定与空气接触的单元边，并称其为自由边．b．基于模型的拓扑信息和节点坐标，即可求得自由边相关参数．c．将特高拱坝三维有限元模型所在空间等分成正方形网格，并按坐标对其编号．依据每一个正方形网格包含的自由边编号，建立链表．d．指定任意的日期和时间，计算当下太阳高度角和太阳表面方位角．e．依据每一条自由边的中心点坐标及步骤d计算出的结果，产生射线．f．计算射线所经过的正方形网格编号．g．当射线与所经过正方形网格的自由边相交时，说明该自由边所受太阳辐射会被遮蔽；否则，不会被遮蔽．h．按步骤e～g即可判断所有自由边所受太阳辐射是否被遮挡，据此得到特高拱坝坝面及库水所受太阳辐射遮蔽情况．1.3　预计算技术当计算特高拱坝及库水各时刻所受太阳辐射强度时，须考虑每一时刻太阳位置、水位变化情况、太阳辐射遮蔽情况和天气条件等信息；同时，特高拱坝模型结构复杂，当采用光线追踪法计算任意时刻太阳辐射强度和日照遮蔽情况时，求解过程十分复杂．为降低计算复杂度，有效简化计算过程，并保证计算结果的准确度，本研究引入预计算技术，预先计算典型日的遮蔽情况和特征水位下的水温修正信息，在实际计算中运用预计算结果，即可插值得到邻近时间点坝面及库水所受太阳辐射强度．基于预计算技术的太阳辐射强度计算流程如下．a．根据特高拱坝所处地区的经纬度计算任意时刻太阳高度角和方位角，代入ASHRAE晴空模型计算太阳辐射强度．b．选取一年内每月的15日作为该月典型日，通过查询当天不同时刻的太阳位置计算24 h内坝面的太阳辐射遮蔽情况，其中被遮蔽区域赋值为1，未被遮蔽区域可赋值为0．上游坝面水位以下不考虑太阳辐射，被遮挡区域不受太阳辐射影响．遮挡区域和水位以下区域的判定均通过有限元软件自带的二次开发程序接口实现．假定任意时刻下坝体遮挡情况随不同日期是线性变化的，则可通过对前后两个典型日的遮蔽区域进行[0，1]区间的插值计算，求解任意时刻不同日期的特高拱坝及库水太阳辐射强度的遮蔽情况．c．以晴朗天气下太阳辐射强度为基本太阳辐射强度，确定晴天、多云及其他(阴、雨、雪)天气情况下对应的太阳辐射折减系数分别为1.0，0.5，0.0．d．根据历史水位数据得到任意时刻上游水位，选择一年冬至、春分、夏至、秋分四个典型日和四个不同的典型水位，计算16个典型场景的水体所受太阳辐射的遮蔽情况，其余场景根据线性插值获得．e．根据三种不同的天气状况确定太阳辐射强度的折减程度，运用有限元数值分析方法计算16个典型场景在不同天气状况下的水温变化情况，对原始的水温数据进行修正．据此可以得到一个三维的查询库，通过代入日期、水位、天气三个参数即可插值得到对应的水温随深度的修正值．f．基于该地区实际气象数据，求解实际天气情况下特高拱坝坝体及库水所受太阳辐射强度折减系数，确定坝面及库水实际受到的太阳辐射强度．2 太阳辐射下特高拱坝温度场及变形数值分析方法2.1　太阳辐射下库水温度数值分析2.1.1　水库表面的热交换由于特高拱坝上游库水面积较大，受太阳照射时间较长，因此须要考虑太阳辐射对库水温的影响以修正库水温度．特高拱坝库水表面的热通量的计算公式为φ=RS+RA-RB-Ew-C，式中：RS为库水所受短波辐射强度；RA为库水所受大气的长波辐射强度；RB为库水长波返回辐射强度；Ew为库水蒸发热损失量；C为实际的热传导通量．一般情况下，上述热交换过程只在水面发生，水面以下水体认为只接收太阳的短波辐射[13]．2.1.2　水温对流扩散方程为了更好地反映特高拱坝狭长河道和地形特征的影响，以及紊动扩散各向异性的特点，引入σ坐标系，垂向σ坐标系用(x,y,σ,Tw)表示．将太阳辐射单独从水温对流扩散方程的热源项中分离出来，则水温对流扩散方程可表示为      ∂Tw∂t+u∂Tw∂x+v∂Tw∂y+ω∂Tw∂σ=∂∂xAH∂Tw∂x+∂∂yAH∂Tw∂y+∂∂σAVH2∂Tw∂σ+1ρwcpH∂RSz∂σ+QT,式中：x,y和σ分别为水平方向和垂向坐标；u,v和ω分别表示x,y和σ方向的流速；h为静水面到水底的水深值；H=h+ς为总水深，其中ς为水面偏移静水面的位移；t为时间；ρw为水体密度；Tw为水温；AH和AV为水平温度扩散系数和垂向温度扩散系数；RSz为太阳辐射强度随库水深的衰减量；QT为源汇项；cp为水体比热容．2.1.3　库水温4类边界条件的确定库水温4类边界条件如图1所示：①为来流边界条件；②为出流边界条件；③为水面边界条件；④为底部边界条件．通过对水动力控制方程和水温对流扩散方程进行数值计算，即可得出库水温度分布情况[14-15]．10.13245/j.hust.230060.F001图1库水温4类边界条件示意图2.2　特高拱坝温度场数值计算采用瞬态温度场来分析坝体温度变化情况，首先基于ASHRAE晴空模型及光线追踪算法计算出坝体所受太阳辐射强度及实际遮蔽情况，之后根据预计算技术求解出的任意时刻太阳辐射遮蔽情况、水温修正值及太阳辐射强度折减系数，考虑太阳辐射、环境温度、库水温度影响，引入热传导方程，通过确定坝体温度场4类边界条件，即可计算出任意时刻坝体瞬态温度场[16-17]．特高拱坝坝体瞬态温度场的有限元形式的表达式为Ctϕ˙+Ktϕ=P，式中：Ct为热容矩阵；Kt为热传导矩阵；P为温度荷载矩阵；ϕ为节点温度列阵；ϕ˙为节点温度对时间的导数列阵，ϕ˙=dϕ/dT．特高拱坝坝体温度场的4类边界条件如图2所示，边界条件S1~S4分别指与库水接触坝面S1、与空气接触坝面S2、地基基底面S3与地基侧壁面S4．10.13245/j.hust.230060.F002图2特高拱坝坝体温度场4类边界条件示意图2.3　太阳辐射影响下特高拱坝变形数值计算特高拱坝在太阳辐射、库水温度、气温等多因素综合影响下，坝体会产生相应的温度变形．基于线弹性本构模型，各节点温度变形等效于体积力与面力作用产生的变形，则变温节点的等效荷载列阵可表示为      YT=αE1-2μ∑i=1nNiTi∙[1,1,1,0,0,0]T,式中：E为坝体弹性模量；α为线膨胀系数；μ为泊松比；Ni为单元形函数；Ti为单元温度；i为单元编号．运用有限元计算坝体温度变形的平衡方程为K⋅δT=YT，式中：K为材料劲度矩阵；δT为温度变形列阵．太阳辐射影响下特高拱坝温度变形数值分析过程如图3所示．10.13245/j.hust.230060.F003图3考虑太阳辐射影响的特高拱坝温度变形数值分析流程3 工程实例分析3.1　有限元模型构建我国某特高拱坝为混凝土双曲拱坝，位于雅砻江大河，其设计坝高达305 m[18]．以该特高拱坝为例，考虑太阳辐射的影响，对其温度变形进行数值计算分析．首先建立该特高拱坝的有限元数值分析模型，该模型以特高拱坝拱冠梁为中心分别向上游、下游、左岸与右岸延伸1 500 m，坝基以下延伸1 200 m．模型利用八节点二次四面体单元，单元总数为9.525 5×104，节点总数为1.106 48×104，其中坝体单元为5.108 3×104．坝体单元网格的水平向与竖向尺寸均为9~10 m，地基单元网格水平向与竖向尺寸由下至上从100 m逐渐变为10 m．坐标原点设在拱冠梁建基面中点，切向指向左岸为正，径向指向下游为正，垂向向上为正．模型精确反映了拱坝的体型和周围山体形状等，以便根据光线追踪法实时计算出遮挡范围．特高拱坝有限元模型如图4所示．因为下游水体面积较小，所以主要考虑上游库水温度受太阳辐射的影响．10.13245/j.hust.230060.F004图4特高拱坝工程有限元模型3.2　计算参数和环境量选取为计算坝面和库水实际接收的太阳辐射，首先须要根据该特高拱坝工程概况计算出任意时刻太阳的高度角和方位角，并带入ASHRAE晴空模型中计算得到任意时刻的太阳辐射强度．已知：该特高拱坝最大坝高为305 m；坝底高程为1 580 m；坝顶高程为1 885 m；拱冠梁顶厚为16 m；拱冠梁底厚为63 m；正常蓄水位为1 880 m；总库容为7.76×109m3；坝址经度为28°10'49″N，纬度为101°37'59″E；拱冠梁上游面法线朝向为25°NE．在ASHRAE晴空模型中，大气清洁度CN根据大坝所处位置实际情况估算为0.9，周围环境的反射率ρg=0.25，混凝土吸收系数ξ=0.65．表1列出了该特高拱坝坝体与地基温度场计算参数的取值，表中：λ为导热系数；c为比热容；a1为表面传热系数(空气)；a2为表面传热系数(水)．10.13245/j.hust.230060.T001表1坝体混凝土及地基温度场计算参数部位λ/(J∙(m∙s∙℃)-1)c/(J∙(kg∙℃)-1)a1/(W∙(m2∙K)-1)a2/(W∙(m2∙K)-1)坝体混凝土1.395902.32500地基1.914766.02500主要考察坝体温度场的分布情况，混凝土和地基全部采用线弹性本构模型，其材料参数如表2所示，表中：E为弹性模量；μ为泊松比；ρ为密度．10.13245/j.hust.230060.T002表2坝体混凝土及地基材料参数部位E/GPaμ线膨胀系数/10-5ρ/(kg∙m-3)混凝土A区30.50.16712 400混凝土B区30.20.16712 400混凝土C区25.00.16712 400地基Ⅱ类岩30.00.200—2 700地基Ⅲ类岩21.00.250—2 650统计坝址多年各月平均气温、水温，以及坝前断面在春、夏、秋、冬典型日库水温度的垂向分布情况．3.3　边界条件和预计算信息设定在该特高拱坝模型底端加载固定约束，周围加载法向约束．水体覆盖的表面施加水温对流传导边界，其余表面施加气温对流传导边界，未被遮挡的区域施加太阳辐射边界，其中水温边界也须考虑太阳辐射的影响．根据前述的计算方法，首先计算出坝面和库水在典型日所受太阳辐射遮蔽情况及折减情况．在计算得到典型日各时刻太阳辐射遮蔽情况后，通过数值模拟可以得到相应的库水温修正信息如图5所示，图中θ1为温度．10.13245/j.hust.230060.F005图5典型日太阳辐射影响下库水温度修正曲线由于太阳辐射对库水温度的影响随着深度的增加而快速衰减，因此一般认为太阳辐射只会影响表层库水温度．由于采集的实测库水温通常为日平均库水温度值，难以准确反映每日每时刻太阳辐射影响下库水温度的实际情况，因此须要根据实测库水温对各时刻库水温度进行修正．其中，根据夏季与冬季某典型日实测垂向库水温度值对某时刻垂向库水温进行修正的结果如图6所示．10.13245/j.hust.230060.F006图6依据垂向库水温度实测值计算出的不同时刻修正值1—库水温度实测值；2—库水温度修正值．因为太阳辐射对坝体温度场的影响相对于环境温度的影响较小，所以可首先计算出常规气温、水温荷载下的基础温度场，根据实测气温和水温数据确定该特高拱坝温度场计算的边界条件S1~S4．其中，夏季与冬季的基础温度场分别如图7(a)和(b)所示．10.13245/j.hust.230060.F007图7特高拱坝基础温度场3.4　特高拱坝温度场及变形数值计算分析通过运用预计算技术求解太阳辐射遮蔽情况可知，冬季坝体下游面只有在中午才在左岸有很小的区域被照射，上游库水只有在下午较短的时间被太阳照射，经计算表明：冬季太阳辐射的影响很小，可以忽略不计．夏季坝体下游面左岸受太阳的照射较为强烈，且上游水面被太阳的直射时间很长，故夏季太阳辐射对坝体温度场及变形性态的影响很大．图8(b)为2018年夏至日晴朗天气下经过一天太阳照射后的坝体温度场剖面图，可见太阳辐射使得坝面温度上升2~3 ℃，局部上升4~5 ℃，太阳辐射主要对坝面温度影响较大．10.13245/j.hust.230060.F008图82018年夏至日拱冠梁剖面温度场图9为考虑太阳辐射影响，进行水温修正前后，夏至日库水温度变化对坝体温度场的影响结果，可以看出：在水温修正后，只会造成接近坝顶的坝面温度有一定的增加．总体而言，太阳辐射造成的上游水温的变化对坝体温度场产生的影响较小．10.13245/j.hust.230060.F009图9夏至日库水温度修正前后的坝体温度场坝体在2018年夏至日晴朗天气太阳照射下升温2~5 ℃，坝体产生向上游最大的变形值为0.42 mm，如图10所示，左岸受到的太阳照射较强烈，所以整体变形中心偏向左岸．10.13245/j.hust.230060.F010图10太阳辐射影响下产生的夏至日径向变形场最终，经过夏季连续晴天太阳辐射连续影响下，并考虑气温、水温共同作用，计算出的坝体径向变形如图11所示，最大变形出现在拱冠梁处，产生的倾向上游变形量为12.61 mm．因此，受夏季连续晴天太阳辐射影响产生的坝体径向变形约占总变形量的1/3．由此可见：当研究特高拱坝坝体变形性态时，不可以忽视太阳辐射的影响．10.13245/j.hust.230060.F011图11夏季考虑气温、水温、太阳辐射连续影响的大坝径向变形场此外，太阳辐射也会对坝体温度应力分布产生影响．由于在上游坝面，太阳辐射仅会对库水面以上及以下一定深度处的坝面温度产生影响，因此如图12(a)和图13(a)所示，上游坝面在库水面附近其应力分布会产生突变．由于山体及拱坝自身对太阳辐射的遮蔽作用，因此导致左岸坝顶附近温度受太阳辐射的影响最为显著，其左岸坝顶附近温度最高．10.13245/j.hust.230060.F012图12受夏季连续晴天太阳辐射影响后坝面最大主应力分布图10.13245/j.hust.230060.F013图13受夏季连续晴天太阳辐射影响后坝面最小主应力分布图由图12(b)和图13(b)可以看出：下游坝面由此产生左右岸不对称的应力分布，靠近左岸坝顶附近的压应力最大．4 结论a．上游坝面水位以下的特高拱坝坝体温度主要受库水温度及其波动的影响，并且在相同高程上具有均匀的温度分布，水位以上坝面主要受太阳辐射和气温影响，同时由于拱坝自身与山体对太阳辐射的遮蔽作用，导致上游坝面所受太阳辐射强度和温度场分布不均．b．本文算例分析选取的特高拱坝的下游坝面温度场分布主要受太阳辐射和气温影响，同时由于山体及坝体自身对太阳辐射的遮蔽作用，导致左岸温升较大，使左岸下游坝面产生较大压应力，产生整体倾向上游的温度变形．针对位于不同位置和不同朝向的特高拱坝，其太阳辐射影响下的温度场、温度应力及变形性态分布规律均可运用本文方法进行分析与研究．c．对于本文算例分析选取的特高拱坝，当考虑太阳辐射影响计算各垂线测点变形值时，计算出的夏季测点径向变形值更接近于实测值，而计算出的冬季考虑与不考虑太阳辐射时计算出的测点径向变形值均与实测值较为接近．这是由于夏季特高拱坝各部位受到的太阳辐射强度较强，而冬季大坝各部位受到的太阳辐射强度较弱，因此在分析夏季特高拱坝变形性态时须要着重考虑太阳辐射的影响，而在分析冬季特高拱坝变形性态时可基本忽略太阳辐射的影响．