高溫混凝土內(nèi)部溫度場計算方法研究

楊海林

（唐山工業(yè)職業(yè)技術(shù)學(xué)院，河北 唐山 063299）

近年來，我國建筑火災(zāi)頻發(fā)，給人民財產(chǎn)造成了不可估量的損失。在危害建筑結(jié)構(gòu)的諸多因素中，火災(zāi)已成為危害性最大的災(zāi)害之一。當(dāng)前我國的建筑大多采用鋼筋混凝土結(jié)構(gòu)，火災(zāi)產(chǎn)生的高溫環(huán)境會使結(jié)構(gòu)內(nèi)部發(fā)生物理化學(xué)變化，從而影響結(jié)構(gòu)的受力性能，嚴(yán)重的甚至?xí)斐山Y(jié)構(gòu)破壞倒塌。結(jié)構(gòu)抗火性能研究是一門復(fù)雜的交叉學(xué)科，混凝土結(jié)構(gòu)應(yīng)用于建筑工程的初期沒有得到足夠的重視，直到1905年美國的ASTM（American Society for Testing and Materials）協(xié)會才研究了墻板構(gòu)件的高溫試驗方法，但并沒有得到廣泛應(yīng)用，直到40年后國際上才對混凝土結(jié)構(gòu)抗火性能重視起來。其中研究較早的國家主要有美國、加拿大、法國、瑞士等，隨著大量理論和試驗研究的展開，混凝土結(jié)構(gòu)的抗火性能研究等進(jìn)入到了快速發(fā)展階段。我國在這方面的研究較晚，20世紀(jì)60年代公安部的消防研究所和冶金部的建筑研究院等單位才開始了混凝土結(jié)構(gòu)耐火性能研究。20世紀(jì)80年代，同濟、清華等高校相繼建立了高溫試驗裝置，建立了高溫下混凝土結(jié)構(gòu)的本構(gòu)模型，取得了大量的研究成果。火災(zāi)發(fā)展一般會經(jīng)歷3個階段，按照先后順序依次為增長階段、全盛階段、衰減階段，研究高溫混凝土的性能首先要研究溫度場。

近年來，許多研究人員對溫度場進(jìn)行了大量研究，成果主要有：劉猛等[1]運用ABAQUS建立了混凝土高溫框架模型，對構(gòu)件受力及內(nèi)部溫度場分布規(guī)律進(jìn)行了研究；王勇等[2]通過對常溫和高溫連續(xù)板進(jìn)行力學(xué)試驗，分析了鋼筋混凝土高溫狀態(tài)下的破壞模式；吳慶良等[3]針對沉管隧道受火前后的變化，應(yīng)力2次分離變量法，研究了RABT火災(zāi)曲線下溫度場的求解方法；宋巖升等[4]應(yīng)用PyroSim分析軟件對某宿舍樓火災(zāi)進(jìn)行了模擬，研究了結(jié)構(gòu)內(nèi)部的溫度場分布及內(nèi)力變形；習(xí)陽等[5]應(yīng)用ANSYS建立了鋼殼混凝土模型，對結(jié)構(gòu)不同受火時間下的溫度場變化進(jìn)行了分析。目前對高溫下混凝土內(nèi)部溫度場的研究已經(jīng)有了一定的基礎(chǔ)，但缺少更加細(xì)致的成果，沒有熱傳導(dǎo)微分方程不同計算方法的對比研究。因此，工程中亟需一種能夠快速準(zhǔn)確地確定構(gòu)件火災(zāi)時溫度場分布的方法。在查閱相關(guān)文獻(xiàn)的基礎(chǔ)上[6-8]，通過大量的理論計算，對建筑火災(zāi)過程中溫度場隨時間的變化規(guī)律進(jìn)行了研究，對3種計算方法的優(yōu)缺點進(jìn)行了對比分析，為今后的研究提供了基礎(chǔ)。

1 火災(zāi)溫度-時間曲線

由于火災(zāi)發(fā)生時溫度場變化非常復(fù)雜，根據(jù)相關(guān)統(tǒng)計，地上建筑發(fā)生的火災(zāi)中，80%在1 h內(nèi)被撲滅，95%在2 h內(nèi)被撲滅。為了便于科學(xué)研究，國內(nèi)外許多機構(gòu)制定了火災(zāi)溫度隨時間變化的標(biāo)準(zhǔn)曲線，這些曲線都是根據(jù)大量試驗數(shù)據(jù)繪制的，曲線形態(tài)都是單調(diào)的升溫曲線，升溫速率隨著時間降低，國際標(biāo)準(zhǔn)化組織ISO（International Organization for Standardization）提出的升溫曲線如圖1所示，其計算公式見式（1）。

