頡志強，付 志，呂興棟

(1. 長江水利委員會 長江科學院，湖北 武漢430010； 2. 水利部水工程安全與病害防治工程技術研究中心，湖北 武漢430010)

0 引 言

日照輻射對混凝土結構溫度、應力時空分布的影響不容忽視，對西藏高海拔地區(qū)的水工混凝土結構尤其如此。研究表明[1-3]，日照輻射能大幅加劇結構溫度場時空分布不均勻性，而后者引起的變形與結構的約束共同導致了混凝土的開裂。因此，準確掌握強日照輻射對混凝土結構的溫度、應力的影響是該地區(qū)水工混凝土結構溫控防裂的前提。當前的水工混凝土溫控研究中，日照輻射影響分析絕大多數采用朱伯芳[1]引入的等效算法。該算法將日照輻射影響等效為氣溫增量，均勻施加于結構邊界面，能在平均意義上反映日照輻射引起的“溫升”效應，因此得到了廣泛的應用并取得了大量的研究成果[4-7]。雖然等效算法能夠近似反映日照輻射對結構溫度場時空分布的影響，但仍有明顯不足，主要表現在其無法準確反映周邊環(huán)境(山體或其他建筑物)、結構的方位、朝向等因素所導致的結構受照不均勻性。也就是說，等效算法在實質上弱化了結構受照的不均勻性。為了進一步提高分析精度，陳拯[8]將ASHRAE日照輻射模型及遮蔽算法引入到拱壩運行期溫度場計算中，取得了良好的效果；蘇衛(wèi)強等[9]對主要的日照輻射模型的適應性進行了探討，為仿真過程中日照輻射模型的選取提供了有意義的建議。通過對現有研究成果的總結，筆者認為水工混凝土結構日照輻射模擬主要涉及3個方面的內容：①日照參數計算，主要是對不同時間、經緯度、海拔高程的日照輻射強度、日照時間、入射角進行計算；②結構受照過程模擬，主要是指對特定地區(qū)、特定形式、方位的結構在日照過程中，表面各部位受日照強度、入射角、遮蔽情況進行模擬；③日照輻射影響下溫度應力仿真，主要將日照輻射影響引入混凝土結構溫度場、應力場仿真計算中，評價日照輻射對結構溫度、應力時空分布特性的影響。為此，筆者系統(tǒng)地梳理了相關計算方法，針對傳統(tǒng)遮蔽判斷算法繁瑣，計算效率低的問題，引入相交快速算法提高了邊界面遮蔽判斷效率，簡化了算法流程。開發(fā)了日照輻射精細化仿真分析程序。通過高海拔地區(qū)某電站開展的立桿實驗及墩墻溫度應力分析，驗證了算法的合理性和準確性。

1 天文參數計算

太陽光穿過大氣層后，除了散射和吸收的部分外，最終達到地面的輻射為直接輻射，大氣層散射或再發(fā)射的輻射為散射輻射，從其它表面反射的輻射為反射輻射。對于混凝土結構溫度場而言，太陽輻射的影響主要與結構受照表面的輻射強度和太陽入射角有關。

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為與計算區(qū)域相關的參數對處于山區(qū)的水利工程而言，一般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為自地面觀察者至太陽的連線在地平面上的投影與正南方向的夾角，一般規(guī)定太陽方位角以由南轉向東為正，向西為負，即上午為正，下午為負，任意時刻的方位角有：

