国家体育场日照非均匀温度效应模拟

2022-10-16沈文爱万华平

土木工程与管理学报 2022年4期

马强，沈文爱，陈潘，万华平

(1. 浙江大学建筑工程学院，浙江杭州 310058；2. 华中科技大学土木与水利工程学院，湖北武汉 430074)

随着我国经济的快速增长，大跨空间钢结构向更大跨度、更复杂造型等方向发展，其温度效应也更为显著。现行规范采用50年重现期的月平均最高气温Tmax和月平均最低气温Tmin，并结合合拢温度，来计算结构均匀升温或降温时的温度效应[1]。由于大跨空间钢结构的温度场呈现出非均匀性和时变性的特点[2]，基于均匀温度场来计算大跨空间钢结构温度效应的做法有待商榷。导致非均匀和时变温度场的主要原因包括：(1)大跨空间钢结构常采用透光性较好的屋面材料(如ETFE膜材)；(2)部分构件直接暴露在外，结构主体受到太阳直射；(3)对流换热等多种因素耦合作用。非均匀和时变温度场不可避免会导致结构产生非均匀温度变形与温度应力，有必要对大跨空间钢结构日照非均匀温度效应开展研究。

目前，构件温度场模拟可以利用ANSYS软件热分析模块中的辐射矩阵生成器(AUX12)或基于空间解析几何关系来考虑遮挡作用，实现构件的瞬态热分析[3]。其中，H型截面构件[4]、圆形截面构件[5]及箱型截面构件[6]日照非均匀温度场的有限元模拟结果和实测值近似一致。由于构件的温度场模拟采用SOLID单元，而大跨空间结构有限元模型大多数采用梁单元建立，导致构件的研究结果难以在大跨空间结构中直接应用。

近年来，大跨空间结构日照非均匀温度场模拟成为了研究热点，较多学者在这方面开展了相关研究。范重等[7]对鸟巢钢结构太阳辐射照度分区，以此来确定构件的太阳辐射温度。刘红波等[8]将箱型截面构件的温度场分布应用到大跨度钢拱结构中。赵中伟等[9]对天津市于家堡交通枢纽站的网状穹顶结构日照非均匀温度场进行了数值模拟。上述研究均未考虑构件间的遮挡，这对于造型复杂、构件之间遮挡较为严重的大跨空间结构，是不可忽略的。鉴于此，陈德珅等[10]采用介于0和1之间的光照系数来考虑构件间的遮挡影响，表达结构的阴影分布情况。周勐等[11]以北京大兴国际机场钢结构为背景，基于消隐算法，实现了圆截面杆系结构日照阴影的计算。虽然大跨空间结构日照非均匀温度场模拟取得了较好的进展，但是考虑遮挡关系的复杂大跨度空间结构温度场及其温度效应的模拟方法还需进一步深入研究。

本文编制了钢板辐射升温的计算程序，并通过太阳与构件的空间解析几何关系，来判断遮挡情况。针对遮挡情况复杂的大跨空间结构，本文忽略了单元内部温度分布的非均匀性，使得结构整体非均匀温度场的模拟过程更高效和实用。以国家体育场主体钢结构为研究对象，通过对结构分区和遮挡判断的参数分析，实现了结构日照非均匀温度场的模拟和基于有限元模型的结构稳态响应计算。本文方法为大跨空间结构的日照非均匀温度效应分析提供了一种有效途径。

1 钢板温度场模型

1.1 热流类型

暴露在空气中的钢板，在太阳辐射的作用下，存在三种热运动，分别是热辐射、热对流、热传导，如图1所示。

图1 钢板的热流类型

1.2 太阳辐射

本文采用美国供暖、制冷和空气调节工程师协会的ASHRAE晴空模型计算太阳辐射强度。在该模型中，假设太阳总辐射为直接辐射、散射辐射和从周围地表反射的太阳辐射之和。

(1)太阳直接辐射强度

(1)

式中：A为大气层外的太阳辐射强度；B为大气消光系数；β为太阳高度角；CN为大气清洁度；θ为入射角。

(2)太阳散射辐射强度

Gdθ=CGNDFws