圖1 ISO834標(biāo)準(zhǔn)升溫曲線

公式（1）中t為火災(zāi)發(fā)生時間，To（℃）為初始溫度，T為火災(zāi)發(fā)生t分鐘后空氣的平均溫度。

歐洲的Eurocode1規(guī)范[9]綜合火災(zāi)現(xiàn)場各種因素給出了室內(nèi)溫度-時間升溫段（t*≤td*）的公式（2）和降溫段（t*＞td*）的公式（3）。

式（2）中T（℃）為某一時刻的溫度，t*（min）為名義升溫時間，td*（min）為名義升溫時間最大值，其余含義同公式（1）。

式（3）中Tmax（℃）為最高溫度，其余含義同公式（2）。

2 溫度場熱傳導(dǎo)微分方程

2.1 基本方程推導(dǎo)

要想研究火災(zāi)對混凝土結(jié)構(gòu)的影響，就要研究火災(zāi)發(fā)生的任一時刻空間各點的溫度分布，但溫度場是動態(tài)發(fā)展的。由于火焰是向上蔓延的，所以建筑受火后樓板和梁的底部為主要受火面，溫度發(fā)展最快，此時可按公式（1）和公式（2）進(jìn)行計算。地面和樓面為背火面，短時間內(nèi)溫度不會明顯升高，為簡化計算，通常按室溫計算。

為了分析混凝土內(nèi)部的熱傳導(dǎo)，可以建立固體物質(zhì)的微分方程，假設(shè)材料各向同性，且λ（熱導(dǎo)率）、c（比熱容）、ρ（表觀密度）都是T的函數(shù)，微元體的熱流分析見圖2。在如圖2所示的直角坐標(biāo)系中，在任一點（x,y,z）附近取一微元體dx, dy, dz，假設(shè)在dx, dy, dz內(nèi)部溫度T均勻分布，某一時刻t的溫度為T（x,y,z,t）。

圖2 微元體熱流分析圖

如圖2所示，微元體與外部存在熱傳遞，首先研究X軸方向吸收的熱量。單位時間單位面積流入微元體的熱量為qx，按公式（4）計算，流出的為qx+dx按公式（5）計算。

由于微元體平行于YOZ平面的面積為dy·dz，所以單位時間內(nèi)X軸方向吸收的熱量為式（6）。

同理在Y軸方向吸收的熱量見式（7），在Z軸方向吸收的熱量見式（8）。

將式（6）、式（7）、式（8）相加可以得到微元體在單位時間內(nèi)吸收的總熱量qA，即公式（9）。

假設(shè)一定體積一定時間內(nèi)材料內(nèi)部放出的熱量為qd，則微元體放熱總量為qddxdydz。吸收熱量會使微元體溫度升高，還可以根據(jù)比熱容的定義進(jìn)行計算，即吸收的總熱量為（c·pdxdydz）。根據(jù)能量守恒，得到式（10）、式（11），即瞬態(tài)熱傳導(dǎo)微分方程。

由于混凝土結(jié)構(gòu)自身產(chǎn)生的熱量遠(yuǎn)小于吸收的熱量，所以可以簡化為qd=0如果構(gòu)件內(nèi)部外部溫度都不隨時間T變化，則公式（11）可簡化為公式（12），即得到了穩(wěn)態(tài)熱傳導(dǎo)微分方程。

2.2 不同構(gòu)件的適用情況

不同的混凝土構(gòu)件試用不同的微分方程，要綜合考慮構(gòu)件類型、維度及狀態(tài)。桿系結(jié)構(gòu)（如梁、柱等）可以按二維溫度場計算，平面結(jié)構(gòu)（如墻、板等）可簡化為一維溫度場，具體結(jié)果見表1。

表1 不同構(gòu)件的熱傳導(dǎo)微分方程適用情況匯總

如果材料的導(dǎo)熱系數(shù)λ不隨外界溫度變化，可引入常數(shù)，其中d(m2/h、m2/s)為材料的熱擴散率，可通過熱工參數(shù)進(jìn)行推導(dǎo)，將其帶入式（11）、式（12），可簡化為式（13）、式（14）。

2.3 定解條件和求解方法

要想求解熱傳導(dǎo)微分方程，還需要得到對應(yīng)的邊界條件，即傳熱過程中材料與周邊熱量相互作用的條件。當(dāng)材料給定時，可以通過試驗測得對應(yīng)的熱工參數(shù)，包括λ（熱導(dǎo)率）、c（比熱容）、p（表觀密度）等，確定初始和邊界條件后就可以對方程進(jìn)行求解，下面討論如何確定初始和邊界條件。

