基于有限單元法寒區(qū)渡槽靜冰壓力分布特征研究

李可也，袁 理

（中水淮河規(guī)劃設(shè)計(jì)研究有限公司，合肥 230601）

高海拔高寒地區(qū)冬季氣溫較低， 且白天太陽(yáng)輻射強(qiáng)烈，近地面的氣溫會(huì)迅速上升，日溫差較大。由于高原低溫期持續(xù)時(shí)間較長(zhǎng)，水庫(kù)、河流、渠道內(nèi)極易結(jié)冰。研究表明，冰在氣溫升高時(shí)體積膨脹［1］，會(huì)對(duì)約束其膨脹的水工建筑物產(chǎn)生溫度膨脹力， 又稱靜冰壓力。 較大的靜冰壓力會(huì)對(duì)水工建筑物產(chǎn)生嚴(yán)重危害， 為研究高寒地區(qū)較大日溫差下的靜冰壓力分布規(guī)律， 本文選取冰蓋形狀和約束形式都比較簡(jiǎn)單的渡槽進(jìn)行仿真分析， 通過(guò)對(duì)其溫度場(chǎng)與應(yīng)力場(chǎng)的耦合計(jì)算，研究渡槽內(nèi)冰蓋在不同溫升率、不同冰厚、不同約束形式下的靜冰壓力分布。

1 靜冰壓力計(jì)算原理

靜冰壓力主要與冰蓋的溫度分布、溫升率、冰蓋厚度等因素有關(guān)， 運(yùn)用溫度加載法對(duì)靜冰壓力進(jìn)行仿真研究時(shí)，由于加載速率較慢，因此冰表現(xiàn)出一定的黏彈性［2］，黏彈性材料的溫度變化與材料變形的關(guān)系可表示為：

式中ε 為冰的線應(yīng)變；α 為冰的溫度膨脹系數(shù)，一般認(rèn)為淡水冰的溫度膨脹系數(shù)為常量；Δθ 是溫度的增量。

對(duì)式（1）的時(shí)間分量求導(dǎo)，得：

由式（2）可知，冰的應(yīng)變速率與溫升率呈線性關(guān)系，靜冰壓力的大小受溫升率的控制。

靜冰壓力求解過(guò)程中首先需要計(jì)算出冰蓋的溫度場(chǎng)， 冰蓋內(nèi)部溫度場(chǎng)的變化主要由冰表面的大氣溫度變化引起的，在有限單元法中，這是一個(gè)瞬態(tài)傳熱問(wèn)題， 根據(jù)能量守恒原理， 瞬態(tài)熱平衡可表達(dá)為（以矩陣形式表示）［3］：

式中[K]為熱傳導(dǎo)矩陣，包含導(dǎo)熱系數(shù)、對(duì)流系數(shù)及輻射率和形狀系數(shù)；[C]為比熱矩陣，考慮系統(tǒng)內(nèi)能的增加；{T}是節(jié)點(diǎn)溫度向量；{T·}為溫度對(duì)時(shí)間的導(dǎo)數(shù)；{Q}為冰點(diǎn)熱流向量，包含熱生率。

計(jì)算渡槽內(nèi)冰蓋溫度場(chǎng)時(shí)， 認(rèn)為淡水冰蓋下表層與水的接觸面的溫度為其融化點(diǎn)0 ℃［4］，冰蓋上表層溫度按第三類邊界條件確定， 即冰層上表面的溫度與大氣溫度相同。第三類邊界條件可表示為：

式中 Tf為熱流介質(zhì)溫度；Γ 為物體邊界；ζ 為熱換系數(shù)。

求得冰蓋溫度場(chǎng)后， 將有限元模型中各節(jié)點(diǎn)的溫度作為荷載引入應(yīng)力場(chǎng)計(jì)算中， 通過(guò)定義材料的熱膨脹系數(shù)α 將溫度場(chǎng)與應(yīng)力場(chǎng)進(jìn)行耦合計(jì)算，最后得到冰蓋內(nèi)的應(yīng)力分布規(guī)律， 作用于邊界上的冰的溫度膨脹力即為冰對(duì)水工建筑物的靜冰壓力［4］。

