水工混凝土结构日照辐射模拟方法

颉志强，付 志，吕兴栋

(1. 长江水利委员会 长江科学院，湖北 武汉430010； 2. 水利部水工程安全与病害防治工程技术研究中心，湖北 武汉430010)

0 引 言

日照辐射对混凝土结构温度、应力时空分布的影响不容忽视，对西藏高海拔地区的水工混凝土结构尤其如此。研究表明[1-3]，日照辐射能大幅加剧结构温度场时空分布不均匀性，而后者引起的变形与结构的约束共同导致了混凝土的开裂。因此，准确掌握强日照辐射对混凝土结构的温度、应力的影响是该地区水工混凝土结构温控防裂的前提。当前的水工混凝土温控研究中，日照辐射影响分析绝大多数采用朱伯芳[1]引入的等效算法。该算法将日照辐射影响等效为气温增量，均匀施加于结构边界面，能在平均意义上反映日照辐射引起的“温升”效应，因此得到了广泛的应用并取得了大量的研究成果[4-7]。虽然等效算法能够近似反映日照辐射对结构温度场时空分布的影响，但仍有明显不足，主要表现在其无法准确反映周边环境(山体或其他建筑物)、结构的方位、朝向等因素所导致的结构受照不均匀性。也就是说，等效算法在实质上弱化了结构受照的不均匀性。为了进一步提高分析精度，陈拯[8]将ASHRAE日照辐射模型及遮蔽算法引入到拱坝运行期温度场计算中，取得了良好的效果；苏卫强等[9]对主要的日照辐射模型的适应性进行了探讨，为仿真过程中日照辐射模型的选取提供了有意义的建议。通过对现有研究成果的总结，笔者认为水工混凝土结构日照辐射模拟主要涉及3个方面的内容：①日照参数计算，主要是对不同时间、经纬度、海拔高程的日照辐射强度、日照时间、入射角进行计算；②结构受照过程模拟，主要是指对特定地区、特定形式、方位的结构在日照过程中，表面各部位受日照强度、入射角、遮蔽情况进行模拟；③日照辐射影响下温度应力仿真，主要将日照辐射影响引入混凝土结构温度场、应力场仿真计算中，评价日照辐射对结构温度、应力时空分布特性的影响。为此，笔者系统地梳理了相关计算方法，针对传统遮蔽判断算法繁琐，计算效率低的问题，引入相交快速算法提高了边界面遮蔽判断效率，简化了算法流程。开发了日照辐射精细化仿真分析程序。通过高海拔地区某电站开展的立杆实验及墩墙温度应力分析，验证了算法的合理性和准确性。

1 天文参数计算

太阳光穿过大气层后，除了散射和吸收的部分外，最终达到地面的辐射为直接辐射，大气层散射或再发射的辐射为散射辐射，从其它表面反射的辐射为反射辐射。对于混凝土结构温度场而言，太阳辐射的影响主要与结构受照表面的辐射强度和太阳入射角有关。

1.1 辐射强度

到达地球的太阳辐射以短波辐射为主，对于地表任意受照平面而言，太阳辐射主要包括：直接辐射、天空散射、地面反射3个部分。

1)直接辐射：是指太阳光经过大气直接到达地球表面的辐射强度，对于地面太阳入射角为φ的结构表面，直接辐射强度可按式(1)计算：

IDφ=IDcosφ

(1)

式中：ID为与太阳光垂直平面的辐射强度，根据KEHLBECK模型：

ID=I00.9tukam

(2)

式中：I0为太阳常数；tu为林克氏混浊度系数，随大气变化而变化；ka为不同海拔高度的相对大气压，m为大气光学质量，按式(3)～式(5)计算或取值。

(3)

(4)

(5)

式中：Atu、Btu为与计算区域相关的参数对处于山区的水利工程而言，一般Atu=2.2，Btu=0.5。

表1 林克氏混浊度系数参数取值[10]Table 1 Parameters of Rink’s turbidity coefficient

2)天空散射：是指达到地球表面的被大气层散射的太阳辐射，天空散射与受照表面的方位角无关，只与表面倾角有关，任意表面天空散射强度为[11]

(6)

式中：βn为受照平面法向与地平面的夹角；IdH为水平表面所受到的天空散射强度，按式(7)计算：

IdH=(0.271I0-0.294ID)sinβs

(7)