混凝土結(jié)構(gòu)在火災(zāi)發(fā)生前一直處于室溫狀態(tài)，內(nèi)部幾乎不存在溫差，所以可以將環(huán)境溫度T0作為初始溫度條件，即T(x,y,z,t=0)=T0。

邊界條件由于受到外部環(huán)境、熱傳導(dǎo)條件、周邊介質(zhì)等的影響會變得復(fù)雜，一般有4種情況。

1）l1（結(jié)構(gòu)邊界）上的溫度T和時間t存在函數(shù)關(guān)系，即。

2）l2（結(jié)構(gòu)邊界）上的即熱流量）和t存在函數(shù)關(guān)系，其中n為邊界的法線方向。

3）結(jié)構(gòu)和流體接觸，已知Ta（外部接觸的流體溫度）和經(jīng)過l3（邊界）上的熱流量，即，其中βT為傳熱系數(shù)，單位W/(m2·K)。

4)結(jié)構(gòu)和固體接觸，已知l4上的熱交換條件，受火面可按式（15）計算，其中α(W·m-2·℃-1)為換熱系數(shù)，為結(jié)構(gòu)表面溫度，Te(℃)為環(huán)境溫度。具體到混凝土結(jié)構(gòu)，可將式（15）換算成式（16），其中hc(W·m-2·℃-1)為對流換熱系數(shù)，hr(W·m-2·℃-1)為輻射換熱系數(shù)， =5.67×10-8W·m-2·℃-4。

實際工程中為了簡化計算，可將hr換算成hc，將輻射影響計算到α（綜合換熱系數(shù)）內(nèi)，即簡化為公式（17）,其中K(W·m-2·℃-1)為傳熱系數(shù)，α可按表2選用，為了便于計算任一溫度下的α值，將表2中數(shù)據(jù)進(jìn)行擬合，見圖3。

表2 綜合換熱系數(shù)α(W·m-2·℃-1)選用表

圖3 換熱系數(shù)α隨火焰溫度變化圖

采用Origin中的ExpDec2進(jìn)行擬合，發(fā)現(xiàn)R2達(dá)到了0.99768，表明擬合相關(guān)性較好。根據(jù)圖3的擬合結(jié)果，可以得到任一T下α的值，用于火災(zāi)發(fā)生時受火面的溫度場計算，為了簡化計算，背火面通常取室溫。

3 求解方法

熱傳導(dǎo)微分方程為二階偏微分方程，對其求解的方法主要有以下3類。

3.1 解析法

為了簡化計算，通常假定混凝土材料內(nèi)部是各向同性、連續(xù)、內(nèi)部不發(fā)熱的，由于內(nèi)部鋼筋占比較少，通常不考慮其影響，工程中的板和梁高跨比較大，通常簡化為無限大板和無限長梁，此時公式（13）就可以簡化為一維或二維問題，即式（18）。

假定初始溫度為20℃，由式（18）和式（15）組成方程組，即式（19）。

對于無限大平板，可進(jìn)行拉普拉斯變換得到式（20），其中y為受火面到內(nèi)部的垂直距離，erf(x)為高斯誤差函數(shù)，可查表獲得，=1/（pc）（m2/s）為導(dǎo)溫系數(shù)。

此外還可以將受火面，即y=0時的溫度變化簡化為正弦函數(shù)，采用分離變量法求解，即式（21），其中T'max為溫度增量的最大值。

值得注意的是，式（20）和式（21）由于做了大量簡化計算，只適用于初步粗算，達(dá)不到較高精度。

3.2 差分法

對于形狀規(guī)則邊界簡單的結(jié)構(gòu)，可以用差分代替微分，同時考慮節(jié)點作用，達(dá)到減小計算量的目的，二維狀態(tài)下的溫度場網(wǎng)格劃分見圖4，顯示差分方程可表示為式（22）。

圖4 二維溫度場差分法網(wǎng)格劃分示意圖

對式（22）進(jìn)行求解可采用泛函數(shù)求極值的方法，用歐拉方程求解得到極值曲線。采用泛函極值條件的二點差分熱傳導(dǎo)方程可表示為式（23），式中[K]、[N]、[P]為溫度矩陣，{T}t為待解溫度列向量，{T}t-△t為上一時刻列向量。應(yīng)用式（23）可依次求得間隔△t各時刻的溫度場。為了簡化計算程序，后期步長可適當(dāng)放大，同時放大背火面步長。

