非均匀日照条件下结构的三维温度场分析

2018-08-20朱宏平

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

高飞，陈潘，翁顺，朱宏平

(华中科技大学土木工程与力学学院，湖北武汉 430074)

随着社会经济、工程材料的迅猛发展，高层建筑日益增多。日照辐射对处于施工期间的高层结构有显著影响[1]。暴露在自然环境中的结构会受到周期性的日照辐射作用，由于结构受局部日照和混凝土导热性较差，结构整体产生不均匀温度场，这会影响结构的水平位移并使结构内部产生较大应力。

结构温度的变化与多种环境因素有关，主要包括太阳辐射、大气温度和风速[2～4]。混凝土温度从外表面到内部呈指数递减分布[5]，具有明显的不均匀性。已有用于温度分析的有限元模型主要是二维平面模型[6]，很难准确模拟实际的热边界条件。本文建立温度分析的三维实体有限元模型，采用既有结构自由度又有热自由度的结构-热耦合单元。通过引入太阳辐射历系统建立虚拟太阳，利用Ansys的辐射矩阵生成器选取日照单元[7]，实现了复杂的日照辐射热荷载的施加，这对于建筑结构温度效应的研究有重要意义。

1 太阳相对位置的确定

太阳和地球的相对位置以年为周期变化。太阳相对于地面上任一物体的位置可以由三个参数确定：太阳高度角β、太阳方位角α和日地距离d(图1)，由于日地距离d在一年内的相对变化值只有3.3%，可被视为常数。

图1 太阳相对位置

传统的用于计算太阳高度角和太阳方位角的公式考虑了地方时和标准时的时差，但实际上它对物体表面的日影分布没有影响[8]。这里引入一种新的计时系统-辐射历计时系统，用于简化传统的太阳高度角和太阳方位角计算公式。

辐射历计时系统用本初子午线绕地轴转过的角度——辐射历时角表示，规定3月21号(春分)0∶00时刻为计时零点，每经过1 h辐射历时角增加15°，因此一年内辐射历时角的最大值为15×24×365=131400°，任意一天的任一时刻有唯一的辐射历时角τ与之对应。例如，3月22号6∶30时刻，对应的τ=1×24×15+6.5×15=457.5°。

通过引入辐射历时角τ，太阳高度角和太阳方位角可以由下列式子表示：

(1)

α=

(2)

式中：k=0,1,2，…，364;A=π/180，表示由弧度制向角度制转化的乘数因子；φ表示纬度；

2 热分析理论

混凝土结构表面温度与其所处地理位置、气候条件、太阳辐射、表面朝向、大气温度、风速等有关。对于一确定结构，其表面热交换主要有三种方式：太阳辐射、热对流和长波辐射。本文分析将用到下列假设条件：

(1)混凝土材料的热性能不因温度的变化而变化；

(2)混凝土是各向同性、线弹性材料；

(3)热对流系数和风速线性相关。

2.1 热传导方程

热传导指完全接触的两个物体之间或同一物体的不同部分之间由于温度梯度而引起的能量交换过程。由著名的傅里叶热传导方程表示为[9]：

(3)

式中：ρ为混凝土密度(kg/m3)；c为混凝土比热(J/(kg·K))；T为结构温度(℃)；kx,ky,kz分别表示x，y，z方向的导热系数(W/(m·K))。

2.2 初始条件和边界条件

热传导方程建立了温度、时间和空间之间的关系。初始条件和边界条件是一个具体问题能求得定解的必要条件。一般来说，日出前1～2 h内结构整体温度分布较为均匀，为简化起见，这里取凌晨5∶00时刻，结构初始温度取为同一值，即

T(x,y,z,0)=T0=常数

(4)

实际上，在初始时刻，结构整体温度分布不为同一值，在后续分析中，将前一天的最后时刻温度场作为初始时刻温度场，即

Tn(x,y,z,tinitial)=Tn-1(x,y,z,tlast)

(5)

对于瞬态热分析，通常有如下三种边界条件。

第一类边界条件，指物体边界上的温度函数已知，即

Ts=Tb(P,t),P∈s,t>t0

(6)

式中：s表示外边界；P=P(x,y,z)，表示空间点的坐标；t0为初始时刻。

第二类边界条件，指物体边界上的热流密度函数已知，即

(7)

式中：n为面法向；qb表示边界上的热流密度函数。

第三类边界条件，指物体边界上的热对流函数已知，即

式中：hc表示结构表面和大气温度之间的热对流系数(W/(m2·K))；Ta表示大气温度(℃)，它可以用正弦函数表示为：

(9)

太阳辐射强度可以用热流密度来衡量，表示单位时间内，通过物体单位横截面上的热量。太阳照射在结构表面上，相当于在结构表面施加一定数量的热流密度值。在有限元分析中，很难同时施加第二类和第三类热边界条件，所以将第二类边界条件转化为第三类边界条件，即综合边界条件可以表示为：