(2)

(3)

式中：C为平面上散射辐射与垂直入射辐射的比值；Fws为构件对天空的角系数；α1为倾斜角度。

(3)太阳反射辐射强度

GR=GtHρgFwg

(4)

(5)

式中：GtH=GND+Gdθ；ρg为地面的反射率；Fwg为构件表面对地面的角系数。

因此，入射到非垂直表面的太阳总辐射为：

Gt=GND+Gdθ+GR

=[max(cosθ,0)+CFws+ρgFwg(sinβ+C)]GND

(6)

入射到垂直表面的太阳总辐射为：

Gt=GND+Gdθ+GR

C)]GND

(7)

(8)

构件实际得到的太阳辐射热流密度qs为：

qs=εGt

(9)

式中：ε为构件表面的太阳辐射吸收率。

1.3 长波辐射

净长波辐射可以表示为：

(10)

式中：ζ为表面长波发生率；σ为斯蒂芬 - 玻尔兹曼常数，取为5.67×10-8W·m-2·K-4；Tg为地表温度；T为钢板温度;Tsky为有效天空温度，夏季可取为Ta-6 ℃[12]，Ta为环境空气温度。

1.4 对流换热

对流换热主要指物体表面和空气流体之间的能量交换。在建筑热分析中，对流换热系数与风速可用线性关系进行模拟，对流换热系数h按下式取值[13]：

h=6.60+3.26v

(11)

式中：v为平均风速。

1.5 热传导

假定钢材的热性能是均匀的和各向同性的，因此钢板的一般瞬态导热控制方程为：

(12)

式中：ρ为钢材密度；c为比热；k为导热系数；t为时间。

1.6 边界条件

在钢板边界表面，需要满足边界面的热平衡。

(13)

式中：λ为钢的导热系数；n为钢板表面的法线方向；qs为太阳辐射热流密度；qr为净长波辐射热流密度。

1.7 钢板温度场

在太阳辐射作用下，若假设钢板内温度均匀分布，则钢板内部导向钢板表面的热流密度为0，由边界面的热平衡方程可得：

h[T-Ta(t)]=qs(t)+qr(t)

(14)

上式是关于钢板温度的复杂函数，可迭代求解钢板温度T。

2 遮挡判断

2.1 太阳相对位置的确定

地球和太阳的相对位置是按照年为周期变化的，可以用日地距离d，太阳方位角α2，太阳高度角β三个参数来描述地球上任一点与太阳的相对位置，如图2所示。

图2 太阳相对位置

引入辐射历计时系统[3]，在该系统中，以3月21日(春分)0:00时刻为计时零点，将地球自转的角度t作为变量。

(15)

(16)

式中：k=0,1,2,…,364；A=π/180；φ表示纬度。

2.2 遮挡判断步骤

对于造型复杂、上下分层、杆件数量众多的大跨空间结构，构件间的遮挡关系较为复杂，可将大跨空间结构简化为由线段组成的集合，进而判断结构中的一点是否被其他构件所遮挡。主要计算步骤如下：

(1)确定太阳相对位置;

(2)建立全部单元的线段坐标库R;

(3)构件中的点A与太阳组成线段L;

(4)判断线段L与坐标库R的关系:若线段L与坐标库有交点，则表示点A被遮挡；反之，则未被遮挡。

3 国家体育场

3.1 工程概况

国家体育场屋面呈双曲马鞍形，长轴为332.3 m，短轴为296.4 m，总建筑面积约为25.8万m2，主体钢结构(钢材物理性能指标见表1)支撑在24根桁架柱上，呈1/2旋转对称，主结构间通过次结构实现交叉支撑，提高整体抗侧刚度的同时形成了“鸟巢”的艺术效果[14]。

表1 钢材物理性能指标

3.2 有限元模型

ANSYS软件用来建立国家体育场的有限元模型(图3)，所有构件采用BEAM4梁单元来模拟，边界条件为固结。BEAM4三维弹性梁单元有2个节点，6个自由度，温度可以作为体力在BEAM4单元的八个“角”输入。

图3 国家体育场钢结构有限元模型

4 国家体育场均匀温度效应