3)地面反射：是指直接辐射和天空散射投射到地面以后被地表反射的太阳辐射，地面反射强度根据式(8)计算[10]：

(8)

式中：re为地面短波反射率，一般地面取值0.2，积雪地面取值0.7。

由式(1)、式(6)、式(8)可知，对于混凝土结构任意受照表面而言，总的太阳辐射为

(9)

1.2 几何参数

对于地面上任意表面在空间的倾向可以由表面外法线n的方位角αn以及n与地平面的夹角βn两者完全确定。而相对于地面上的观察者，太阳在天空半球内的位置可以由太阳方位角αs和高度角βs两者所完全确定。

建立图1所示的空间直角坐标系OXYZ。坐标原点置于受照面上，OXY位于地平面内，X轴指向正南(S)，Y轴指向正东(E)，Z轴垂直地平面指向天顶。设n和s分别表示表面外法线n方向和太阳射线方向(从原点指向太阳)的单位向量。由图1，受照面的入射角φ即为向量n与向量s的夹角。根据几何关系可知，两向量可分别表示为

n=icosβncosαn+jcosβnsinαn+ksinβn

(10)

s=icosβscosαs+jcosβssinαs+ksinβs

(11)

因此，n和s之间的夹角即太阳入射角为

cosφ=n·s=cosβncosβscos(αn-αs)+ sinβnsinβs

(12)

图1 向量关系Fig. 1 Relationship of vectors

采用式(12)计算入射角时，还需确定如下参数：

1)太阳高度角βs，是指自地面观察者至太阳的连线与观察者所在的地平面的夹角，按式(13)计算：

sinβs=cosφcosδcosτ+ sinφsinδ

(13)

2)太阳方位角αs：太阳方位角αs为自地面观察者至太阳的连线在地平面上的投影与正南方向的夹角，一般规定太阳方位角以由南转向东为正，向西为负，即上午为正，下午为负，任意时刻的方位角有：

(14)

(15)

3)太阳倾角δ：亦称太阳赤纬，是指由地心指向日心的连线与地球赤道平面的夹角。太阳倾角在春分和秋分时刻等于0，夏至和冬至时刻出现极值，分别为+23.45°和-23.45°，其在一年中的变化可按经验公式(16)计算。

(16)

式中：N为日序数。

实际计算中，首先根据太阳日计算太阳倾角δ，再根据δ、计算时间τ及受照面地理纬度φ，依次计算太阳高度角βs和方位角αs，通过简单推导有[10]：

cosφ=(sinφsinβn- cosφcosαncosβn)sinδ+

(cosφcosβn- sinφcosαncosβn)cosδcosτ+

sinαncosβncosδsinτ

(17)

根据式(17)可知，在已知结构方位、角度、地理位置(经度、纬度)、日期、时间情况下，可以方便地计算出结构任意表面(非遮蔽)日照辐射入射角。

1.3 受照时段

在日照辐射精细化仿真分析中，需准确考虑日出日落时间，根据公式(13)令高度角βs=0，则特定经纬度地区日出日落时间为

(18)

基于此，可求得日出日落真太阳时tt和ts分别为

(19)

2 结构受照过程模拟算法

水利工程一般地处高山峡谷地区，在受照过程中不可避免地受到周边山体或自身不同部位的遮蔽。遮蔽效应进一步加剧了日照辐射的空间不均匀性。因此，在结构受照过程模拟中，必须对有限元模型进行实时遮蔽判断分析。

2.1 边界面遮蔽分析方法

天文参数计算完成后，即可对于有限元模型进行遮蔽判断分析，算法思路如图2，对于某待分析面，从该面沿入射光线的反向引出一条射线，判断该射线是否与模型其他面相交，若相交则该面被遮蔽，若不相交则该面未被遮蔽。例如，图2中i号面引出的一条与入射光线平行的射线，该射线与岸坡相交，说明i号面被遮蔽，而j面未被遮蔽因此可以被光线直接照射。需要注意的是，实际分析中亦可对边界节点进行逐个分析，边界节点所在引出的射线方向为其所在所有面的入射向量的均值。

图2 遮蔽效应Fig. 2 Shadow effect