(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，則特定經緯度地區(qū)日出日落時間為

(18)

基于此，可求得日出日落真太陽時tt和ts分別為

(19)

2 結構受照過程模擬算法

水利工程一般地處高山峽谷地區(qū)，在受照過程中不可避免地受到周邊山體或自身不同部位的遮蔽。遮蔽效應進一步加劇了日照輻射的空間不均勻性。因此，在結構受照過程模擬中，必須對有限元模型進行實時遮蔽判斷分析。

2.1 邊界面遮蔽分析方法

天文參數計算完成后，即可對于有限元模型進行遮蔽判斷分析，算法思路如圖2，對于某待分析面，從該面沿入射光線的反向引出一條射線，判斷該射線是否與模型其他面相交，若相交則該面被遮蔽，若不相交則該面未被遮蔽。例如，圖2中i號面引出的一條與入射光線平行的射線，該射線與岸坡相交，說明i號面被遮蔽，而j面未被遮蔽因此可以被光線直接照射。需要注意的是，實際分析中亦可對邊界節(jié)點進行逐個分析，邊界節(jié)點所在引出的射線方向為其所在所有面的入射向量的均值。

圖2 遮蔽效應Fig. 2 Shadow effect

結構遮蔽分析核心是線、面相交判斷，計算量巨大，以圖2中的小拱壩為例，邊界面總數為6 086個，溫度場計算中，每個時步都需要分別對6 086個臨空面引出的射線與其他邊界面(共6 085個)進行至少3.6億次以上的相交判斷，才能夠準確分析表面受照情況。顯然，在原本耗時巨大的仿真分析中，加入巨量的相交分析，既不經濟也不合理。為此，文獻[8]引入了柵格加速算法減少判斷次數。主要是建立如圖3的虛擬柵格，特定射線只需與所經柵格內的臨空面進行相交判斷。這種處理方式能夠有效減少射線與臨空面的相交判斷次數，但仍然需要大量的相交計算。

圖3 柵格分區(qū)Fig. 3 Grid partition

2.2 射線空間面相交快速判斷

三維有限元分析絕大多數為六面體或四面體單元，因此應用中可以統(tǒng)一簡化為射線與空間三角形的相交判斷。常規(guī)做法是直接判斷射線與三角形面相交，再判斷交點是否在三角形內，則需要二次判斷，計算過程繁瑣，且效率不高?，F引入一種新的射線與空間三角形相交快速判斷方法，進一步提高遮蔽判斷效率。

設任意射線的方程為

(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開發(fā)了大體積混凝土結構日照輻射過程精細化仿真模擬計算程序。采用PC機(AMD Athlon(tm) II X64 2.8GHz)對圖2所示拱壩進行遮蔽分析，單個時間步天文參數計算及遮蔽判斷總耗時僅0.01 s左右。

2.3 試驗驗證

為進一步驗證上述算法和程序的正確性，依托高海拔地區(qū)在建工程營地現場，開展了現場立桿試驗。試驗裝置如圖4，立桿高度2.0 m，立桿底部經緯度為北緯29°15′28″，東經92°25′12″，海拔高程3 377 m，立桿高度為2.0 m，試驗時間為2017年11月3日。根據立桿尺寸，建立如圖5的有限元網格，該網格節(jié)點總數為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，其中紅色短線為實測陰影，黑色區(qū)域為計算得到的陰影。對比可知，除了上午09:00計算陰影與實測陰影在長度上略有差異外，其他各整點計算結果與實測結果無論陰影角度還是陰影長度均符合良好。09:00誤差主要是由陰影測量誤差(上午陰影不夠清晰)所致，總體而言，所述算法及開發(fā)的日照輻射模擬程序能夠有效地對日照輻射過程及受照情況進行精細化仿真模擬。

3 考慮日照輻射影響的混凝土結構溫度場算法

通過天文參數計算、受照過程模擬分析，確定結構的受照情況(即表面日照輻射熱流密度等)后，可進一步在溫度場計算中考慮日照輻射的影響。大體積混凝土溫度場、應力場仿真相關理論在文獻[2]中有詳細的介紹，不再贅述，此處僅對溫度場計算中日照輻射的精細化模擬方法進行介紹。根據熱傳導理論，對于結構任意邊界面，熱流密度可根據式(28)計算：

(28)

式中：λ為混凝土導熱系數；T為混凝土溫度；qs、qc、qr分別為通過日照短波輻射、表面對流換熱、物體間長波輻射進入混凝土內的熱流密度，對于水工混凝土結構而言，物體間的長波輻射對混凝土溫度場影響很小，一般不予考慮。而對流換熱引起的熱流量為

qc=β(T-Ta)

(29)

考慮表面日照輻射吸收率為αs，則由于日照輻射引起的表面熱流量為

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

(30)

進一步，在推導有限元計算格式時，對于受照邊界面，利用式(31)即可。

(31)

基于此，筆者將開發(fā)的日照輻射模塊引入長江科學院現有大體積混凝土溫控仿真計算程序中，開發(fā)了能夠充分考慮日照輻射影響的大體積混凝土溫度場、應力場仿真計算程序。

4 高海拔地區(qū)混凝土結構日照響應分析

4.1 工程概況

以前述工程中某南北走向墩墻為研究對象，計算2017年7月20日全天墩墻溫度場、應力場變化過程。該墻高6.0 m、墻體厚度為3.0 m、長10.0 m，墻體及地基有限元模型如圖7，模型單元總數為11 590個、節(jié)點總數為13 560個?；炷翉椥阅Ａ繛?7.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以上。

計算成果表明，日照輻射明顯改變了結構溫度場時空分布特性，相應地，對結構溫度應力時空分布也有顯著影響。日照輻射的存在大幅加速了結構表面的開裂風險，對結構工作性態(tài)影響明顯，因此進行日照輻射精細化模擬是必要的。

圖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)準確分析了日照輻射對結構的影響，對于推進仿真分析精細化和掌握結構真實工作性態(tài)至關重要，也是當前智能建造的迫切需求。筆者系統(tǒng)地梳理了水工混凝土結構日照輻射仿真模擬方法，給出了相應的計算公式，并開發(fā)了程序，實現了對水工混凝土結構日照響應的精細化模擬。

2)針對傳統(tǒng)遮蔽判斷算法繁瑣、效率低下的問題，引入相交判斷快速算法，避免了二次判斷，有效提高了遮蔽分析效率，滿足了大型工程實時日照輻射響應分析對計算效率的要求。

3)算例表明，對高海拔地區(qū)水工混凝土結構而言，日照輻射會大幅加劇溫度場、應力場時空分布的不均勻性，增加混凝土結構開裂風險。因此在高海拔地區(qū)水工混凝土結構溫控分析中，對日照輻射過程進行精細化仿真模擬是必要的。