寒區(qū)冰膨脹力的有限元分析

李 達(dá)，蘇安雙，孫穎娜，閻 爍

(1.黑龍江大學(xué) 水利電力學(xué)院，黑龍江 哈爾濱 150080；2．黑龍江省水利科學(xué)研究院，黑龍江 哈爾濱 150080)

在對寒區(qū)水工建筑物的影響研究中，多以研究凍土對水工建筑物凍害的研究為主，并且較為完善。魏錦等[1]總結(jié)了水工建筑物由于凍土的原因產(chǎn)生的凍害類別，以及其凍害機(jī)理，并簡述了處理凍害的防護(hù)措施；蔡梅等[2]增加了在凍融循環(huán)的過程中，水工建筑物混凝土老化引起的冰害。而對于冰對水工建筑物的損害影響研究頗少，研究課題也大多圍繞冰害防治方面；蘇盛奎[3]就冰荷載問題進(jìn)行了分析，從冰的力學(xué)性能、冰荷載的理論和計(jì)算公式以及計(jì)算實(shí)例上較為完善的闡述總結(jié)；張?jiān)衃4]就黑河市水工建筑物的冰害問題展開研究，提出利用聚苯乙烯硬質(zhì)泡沫板作為保溫材料減少靜冰壓力對水工建筑物的損害；徐秉衡[5]總結(jié)了以往對水工建筑物冰害問題的處理經(jīng)驗(yàn)，并著重對工程設(shè)計(jì)和運(yùn)行管理方面提出了有效治理冰害的措施。水工建筑物冰害的造成主要原因是冰膨脹力、動(dòng)冰撞擊力以及冰與水工建筑物界面間黏結(jié)力綜合產(chǎn)生的。因此在水工建筑物的設(shè)計(jì)過程中常常依靠經(jīng)驗(yàn)公式的計(jì)算和按照水工建筑物抗冰凍規(guī)范對冰壓力的取值和計(jì)算進(jìn)行設(shè)計(jì)，沒有從分析冰害機(jī)理層面進(jìn)行分析，以至于在復(fù)雜環(huán)境下的冰害問題上并沒有得到廣泛性的解決。

在信息技術(shù)高度發(fā)達(dá)的21世紀(jì)，利用計(jì)算機(jī)技術(shù)能夠模擬越來越多真實(shí)的情況，解決實(shí)際問題，在冰膨脹數(shù)值分析的問題上，黃焱等[6]就冰膨脹力問題展開過數(shù)值分析，但是其主要內(nèi)容主要針對就冰膨脹原理及約束條件對冰膨脹力的影響，其他因素并未詳細(xì)闡述；劉荔銘[7]以ANSYS參數(shù)化語言對冰荷載進(jìn)行數(shù)值分析，分析過程中，完善了對初始溫度、溫升率、冰蓋厚度和邊界條件四個(gè)方面進(jìn)行單獨(dú)分析，其反映溫度場變化率的過程中，圍繞的是冰蓋厚度同溫度的關(guān)系；另外李夢姌[8]借用ABAQUS有限元分析軟件對冰荷載進(jìn)行分析，除了對冰荷載影響因素單獨(dú)分析外，還利用軟件特點(diǎn)研究了具有裂紋的冰蓋在水位下降過程中對水工建筑物產(chǎn)生的危害。在以上數(shù)值分析的過程中集中對冰荷載影響因素進(jìn)行了重點(diǎn)分析，然而數(shù)值分析對實(shí)際工況的模擬還具有不完善性，以及借用有限元分析的手段來研究冰問題的精準(zhǔn)性是需要論證的。本文在介紹冰膨脹力有限元分析的同時(shí)，著重進(jìn)行數(shù)值分析。

1 冰溫度膨脹力有限元分析計(jì)算

本文使用由ANSYS公司提供的ANSYS Academic Student(ANSYS學(xué)生版)以及圖像，其中利用ANSYS WORKBENCH 19.2有限元軟件對冰溫度膨脹力(靜冰壓力)進(jìn)行分析，ANSYS WORKBENCH有限元數(shù)值分析軟件相較于ANSYS CLASSIC，具有更完善的圖形軟件接口，能夠建立較為復(fù)雜的制圖，而不出現(xiàn)丟線、丟面的狀況；同時(shí)，在做網(wǎng)格劃分、材料屬性填充和耦合場計(jì)算等數(shù)值分析時(shí)，能夠自動(dòng)選取單元、優(yōu)化劃分、操作規(guī)范化和快捷化，有效的增強(qiáng)計(jì)算能力。

1.1 有限元模型仿真

根據(jù)氣象資料，黑龍江省寒冷期極其長，多達(dá)5個(gè)月的封凍期，冬季室外溫度多數(shù)在-20 ℃左右，加上黑龍江省河流眾多，封凍期冰厚一般都在1 m左右，因此對于冰蓋的模擬仿真，為確保具有普適性，建立長100 m、寬30 m、高1 m的冰模型，將模型劃分網(wǎng)格，將模型沿厚度方向分10層,模型共劃分27 840個(gè)單元,123 733個(gè)節(jié)點(diǎn)，如圖1；另外添加冰材料屬性，見表1。