通過(guò)原型觀測(cè)，計(jì)算冰蓋應(yīng)力場(chǎng)時(shí)，做以下兩點(diǎn)假定［5］：

（1）計(jì)算初始時(shí)，冰蓋內(nèi)部無(wú)溫度膨脹產(chǎn)生的應(yīng)力，邊界上的靜冰壓力為零。

（2）認(rèn)為冰蓋邊界與渡槽凍結(jié)，仿真計(jì)算時(shí)，約束類型為四周全約束。

2 直立墻渡槽內(nèi)的靜冰壓力

為研究高寒地區(qū)靜冰壓力分布規(guī)律， 本文以地處高原地區(qū)平均海拔2300 m 的某縣某直立墻渡槽為例進(jìn)行了靜冰壓力的仿真研究，據(jù)氣象資料查得，該縣近年來(lái)測(cè)得的年最低溫度為-20 ℃，最大日溫升達(dá)2 h 升高10 ℃。冰蓋各項(xiàng)參數(shù)取值如下：

表1 冰蓋計(jì)算參數(shù)

2.1 冰蓋溫度場(chǎng)分布

取直立墻渡槽內(nèi)2 m×8 m×0.4 m 的冰蓋為研究對(duì)象，為取得較精確的結(jié)果，將冰蓋沿厚度分為20層，模型共劃分單元32000 個(gè)，節(jié)點(diǎn)139041 個(gè)，冰蓋有限元模型如圖1。首先計(jì)算冰蓋溫度場(chǎng)分布，設(shè)冰蓋上表面度為-20 ℃，且在2 h 升高10 ℃，冰蓋下表面溫度為0 ℃，計(jì)算初始時(shí)，認(rèn)為冰蓋已形成穩(wěn)定的溫度場(chǎng)。將穩(wěn)態(tài)溫度場(chǎng)作為初始條件進(jìn)行迭代，最終計(jì)算得到溫升2 h 過(guò)程的瞬態(tài)溫度場(chǎng)分布， 穩(wěn)態(tài)溫度場(chǎng)云圖和沿程瞬態(tài)溫度場(chǎng)曲線如圖2、圖3。

圖1 冰蓋三維有限元計(jì)算模型

圖2 冰蓋穩(wěn)態(tài)溫度場(chǎng)分布云圖

圖3 冰蓋沿程瞬態(tài)溫度場(chǎng)分布曲線

由圖2 和圖3 可看出，渡槽內(nèi)的穩(wěn)態(tài)溫度場(chǎng)呈線性分布，隨著時(shí)間的增加，冰蓋上表面溫度的升高，冰蓋溫度場(chǎng)產(chǎn)生非線性變化，但由于冰的導(dǎo)熱性較差，溫差的升高集中在冰蓋上部1/3以上處，且從表層向下變化程度逐漸減小，瞬態(tài)溫度場(chǎng)在約15 cm 處以下基本與穩(wěn)態(tài)溫度場(chǎng)重合。

2.2 靜冰壓力分布

得到冰蓋的溫度場(chǎng)之后， 將熱單元轉(zhuǎn)換為結(jié)構(gòu)單元。由于假定冰層內(nèi)部初始應(yīng)力為零，這就要求給材料賦予線膨脹系數(shù)，設(shè)定材料參考溫度時(shí)，將穩(wěn)態(tài)溫度場(chǎng)計(jì)算時(shí)求出的單元溫度值分布賦予應(yīng)力場(chǎng)分析時(shí)對(duì)應(yīng)的單元。 然后將瞬態(tài)溫度場(chǎng)計(jì)算得到的各步溫度作為荷載施加到應(yīng)力計(jì)算的模型中， 經(jīng)過(guò)迭代得出各節(jié)點(diǎn)的應(yīng)力分量［6］。由于所求結(jié)果為壓應(yīng)力分布，將Y 方向的應(yīng)力作為輸出結(jié)果進(jìn)行分析，直立墻渡槽內(nèi)0.4 m 厚冰蓋的靜冰壓力云圖分布如圖4～圖7，冰蓋沿程瞬態(tài)應(yīng)力場(chǎng)分布曲線如圖8。