北京地区夏季温度很高，且湿度较大，冬季温度很低，较为干燥，是典型的温带大陆性气候。本文均匀温度作用的最大正温差取50.6 ℃，最大负温差取-45.4 ℃，合拢温度取14±4 ℃[7]。

在结构自重作用下，内环桁架的变形最为显著，故选取如图4所示的4个节点，来研究均匀温度作用对结构变形的影响。同时，选取节点2和节点5附近的八个单元，来研究均匀温度作用下的单元应力变化情况，如图5所示。轴向应力为沿杆件轴线方向的应力，而最大应力为轴向应力加上弯曲应力，即杆件横截面边缘处的应力。

图4 节点选取

图5 单元选取

在均匀温度作用下，图6为节点1～4的位移随温差变化情况，图7为单元①～⑧的应力随温差变化情况。

图6 均匀温度作用下的节点位移变化

图7 均匀温度作用下的单元应力变化

在均匀温度作用下，国家体育场钢结构的节点位移、单元应力变化均是线性的：温度升高，则结构向外膨胀；温度降低，则结构向内收缩。部分单元对均匀温度作用敏感，如单元⑦，最大负温差下的应力和自重作用下的应力处于同一数量级。不是所有构件都对均匀的温度变化敏感，如上述单元④，单元应力基本不随温度荷载的变化而变化。对均匀温度效应的研究可以用来模拟主体钢结构在无光照情况下，或冬季严寒气候下，或以全年为时间跨度的结构性能变化。

5 国家体育场日照非均匀温度效应

选取2020年夏至(即6月21日)来研究国家体育场钢结构的日照非均匀温度效应。

5.1 钢板辐射升温

(1)温度与风速数据

从WheatA小麦芽 - 农业气象大数据系统，可获取2020年夏至日北京地区大气温度Ta和风速v的变化情况，如图8，9所示。其中,地表温度Tg采用正弦曲线近似模拟,Tg的最低值Tg,min取为5:00时的气温，最高值Tg,max取为[15]：

图8 温度变化情况

图9 风速变化情况

Tg,max=8.67+0.728×Ta,max+0.0278×Gt

(17)

(2)其他参数

国家体育场钢结构面漆为金属浅银灰色，太阳辐射吸收率为0.34，涂装防护设计使用年限25年[7]。在设计使用年限内，太阳辐射吸收率不大于0.45[16]，考虑到材料的老化，本文采用线性插值，取ε=0.40。系数A，B，C和其余各参数取值如表2所示。

表2 太阳辐射参数取值汇总

选取水平放置、倾角为45°和竖直放置3种钢板来研究其一天之中的辐射升温。钢板的辐射升温还与其朝向有关，因钢板的温度要在单元遮挡类型判断后赋值给一类单元，这一类单元外表面的朝向不尽相同，故从保守的角度考虑取各个时刻钢板的朝向始终指向太阳，此时钢板接受到辐射量最多，所计算出的钢板温度是最高的。由此可以得到不同倾角的钢板一天之中的温度变化，如图10所示。

图10 不同倾角钢板温度变化

5.2 结构分区

国家体育场有限元模型共由9121个梁单元组成，现将其分为如图11所示的五部分：(1)上弦杆+顶部次结构部分，共1492个梁单元；(2)腹杆部分，共780个梁单元；(3)下弦杆部分，共596个梁单元；(4)外围次结构+桁架柱外围部分，共4669个梁单元；(5)楼梯+桁架柱内侧部分，共1584个梁单元。

图11 结构分区

5.3 热分析基本假设

(1)在太阳辐射作用下，国家体育场箱型构件不同表面温度变化趋势一致[17]，且钢材具有良好的导热性，故本文假设单元内部温度均匀分布。

(2)钢结构所有单元的遮挡类型仅分为未被遮挡的阳面单元和被遮挡的阴面单元两类。在太阳辐射的作用下，阳面单元所施加的温度采用相同条件下钢板的温度值，阴面单元则施加气温。屋面单元和下弦杆单元中的阳面单元对应平放钢板；腹杆单元中的阳面单元对应倾角为45°的钢板；外围结构中的阳面单元对应倾角为90°的钢板。