结构遮蔽分析核心是线、面相交判断，计算量巨大，以图2中的小拱坝为例，边界面总数为6 086个，温度场计算中，每个时步都需要分别对6 086个临空面引出的射线与其他边界面(共6 085个)进行至少3.6亿次以上的相交判断，才能够准确分析表面受照情况。显然，在原本耗时巨大的仿真分析中，加入巨量的相交分析，既不经济也不合理。为此，文献[8]引入了栅格加速算法减少判断次数。主要是建立如图3的虚拟栅格，特定射线只需与所经栅格内的临空面进行相交判断。这种处理方式能够有效减少射线与临空面的相交判断次数，但仍然需要大量的相交计算。

图3 栅格分区Fig. 3 Grid partition

2.2 射线空间面相交快速判断

三维有限元分析绝大多数为六面体或四面体单元，因此应用中可以统一简化为射线与空间三角形的相交判断。常规做法是直接判断射线与三角形面相交，再判断交点是否在三角形内，则需要二次判断，计算过程繁琐，且效率不高。现引入一种新的射线与空间三角形相交快速判断方法，进一步提高遮蔽判断效率。

设任意射线的方程为

(20)

式中：(x0,y0,z0)为射线的起点，向量(m,n,p)为射线的方向向量。令：

则，射线方程可表示为

X=Dt+O

(21)

对于空间三角形，其内部任意点可表示为

(22)

式中：(xi,yi,zi)(i=1,3)为空间三角形的3个角点，令：

则，空间三角形可表示为

X=(1-u-v)V1+uV2+vV3

(23)

若射线与空间三角形相交，则有：

Dt+O=(1-u-v)V1+uV2+vV3

(24)

即：

Dt+u(V1-V2)+v(V1-V3)=V1-O

(25)

令，P0=V1-O,P1=V1-V2,P2=V1-V3，则有:

Dt+uP1+vP2=P0

(26)

显然，式(26)为简单的三阶线性方程组，易知式(26)的解为

(27)

利用式(27)，可以快速地判断射线是否与三角形相交。

基于上述算法，笔者采用Fortran开发了大体积混凝土结构日照辐射过程精细化仿真模拟计算程序。采用PC机(AMD Athlon(tm) II X64 2.8GHz)对图2所示拱坝进行遮蔽分析，单个时间步天文参数计算及遮蔽判断总耗时仅0.01 s左右。

2.3 试验验证

为进一步验证上述算法和程序的正确性，依托高海拔地区在建工程营地现场，开展了现场立杆试验。试验装置如图4，立杆高度2.0 m，立杆底部经纬度为北纬29°15′28″，东经92°25′12″，海拔高程3 377 m，立杆高度为2.0 m，试验时间为2017年11月3日。根据立杆尺寸，建立如图5的有限元网格，该网格节点总数为5 242，单元总数2 510，x轴正方向为正南方向。本试验主要通过模拟立杆阴影的变化验证天文参数、遮蔽算法和程序的准确性。

图4 阴影试验装置Fig. 4 Shadow test device

图6 试验与计算结果比较Fig. 6 Comparison of test and calculation results

根据计算，试验地点当天的理论日出时间为06:03，日落时间为17:23，由于两侧山体遮蔽，测点在上午09:00左右才受到日照作用，下午16:00后太阳被山体遮蔽，也无法测得阴影数据。试验测得的立杆阴影与计算得到的阴影对比见图6，其中红色短线为实测阴影，黑色区域为计算得到的阴影。对比可知，除了上午09:00计算阴影与实测阴影在长度上略有差异外，其他各整点计算结果与实测结果无论阴影角度还是阴影长度均符合良好。09:00误差主要是由阴影测量误差(上午阴影不够清晰)所致，总体而言，所述算法及开发的日照辐射模拟程序能够有效地对日照辐射过程及受照情况进行精细化仿真模拟。

3 考虑日照辐射影响的混凝土结构温度场算法

通过天文参数计算、受照过程模拟分析，确定结构的受照情况(即表面日照辐射热流密度等)后，可进一步在温度场计算中考虑日照辐射的影响。大体积混凝土温度场、应力场仿真相关理论在文献[2]中有详细的介绍，不再赘述，此处仅对温度场计算中日照辐射的精细化模拟方法进行介绍。根据热传导理论，对于结构任意边界面，热流密度可根据式(28)计算：

(28)

式中：λ为混凝土导热系数；T为混凝土温度；qs、qc、qr分别为通过日照短波辐射、表面对流换热、物体间长波辐射进入混凝土内的热流密度，对于水工混凝土结构而言，物体间的长波辐射对混凝土温度场影响很小，一般不予考虑。而对流换热引起的热流量为