圖4 0.5h 時(shí)冰蓋的應(yīng)力場(chǎng)分布云圖

圖5 1h 時(shí)冰蓋的應(yīng)力場(chǎng)分布云圖

圖6 1.5h 時(shí)冰蓋的應(yīng)力場(chǎng)分布云圖

圖7 2h 時(shí)冰蓋的應(yīng)力場(chǎng)分布云圖

圖8 冰蓋沿程瞬態(tài)應(yīng)力場(chǎng)分布曲線

由圖8 可看出，隨著冰蓋表層溫度的緩慢上升，冰蓋產(chǎn)生的溫度膨脹力受到周圍約束的作用產(chǎn)生靜冰壓力。靜冰壓力沿程分布不均，集中在冰蓋上部1/3以上處，下部趨近于零，表層靜冰壓力最大，靜冰壓力隨著溫度的不斷升高持續(xù)增大，2 h 的最大靜冰壓力值為268.7 kPa，這與溫度場(chǎng)仿真的結(jié)果相對(duì)應(yīng)，即溫度變化最大處的靜冰壓力也最大。

2.3 溫升率對(duì)靜冰壓力的影響

為了研究溫升率對(duì)靜冰壓力的影響，假定大氣溫度在 2 h 內(nèi)從-20 ℃升高至-15 ℃，其他條件不變， 得到冰蓋的應(yīng)力場(chǎng)仿真結(jié)果，2 h時(shí)兩種溫升率下的Y 方向靜冰壓力沿程分布曲線如圖9。

圖9 兩種溫升率下2h 時(shí)冰蓋沿程靜冰壓力分布曲線

由圖9 可看出，兩種溫升率下的靜冰壓力均集中于冰蓋1/3 以上處，靜冰壓力最大處同樣位于冰蓋表層。溫升率的大小決定了靜冰壓力的大小，隨著溫升率的增大，靜冰壓力也變大，符合之前的理論分析。

2.4 冰厚對(duì)靜冰壓力的影響

為研究不同冰蓋厚度對(duì)靜冰壓力的影響， 根據(jù)該縣當(dāng)?shù)卣鎸?shí)的冰清資料進(jìn)行假定， 當(dāng)冰蓋厚度分別為0.3，0.4，0.5，0.6，0.7，0.8 m， 按2 h 溫升10 ℃，其他條件不變時(shí)，冰蓋表層最大靜冰壓力值分布如圖10。

圖10 不同厚度冰蓋表層靜冰壓力分布曲線

由圖10 可看出，冰蓋表層靜冰壓力的最大值隨著冰蓋厚度的增大而增大。且由仿真分析可知，不同厚度的冰蓋，其靜冰壓力仍集中作用于冰蓋表層值冰厚15cm處，說(shuō)明冰蓋表層的溫升產(chǎn)生的靜冰壓力僅對(duì)冰蓋上部產(chǎn)生比較大的影響，對(duì)冰蓋下部影響很小。

3 結(jié)語(yǔ)

（1）采用有限元法對(duì)冰蓋的靜冰壓力進(jìn)行仿真計(jì)算具有相當(dāng)精度， 在進(jìn)行高寒地區(qū)水工建筑物設(shè)計(jì)時(shí)， 其表層最大靜冰壓力仿真結(jié)果可作為估算冰蓋靜冰壓力極值的參考。

（2）冰蓋瞬態(tài)溫度場(chǎng)變化集中于冰蓋上部15cm以上處。

（3）溫升率對(duì)冰蓋表面靜冰壓力有較大影響，靜冰壓力隨溫升率增加而增加。

（4）冰蓋厚度對(duì)表層最大靜冰壓力值有一定影響，但靜冰壓力仍集中于冰蓋上部15cm 以上處。