式中：hr为辐射热交换系数(W/(m2·K))；qs为结构外表面获得的总的太阳辐射热流密度值，由式(11)确定；qra为长波辐射热流密度值，由式(12)计算得到；k为导热系数(W/(m·K))。

qs=αs(IDf+IDβ+If)

(11)

(12)

式中：αs为结构外表面的短波辐射吸收系数，通常取0.50～0.77；IDf，IDβ和If分别表示太阳直接辐射热流密度值、太阳散射辐射热流密度值和地面反射辐射热流密度值(W/m2)；γ为表面和水平地面之间的倾斜角；εa为长波辐射吸收系数，通常取为0.82；ε表示材料的辐射发射率，取0.85～0.95；C为斯蒂芬-玻尔兹曼常数，C=5.67×10-8W/(m2·K4)。

3 有限元模型

本文以在建的长江航运中心工程为背景，建立三维有限元实体模型，研究日照辐射对结构的影响规律。武汉长江航运中心(图2)位于汉口民权路与民生路交汇处，主体结构地下4层，地上62层，建筑高度330 m，为框架-核心筒结构。建筑整体呈正方形，平面尺寸约为52.9 m×52.9 m，结构由外围20根型钢巨型柱和中间核心筒组成，核心筒和巨型柱通过混凝土梁连接。

图2 武汉长江航运中心

ANSYS具有强大的耦合场分析功能[10]。长江航运中心有限元实体模型(图3)选用Solid226单元。Solid226单元是三维热-结构耦合单元，每个节点具有ux，uy，uz和temp四个自由度。在建立有限元实体模型时，忽略次要构件的影响，主要研究柱、楼板和核心筒的温度分布。在划分网格时，对不同的构件采用不同的网格尺寸(图4)。柱单元尺寸最小，核心筒单元尺寸次之，楼板单元尺寸最大。

图3 全楼有限元模型图4 网格划分(局部)

热分析的难点是如何将真实的太阳辐射边界条件施加到有限元模型上去。实际上，太阳光沿直线投射到结构表面，结构表面有太阳光照射的区域称为阳面，没有太阳光照射的区域称为阴面，一天中随着太阳的移动，阳面和阴面不断变化。利用ANSYS中的辐射矩阵生成器可以实现阳面和阴面的区分，辐射矩阵生成器用于求解空间中面与面之间的形状系数，如果形状系数不为零，表明这两个面之间可以相互“看见”对方，反之，表明两个面之间有遮挡。

首先利用辐射历计时系统，在ANSYS中建立虚拟太阳(图5，6)。虚拟太阳为一个立方体单元，单元大小与分析模型的最小单元尺寸相当，为0.6 m×0.6 m×0.06 m，它与结构模型中心的距离为500倍的最小单元尺寸。这样可以更真实地模拟平行太阳光。根据式(1)，(2)计算得到一天中任一时刻太阳的相对位置。

图5 虚拟太阳图6 虚拟太阳运动轨迹

将虚拟太阳单元的正对面与分析模型的所有面组成辐射对，然后进入Ansys角系数矩阵生成器(/AUX12)，利用VFCALC命令计算得到虚拟太阳单元与结构单元的角系数矩阵，矩阵的维数为1×n，即表示虚拟太阳-结构单元1、虚拟太阳-结构单元2、……、虚拟太阳-结构单元n(假设模型结构有n个单元)的角系数F11，F12，…，F1n。如果F1i=0，则表示虚拟太阳“看不见”i结构单元，即此单元为阴面单元；如果F1j≠0，则表示虚拟太阳“看得见”j结构单元，即此单元为阳面单元。最后通过条件语句实现阴面单元和阳面单元的区分。图7，8中紫色部分为选取的阳面单元。

图7 9∶00时刻阳面单元图8 12∶00时刻阳面单元

选取2017年3月6号这一天(晴朗天气)，每隔一小时输入一次边界条件，共计分析24 h的温度场。当太阳位于地表面以上时，结构得到照射，即太阳高度角大于0时，才计算阳面单元，这样能大大节省计算时间。为了清晰表达分析结果，这里仅截取模型0～50 m的高度范围。

气象站可查询到武汉市2017年3月6号的气象资料，大气温度、太阳辐射热流密度值和风速的变化如表1所示。混凝土参数取值如表2所示。

表1 气象数据

表2 混凝土参数取值

图9所示为3月6号不同时刻的温度场分布图。6∶30时刻模型东侧温度开始升高，随着太阳的移动，南侧和西侧温度也慢慢升高；在15∶30时刻，西南侧温度达到最大，约为21 ℃；在18∶30之后，结构温度缓慢降低，整体温度场趋于均匀，在计算终止时刻，模型整体温差约为3.8 ℃，为非均匀分布，在分析完一天的温度场时，将终止时刻的温度分布作为后一天的初始温度场。阳面随着太阳的转动而移动，并且阳面温度要高于阴面的温度，真实地模拟了太阳的动态变化过程。