qc=β(T-Ta)

(29)

考虑表面日照辐射吸收率为αs，则由于日照辐射引起的表面热流量为

qs=-αs×(IDφ+Idβ+Irβ)

(30)

进一步，在推导有限元计算格式时，对于受照边界面，利用式(31)即可。

(31)

基于此，笔者将开发的日照辐射模块引入长江科学院现有大体积混凝土温控仿真计算程序中，开发了能够充分考虑日照辐射影响的大体积混凝土温度场、应力场仿真计算程序。

4 高海拔地区混凝土结构日照响应分析

4.1 工程概况

以前述工程中某南北走向墩墙为研究对象，计算2017年7月20日全天墩墙温度场、应力场变化过程。该墙高6.0 m、墙体厚度为3.0 m、长10.0 m，墙体及地基有限元模型如图7，模型单元总数为11 590个、节点总数为13 560个。混凝土弹性模量为37.5 GPa，线膨胀系数为8.0×10-6/℃，基岩弹性模量为30 GPa。

图7 墩墙有限元模型Fig. 7 Finite element model of pier wall

为便于比较，分别计算考虑和不考虑日照辐射影响情况下，墙体一天内的温度、应力时空分布。在墙体中截面上，选取T1～T5 5个特征点，其中T1为墩墙东侧表面点、T2距东侧表面0.4 m，T3为西侧表面点，T4距西侧表面0.4 m，T5为墩墙中心点，见图8。

图8 特征点Fig. 8 Typical points

4.2 结果分析

图9～图11为考虑和不考虑日照辐射影响下，墩墙中部特征点T1～T5温度、应力历程。

图9 阳面特征点、温度应力Fig. 9 Typical points and temperature stresses in sunward side

温度方面，根据计算，考虑日照辐射以后，上午10点东侧表面点T1温度由于日照辐射引起的温差达到了16.6 ℃，下午17点左右，墩墙西侧表面点T3由于日照辐射影响温度同样达到了16.7 ℃。考虑日照影响后，两侧表面附近(T2、T4)均明显高于无日照工况，T2点最大温差2.0 ℃左右，但与表面不同，温差最大值出现的时间在20点左右，有明显的延迟。考虑日照辐射影响后，由于墙体较厚(3.0 m)墙体中心点温度受日照辐射影响不大。

应力方面，在墙体两侧(即T1、T3)，由于日照辐射导致的温度升高，表面产生了一定程度的压应力。由于日照辐射导致内外温差增大，表面附近(即T2、T4)拉应力大幅增加，以东侧T2点为例，不考虑日照影响时T2点第一主应力最大值为0.25 MPa左右，考虑日照辐射影响后，第一主应力达到了0.9 MPa左右。说明对于实际工程而言，日照辐射能够大幅增加表面混凝土开裂风险。同样，对于墩墙中心点T5，由于日照辐射影响，第一主应力最大值由原来的0.2 MPa增大至0.9 MPa以上。

计算成果表明，日照辐射明显改变了结构温度场时空分布特性，相应地，对结构温度应力时空分布也有显著影响。日照辐射的存在大幅加速了结构表面的开裂风险，对结构工作性态影响明显，因此进行日照辐射精细化模拟是必要的。

图10 阴面特征点温度、应力历程Fig. 10 Temperature and stress history of typical points in shaded side

图11 T5特征点温度、应力历程Fig. 11 Temperature and stress history of T5 typical points

5 结 论

1)准确分析了日照辐射对结构的影响，对于推进仿真分析精细化和掌握结构真实工作性态至关重要，也是当前智能建造的迫切需求。笔者系统地梳理了水工混凝土结构日照辐射仿真模拟方法，给出了相应的计算公式，并开发了程序，实现了对水工混凝土结构日照响应的精细化模拟。

2)针对传统遮蔽判断算法繁琐、效率低下的问题，引入相交判断快速算法，避免了二次判断，有效提高了遮蔽分析效率，满足了大型工程实时日照辐射响应分析对计算效率的要求。

3)算例表明，对高海拔地区水工混凝土结构而言，日照辐射会大幅加剧温度场、应力场时空分布的不均匀性，增加混凝土结构开裂风险。因此在高海拔地区水工混凝土结构温控分析中，对日照辐射过程进行精细化仿真模拟是必要的。