圖1 冰模型圖和有限元分析模型圖

彈性模量/MPa熱傳導(dǎo)系數(shù)/(W·m-1· ℃)比熱/(J·kg-1· ℃)密度/(kg·m-3)熱膨脹系數(shù)/(×10-6/ ℃)23000.612096.1+7.116 t(t為冰溫)916.715.2

1.2 有限元模型仿真計(jì)算

溫度場的變化和邊界條件是導(dǎo)致冰膨脹力產(chǎn)生的直接因素，若要探究冰膨脹力的影響必須從這兩方面共同入手。

1.2.1 穩(wěn)態(tài)溫度場計(jì)算

對穩(wěn)態(tài)溫度場計(jì)算是為了完成對冰膨脹力影響的最初始物質(zhì)狀態(tài)條件，其他因素的改變必須要以此為基礎(chǔ)，尤其是初始溫度改變時(shí)應(yīng)重新計(jì)算穩(wěn)態(tài)溫度場，再加以分析。分析過程如下：

設(shè)置環(huán)境溫度和上表面溫度為-20 ℃，由于下表面與水直接接觸設(shè)置為0 ℃，由此可以認(rèn)為在模型內(nèi)部形成了穩(wěn)定的溫度場，冰蓋溫度場分布，如圖2，穩(wěn)態(tài)場直線圖，如圖3，得知其分布狀態(tài)是線性分布的，而且每一平面上的溫度是相同的。

圖2 冰蓋穩(wěn)態(tài)溫度場分布效果圖

圖3 冰蓋穩(wěn)態(tài)溫度場分布直線圖

1.2.2 瞬態(tài)溫度場計(jì)算

根據(jù)黑龍江省封凍期的河流情況，在3、4月份氣溫日較差最高，達(dá)到10 ℃～24 ℃，2 h內(nèi)最高溫差也可達(dá)到10 ℃，沿用上節(jié)穩(wěn)態(tài)溫度場結(jié)果，施加動(dòng)態(tài)環(huán)境溫度和上表面為在初始溫度-20 ℃，2 h升高溫度10 ℃，下表面與水直接接觸設(shè)置為0 ℃的條件，由此可以認(rèn)為在模型內(nèi)部形成了動(dòng)態(tài)的溫度場，得到冰蓋溫度場的效果圖，如圖4所示，得知其分布狀態(tài)是非線性分布的，同樣每一平面上的溫度是相同的，由圖像可以看出，因溫升作用對冰蓋溫度場的影響區(qū)域主要分布在冰蓋上部15 cm左右的范圍，另外冰蓋內(nèi)部溫度場的溫升率變化也是由外向里逐漸減弱，直到15 cm處變化率接近于0，其他厚度下的溫度場與穩(wěn)態(tài)溫度場數(shù)值保持一致。

圖4 2.0 h冰蓋瞬態(tài)溫度場分布效果圖

1.2.3 瞬態(tài)應(yīng)力場計(jì)算

在變化的溫度場中，隨著溫度的上升，冰內(nèi)部溫度場的變化會導(dǎo)致冰體膨脹，由于邊界條件的約束而產(chǎn)生膨脹力，所以對于穩(wěn)態(tài)溫度場，內(nèi)部溫度場恒定不變，所以不產(chǎn)生溫度應(yīng)力，只有瞬態(tài)溫度場中才會產(chǎn)生溫度應(yīng)力。

直立墻約束是在冰蓋和水工建筑物作用中最常見的約束，在數(shù)值建模中設(shè)置約束條件為假定冰蓋與水工建筑物的接觸一面做全約束，然后對整個(gè)冰蓋長度和寬度方向做位移約束，厚度方向自由，在加載上面模型分析的初始溫度-20 ℃、2 h升溫10 ℃的溫度場，得到應(yīng)力場分布情況，如圖5所示，分析得知，冰蓋在厚度方向發(fā)生了較大的膨脹變形，并由于應(yīng)力集中最大主應(yīng)力發(fā)生在距離冰面1/3處，且在邊角處達(dá)到最大值，最大主應(yīng)力值為-0.271 67 MPa，為壓應(yīng)力。

圖5 直立墻約束冰蓋應(yīng)力場分布效果圖

2 結(jié)果分析

在以往的研究中，對有限單元法的適用性問題上的探究較少；就有限單元法利用有限元軟件模擬計(jì)算的精準(zhǔn)性，本文將以通過規(guī)范和經(jīng)驗(yàn)公式的方式論證其可靠性。