(3)本文单元所施加的温度荷载已是钢板稳态温度，同时考虑到计算的简便性，故不考虑单元之间的热传导。

(4)上弦膜材的透光率为94%[7]，故本文不考虑膜材对钢板辐射升温的影响。

5.4 单元遮挡类型判断

太阳高度角<0°时为夜晚，全部视为阴面单元。当太阳高度角>0°时，需要进行阳面单元的判断。国家体育场屋面为鞍形曲面，总体而言相对平缓，现将全部屋面单元视为阳面单元。楼梯和桁架柱内侧的单元被其他部分遮挡较为严重，将其视作阴面单元。对于外围次结构和桁架柱外围部分，则通过单元与太阳的相对位置进行判断，靠近太阳的一半视为阳面单元，背离太阳的一半视为阴面单元。

对于腹杆单元和下弦杆单元，遮挡情况较为复杂，采用2.2节所述的遮挡判断步骤。其中，太阳距离d要足够大，来模拟平行光线，徐丰[3]建议取结构几何尺寸的5～10倍，本文太阳距离d取2500 m。

随机选择100个腹杆单元和100个下弦杆单元做单元分段数n的参数分析(图12，13)，时刻选为2020年6月21日14:00。单元分段数为n，表示将单元平分为n份，并取中间的n-1个点分别与太阳组成线段L。

图12 单元分段数n参数分析(腹杆)

图13 单元分段数n参数分析(下弦杆)

由图12，13可知，随着单元分段数n的增加，0个交点的单元数量逐渐减少，但在分段数n<50时，0个交点的单元数量减少得较快，在分段数n>50时，0个交点的单元数量减少的趋势变缓。随着单元分段数n的增加，有交点的单元数量增加，同时交点的个数也在增加，因此将有无交点作为阳面单元和阴面单元的判断依据略有失偏颇，需要同时考虑单元分段数的影响。本文采用单元分段数n=50时，有无交点作为阴面单元和阳面单元的判断依据。

国家体育场钢结构所有单元在夏至日各个时刻的遮挡类型均可以通过遮挡判断算法予以区分，由此可建立结构的日照非均匀温度场。图14为夏至日14:00钢结构阳面单元与阴面单元的分布情况。

图14 14:00钢结构阳面单元与阴面单元的分布

5.5 日照非均匀温度作用下的结构响应

在日照非均匀温度作用下，结构响应如图15，16所示。由图可知，节点位移和单元应力的最值出现在7:00和13:00左右。钢结构对称，若温度分布对称，则节点1和节点2、节点3和节点4的节点位移变化趋势是一致的。由图16中节点1和节点2、节点3和节点4的节点位移变化差异可见温度分布的非对称性。同样，应力水平变化应一致的单元，如单元⑥和单元⑦，在日照非均匀温度作用下的应力水平也存在差异。

图15 非均匀温度作用下的节点位移变化

图16 非均匀温度作用下的单元应力变化

在均匀温度作用下，节点位移和单元应力可认为是线性变化的。现考虑国家体育场钢结构在温差等于2020年夏至日气温变化幅度(11 ℃)和钢板温度变化幅度(23.5 ℃)的均匀温度作用下，四个节点的位移变化幅度，将其与本文模拟的日照非均匀温度场下的位移变化幅度进行对比，对比结果如图17所示。由图可知，日照非均匀温度作用下四个节点各个方向的位移变化幅度均大于气温温差(11 ℃)情况下的位移变化幅度，较为接近钢板温差(23.5 ℃)情况下的位移变化幅度。

图17 节点位移变化幅度对比

同样地，单元应力变化幅度对比结果如图18所示。由图可知，八个单元的轴向应力变化幅度和最大应力变化幅度均高于均匀温度作用下的应力变化幅度，单元⑥～⑧在非均匀温度场下的应力变化幅度为气温温差(11 ℃)作用下的3倍左右，钢板温差(23.5 ℃)作用下的2倍左右。下弦单元⑥～⑧比上弦单元①～③的应力变化幅度更大，其中，单元⑥的最大应力变化幅度接近30 MPa。