在實際工程中，還要綜合考慮火災(zāi)發(fā)生時混凝土內(nèi)部水分的影響，為了簡化分析，通常將這些水分的影響綜合到α、λ、c、p中進(jìn)行計算。綜上，給定邊界條件的情況下，截面形狀規(guī)則的構(gòu)件（如三角形、矩形、正方形）應(yīng)用差分法求解溫度場具有較高的精度，此方法計算靈活，速度較快，但不適于處理復(fù)雜形狀的構(gòu)件。

3.3 有限元-差分法

有限元-差分法目前應(yīng)用廣泛，可用于求解任意時刻的溫度場問題，包含兩個部分，即空域采用有限元法，時域采用差分法。目前的大型軟件如ANSYS、ABAQUS等都可進(jìn)行有限元-差分法的運算，計算時沿時間方向遞推進(jìn)行迭代運算，分析軟件會自動根據(jù)迭代情況進(jìn)行收斂。為了簡化計算，將瞬態(tài)熱傳導(dǎo)微分方程構(gòu)造成公式（24）所示的線性函數(shù)，其中{T}為溫度向量，{Pβ}為熱載荷向量，[c]為熱容矩陣，[K]為總導(dǎo)熱矩陣。進(jìn)行計算時，可先給定收斂容許值Tg，當(dāng)時運算停止。

式（24）中，Nk=1-ξ，Nk+1=ξ，ξ=τ/△tk且0≤ξ≤1。對式（24）中的{T}求時間t的偏導(dǎo)可得式（25）。

如圖1所示，標(biāo)準(zhǔn)升溫曲線為對數(shù)函數(shù)，溫度升高速率逐漸降低，這與構(gòu)件溫度場的升溫規(guī)律相近。當(dāng)升溫速率較大時（即升溫初期），運算過程中的振蕩劇烈，反之（即升溫末期）振蕩趨于平緩，為了提高運算效率，可以隨溫度升高適當(dāng)延長步長，進(jìn)行變步長運算。

式（24）與式（25）可以用于瞬態(tài)熱傳導(dǎo)的求解，除此之外還有穩(wěn)態(tài)熱傳導(dǎo)，下面介紹穩(wěn)態(tài)熱傳導(dǎo)微分方程的解法，方程見式（26）。

式（26）為一個橢圓形的偏微分方程，穩(wěn)態(tài)熱傳導(dǎo)問題本質(zhì)上可以看作瞬態(tài)熱傳導(dǎo)問題的特殊情況，即時間無限長且溫度恒定時的情況。這類方程可以采用多種方法進(jìn)行求解，過程較瞬態(tài)要簡單得多，常用方法是通過求解一系列線性方程組來逼近真實解，工程中應(yīng)用較多的是Newton-Raphson法和修正Newton-Raphson法。

綜上，有限元-差分法綜合了兩種計算方法，適用性更強，可用于求解非線性偏微分方程，即瞬態(tài)溫度場分布情況，可以根據(jù)計算機配置情況及工程需要合理選擇計算精度，可用于求解高溫混凝土內(nèi)部溫度場。但這種計算方法也是具有一定的局限性的，很難同時保證穩(wěn)定性和精度要求。運用兩點后差分方法進(jìn)行計算時能夠保證計算過程的穩(wěn)定性，可以避免無效振蕩，但計算結(jié)果精度不高。Crank-Nicolson方法、Galerkin方法以及三點后差分法可以在滿足穩(wěn)定性的同時提高計算精度，缺點是對單元尺寸和時間步長的選擇要求較高，求解過程容易發(fā)生振蕩，為了克服無效振蕩，通常需要進(jìn)行變步長運算，過程較復(fù)雜。

4 結(jié)論

運用理論計算研究了高溫狀態(tài)下混凝土內(nèi)部的溫度場分布情況，總結(jié)了火災(zāi)過程中的升溫曲線變化規(guī)律，推導(dǎo)了熱傳導(dǎo)微分方程，比較了求解的三種方法，得出以下結(jié)論:1）升溫曲線為對數(shù)函數(shù)，初期升溫速率較快，后期趨于平緩；2）桿系結(jié)構(gòu)可簡化為二維溫度場進(jìn)行計算，平面結(jié)構(gòu)可簡化為一維溫度場；3）實際工程中通常將各復(fù)雜情況考慮到綜合換熱系數(shù)中進(jìn)行計算，通過軟件擬合，得到了綜合換熱系數(shù)和火焰溫度的關(guān)系曲線；4）對比三種求解方法，發(fā)現(xiàn)有限元-差分法適用各種復(fù)雜工況條件，可用于高溫狀態(tài)下混凝土內(nèi)部溫度場計算。