(1)在北方地區(qū)對水工建筑物的設(shè)計(jì)和維護(hù)中，常常無法避免對冰工程的探索，根據(jù)大量的研究試驗(yàn)和工況監(jiān)測，我國制定了中華人民共和國國家標(biāo)準(zhǔn)《水工建筑物抗冰凍設(shè)計(jì)規(guī)范》(GB/T 50662—2011)，用以對預(yù)防和設(shè)計(jì)寒區(qū)水工建筑物，本規(guī)范極具代表性和準(zhǔn)確性，規(guī)范中對于靜冰壓力部分的數(shù)值規(guī)定如下，見表2，為了驗(yàn)證對靜冰壓力有限元分析準(zhǔn)確性，應(yīng)選取表中5種厚度的冰進(jìn)行檢驗(yàn)。

(2)徐秉衡[5]在對冰害的研究過程中，總結(jié)出最大靜冰壓力經(jīng)驗(yàn)公式為：

Fd=137.3k·δ·Cδ·mt

(1)

式中：Fd為靜冰壓力，kN/m；137.3為應(yīng)力常數(shù)，kPa;δ為冰厚，m；mt為與升溫持續(xù)時(shí)間t有關(guān)的時(shí)間系數(shù)，一般天氣，δ<50 cm，t=6 h，mt=1.0，δ>50 cm，2 d連續(xù)升溫t=30 h，mt=1.82；k為綜合影響系數(shù)，一般處取4.0～4.5，庫面很大的平原水庫取5.0，小型水庫取3.5左右；Cδ為與冰厚有關(guān)的換算系數(shù)，見表3。

計(jì)算得靜冰壓力隨冰厚變化，見表4。

表2 靜冰壓力Fd

注：1.對于庫面狹小的水庫，將表中靜冰壓力值乘以0.87后采用；對于庫面開闊的大型平原水庫，將表中靜冰壓力值乘以1.25后采用。2.冰厚取多年平均最大值。3.所列靜冰壓力值系水庫在結(jié)冰期內(nèi)水位基本不變情況下的壓力，在此期間水位變動(dòng)情況下的冰壓力需做專門研究。4.靜冰壓力值按冰厚內(nèi)插確定。

表3 換算系數(shù)Cδ

表4 經(jīng)驗(yàn)公式計(jì)算靜冰壓力隨冰厚變化表

(3)依次按(1)中對其他冰厚進(jìn)行有限元分析，生成應(yīng)力場分布效果圖，如圖6～圖9，不同冰厚下最大主應(yīng)力值、轉(zhuǎn)換靜冰壓力值見表5。通過和規(guī)范數(shù)據(jù)、經(jīng)驗(yàn)公式計(jì)算對比得出，見表6，有限元分析具有一定的適用性和準(zhǔn)確性。

圖6 0.4 m冰蓋應(yīng)力場分布效果圖

圖7 0.6 m冰蓋應(yīng)力場分布效果圖

圖8 0.8 m冰蓋應(yīng)力場分布效果圖

圖9 1.2 m冰蓋應(yīng)力場分布效果圖

冰厚δ/m0.40.60.81.01.2最大主應(yīng)力σ/MPa-0.208 80-0.232 38-0.247 96-0.271 67-0.289 37靜冰壓力Fd/(kN·m-1)-83.520-139.428-198.368-271.670-347.244

表6 規(guī)范、經(jīng)驗(yàn)公式、數(shù)值分析法靜冰壓力隨冰厚變化對比表

3 結(jié) 論

(1)本文通過總結(jié)冰荷載對水工建筑物破壞的原因以及破壞的各類研究，著重就冰溫度膨脹力有限元數(shù)值分析適用性問題展開模擬試驗(yàn)和論證分析，分析得出在寒冷升溫的氣候下，冰蓋會產(chǎn)生升溫膨脹，由于約束作用，冰溫度膨脹力的應(yīng)力場隨著瞬態(tài)溫度場的變化而產(chǎn)生，數(shù)值分析得到1 m厚的冰蓋在初始溫度-20 ℃、2 h升溫10 ℃的溫度場，最大主應(yīng)力值為-0.271 67 MPa，應(yīng)力值較大，極可能造成水工建筑物的損傷和破壞。

(2)最后對不同冰厚下進(jìn)行數(shù)值分析，并對其結(jié)果同國標(biāo)規(guī)范和經(jīng)驗(yàn)公式進(jìn)行對比，事實(shí)證明有限元分析足夠作為對日后性能劣化分析的工具，但是由于實(shí)際工況的復(fù)雜性，不能夠僅借用有限元手段對工程設(shè)計(jì)的研究，仍需要利用規(guī)范設(shè)計(jì)—經(jīng)驗(yàn)公式計(jì)算—簡化試驗(yàn)—數(shù)值分析相結(jié)合的方法來解決問題。

(3)本文的數(shù)值分析及論證，一方面證明了有限元法數(shù)值分析可以充分作為工程研究的方法之一；另一方面，也提供了寒區(qū)水工建筑物在不同冰厚的情況下，冰對建筑的應(yīng)力影響范圍，可為日后水工建筑物的冰害研究、建造設(shè)計(jì)、防護(hù)維修提供有效的理論依據(jù)。