基于大涡模拟的平屋盖锥形涡数值分析研究*

李秋胜, 刘 顺

(1.湖南大学 土木工程学院，湖南 长沙 410082；2.香港城市大学，香港 100013)

基于大涡模拟的平屋盖锥形涡数值分析研究*

李秋胜1,2†, 刘 顺1

(1.湖南大学 土木工程学院，湖南 长沙 410082；2.香港城市大学，香港 100013)

采用大涡模拟(LES)对平屋盖建筑受45°风向角作用下的表面风荷载问题进行了非稳态数值模拟分析.通过与风洞试验结果的对比得出，大涡模拟能较好地捕捉到建筑物顶面出现的锥形涡及其特性.在此基础上，研究了锥形涡作用下建筑物顶面平均风压与脉动风压的分布，以及加设分隔挡板和不同高度的女儿墙对屋面风压分布和旋涡强度的影响.研究结果表明，基于Q准则的旋涡判别法可以较好地识别斜风向下屋面形成的锥形旋涡；在背风区锥形涡与侧面脱体涡相互作用并脱落，其影响将反馈至屋面旋涡上导致屋盖两个锥形涡强度以屋面对角线为轴交替波动，此消彼长；屋面女儿墙的存在使得两个锥形涡之间的间隙变窄,旋涡足迹变阔，且屋面峰值吸力随女儿墙高度的增加而迅速减小.

大涡模拟；锥形涡；Q准则；交替波动；女儿墙

当风以一定的角度吹过建筑物表面时，会产生复杂的流动结构，包括在屋顶角部的锥形涡和在背风侧面的脱体涡等.在建筑物屋面产生的锥形涡会产生很大的负压区，从而使建筑物屋面等部位承受很大的压差力.相关实测[1-2]及风洞试验[3]研究表明，低矮房屋在屋盖迎风角部和迎风前缘会遭受强风吸力作用.而风灾调查[4]也显示，强风造成的房屋破坏主要集中在低矮房屋的屋面角部、屋檐边缘和屋脊等部位.综上可知，锥形涡的存在是强风地区建筑物受破坏的主要原因之一.

考虑到建筑物屋顶锥形涡的重要性，国内外很多学者基于风洞实验对锥形涡进行了研究.Kawai[5]利用速度测量得出了建筑物顶部锥形涡的具体结构(45°风向角下)，发现均匀层流下的锥形涡强度强于湍流下锥形涡的强度，两个锥形涡交替生成、耗散引起了表面压力沿对角线不对称的脉动；Banks等[6]通过风洞试验和对TTU建筑的现场实测，运用流场可视化技术研究发现在均匀层流作用下涡核处最大吸力的大小与锥形涡的大小成反比，而对于湍流作用下的屋顶最大吸力与锥形涡的大小并没有类似的关系；Kawai[7]通过风洞试验指出屋面局部负压峰值的出现和在一定风向角下屋面形成巨大强烈的锥形涡有关，并分析了在湍流作用下产生局部负压峰值的条件，同时还探讨了改变屋檐结构形状来减少负压峰值的方法.国内方面，陈学锐等[8]通过风洞试验研究了在锥形涡诱导下建筑物顶面风荷载的特性，给出了在不同风向角下压力分布的结果，分析了产生的原因和流动机理以及建筑物顶面的分离流动结构，并指出锥形涡的出现是建筑物顶面局部出现峰值负压的主要原因.

相对于诸多锥形涡的风洞试验研究，有关锥形涡的数值模拟研究较少，李鹏年等[9]以及陈青松[10]利用流体力学计算软件FLUENT，选择v2-f湍流模型对40°风向角下建筑物顶面锥形涡的演化、强度和位置与建筑物表面压力分布进行了分析.此外还模拟分析了风向角和建筑物高度对屋顶锥形涡的影响.随着计算机技术、数值计算和湍流模拟技术的发展，采用数值方法对建筑物绕流进行数值模拟更为简捷、经济，同时还可以得到某些风洞实验和现场实测不能观测到的结果.本文利用数值模拟的优势，对长宽高比为1∶1∶0.5的建筑模型进行了大涡模拟(LES)研究，通过数值的可视化处理模拟了以Q准则识别的锥形涡结构，并重点分析了45°风向角下屋顶两锥形涡强度的非稳定的波动特性.此外，本文也进行了有女儿墙的平屋盖模型的大涡模拟，探讨了女儿墙的存在及高度变化对屋面风荷载分布的影响以及对屋面锥形涡的结构和其他特性的影响.

1 数值风洞

1.1 大涡模拟方法

本文选用大涡模拟进行CFD数值计算[11-12]，其原理是将流动中的旋涡分成大涡和小涡，对大涡进行直接求解，对小涡采用亚格子尺度模型进行计算.大多数亚格子模型都是在涡粘性的基础上，把脉动的影响用一个湍流粘性系数μt来表示.根据各亚格子模型的特点，本文采用一方程亚格子模型来求解.

(1)

亚格子涡粘性系数μt可由亚格子湍动能表示为：

μt=Ckksgs1/2Δf．

(2)

式中:Δf表示过滤的尺度，由Δf=V1/3确定.

一方程亚格子涡粘模型可以考虑亚格子动能的历史影响，通过其生成项，耗散项以及扩散项得以体现.如下式：

(3)

1.2 数值计算模型

本文选用长宽高之比为2h∶2h∶h(h=60 mm)的模型进行数值模拟，计算域、模型位置及边界条件的设置如图1所示.来流风向角为45°，来流攻角取零度.

(a)

(b)

本文计算流域网格划分采用Hexcore型非结构化网格，由TGrid网格软件划分而成.建筑模型表面及计算域地面附加边界层，以便更加准确地模拟近壁面区的流动.通过多次试算，最终确定的网格最小尺度为0.000 5h，网格总数在90万左右.经计算壁面网格无量纲高度y+≤3(y+=ρuy/μ).对LES计算来说，近壁面网格的疏密对于模拟的计算结果影响相对较大，近壁面网格越密对壁面流动的描述越好，而亚格子模型的影响相对较小.LES湍流模型要对壁面边界层进行完全求解，网格要求是y+≈1.本文的y+虽略大于1，但通过采用增强型壁面函数，LES湍流模型的结果能满足壁面湍流的处理要求,可保证结果的可靠性.

本文模拟均匀流场，入口切向速度为零，只有法向速度.为了方便与前人的风洞试验结果做对比，入口速度设为15 m/s(建筑雷诺数约为1.2×105)，保证了与风洞试验的雷诺数相一致，可以避免雷诺数对结果的影响.因为本文重点不在分析雷诺数对锥形涡的影响，故暂未考虑不同雷诺数的变化.为了使计算更快更稳定的收敛，在大涡模拟计算之前先进行了RANS模型的计算，将RANS模型计算的结果通过瞬态化处理作为大涡模拟计算的初始流场.至于入口速度的脉动成分则采用Fluent中的Spectral Synthesizer[13]法生成.为了切实地模拟风场，在入口上加入少许湍流度(Iu=0.5%).建筑物表面采用无滑移壁面，地面采用自由滑移壁面，对压力和速度场的耦合采用SIMPLEC算法求解，流体的空间离散采用二阶迎风格式(Second Order Upwind)，时间步长经试算对比，取为0.001 s，配合一方程亚格子模型进行模拟.通过多次试算并和已有风洞试验结果对比发现，本文所用方法能较好地模拟出屋面平均风压及脉动风压的分布.

1.3 分析工况

本文共设置了六组工况，为了方便对比，又将各工况归纳为3类：

A类：无分隔板，无女儿墙(设为A0工况)；

B类：有分隔板，无女儿墙(详见图2)；

C类：无分隔板，带女儿墙(详见图3).

其中，为了避免分隔板厚度对数值模拟的结果造成影响，对分隔板进行了零厚度处理.对C类工况，有C1，C2和C3组工况，对应女儿墙的高度h0分别设为0.05h，0.1h，0.3h(h为模型的高度).

(a)B1：屋面沿对角线设竖直分隔板(高h)

(b)B2：尾部沿流向设竖直分隔板(长6 h)

图3 C类设女儿墙工况(女儿墙高h0)

为准确地追踪到锥形涡的实时特性，根据试算结果和前人的分析研究，特在锥形涡范围内与迎风边沿的夹角约为θ=14°的角线上依次布置了一系列监测点，各监测点沿屋面对角线对称布置，如图4所示.

图4 屋面监测点布置

2 计算结果及分析

2.1 涡结构可视化处理

本文依据Q准则来判别漩涡区域[14]，从而识别锥形涡及其他旋涡结构，达到可视化的目的.Q准则是由Hunt等[15]于1988年提出，他们定义流场中速度梯度张量V的第二矩阵不变量Q具有正值的区域为旋涡.另外，它要求旋涡区域的压强要低于周围的压强.对于不可压缩流动Q可定义为：

Q=(‖Ω‖2-‖S‖2)/2．

(4)

(5)

(6)

S和Ω分别代表了流场中一点的变形和旋转.因此，Q准则反映了流场中一个流体微团旋转和变形之间的一种平衡.Q>0则反映了旋转在流动中占据统治地位.通过对模拟结果数据的转化，可得各类工况的旋涡结构，如图5所示.

由于C3工况模型附近旋涡结构图与C2工况的分布类似，仅锥形涡在屋面分布范围更大，尾部涡流更复杂，故在此省略未画出.对比上述各工况的旋涡结构图，可看出A0和B1工况各旋涡结构的分布基本相同，说明沿对角线布置的竖向分隔板对屋面锥形涡的尺寸及分布的影响较小，可忽略.而B2,C2工况(包括C3工况)的涡结构与无分隔板、无女儿墙的A0工况相比差别较大，所造成的影响不可忽略.具体来说：B2工况中的竖直分隔板影响了尾流中两侧旋涡的交替脱落及涡流间的相互作用，使得屋面锥形涡尾部涡流有所聚集；而通过比较A0,C1,C2工况中的锥形涡结构，可发现女儿墙的存在抬高了锥形涡的位置，扩大了屋面锥形涡的作用范围，且随着女儿墙高度的增加锥形涡的尺寸也随之增大，并逐渐覆盖整个屋面.

(a)A0工况模型附近旋涡结构图

(b)B1工况模型附近旋涡结构图

(c)B2工况模型附近旋涡结构图

(d)C1工况模型附近旋涡结构图

(e)C2工况模型附近旋涡结构图

2.2 风压分布特性与分析

为便于对比分析，本文将无量纲的平均风压系数及脉动风压系数分别定义为：

(7)

Cp′=2P′/(ρUH2)．

(8)

式中:Pm是指平均风压；P′为脉动风压的均方根值；UH取来流风速，为15 m/s；ρ=1.225 kg/m3.

日本建筑协会(AIJ)、东京工艺大学(TPU)及多位学者都进行了长∶宽∶高=1∶1∶0.5的平屋盖建筑模型的风洞试验.本文以Nishimura和Kawai[16]的风洞试验结果作为对比(对应本文中的A0,B1,B2工况)，以验证本文数值模拟结果的可靠性.具体的试验结果和模拟结果分别如图6及图7所示.

图6 屋面平均(上)及脉动(下)

图7 屋面平均(上)及脉动(下)风压系数模拟结果

通过对比可发现，本文相应工况的数值模拟结果与风洞试验结果吻合较好，且风压分布规律基本一致.图8分别给出了A0，B1及B2工况中的屋面最小平均风压系数及最大脉动风压系数CFD模拟结果与风洞试验结果的对比.可发现最大相对误差(为6.25%)发生在B1工况最小平均风压系数的对比上，误差相对较小且各工况的模拟结果与试验结果在整体趋势上是一致的.这也再次证明本文采用的大涡模拟能够较好地反映锥形涡下屋面平均及脉动风压的分布特性，同时也说明了本文其他工况及分析结果的可靠性.

(a)最小平均风压系数的比较

(b)最大脉动风压系数的比较

由A0,B1,B2工况的风压系数分布云图可知，A0和B1工况屋面风压分布基本一致，沿屋面对角设分隔挡板后对锥形涡的尺寸、形状及平均、脉动风压分布产生的影响较小，最小平均风压系数及最大脉动风压系数分别存在0.05和0.01的差值.在模型尾流区加设竖直分隔板后(B2工况)，最小平均风压系数下降了0.05，而屋面上的平均风压系数整体上是有所提高的，即屋面吸力有所减小.同时可看出B2工况中的脉动风压系数的大小较A0工况整体上是减小的，且相应区域最大有0.05的降幅，不可忽略.初步推断这是尾流区加设的竖向分隔板阻碍了背风区锥形涡与侧面脱体涡相互作用后的脱落及两侧涡流的相互影响所造成的.

对比A0工况与加设女儿墙的各工况(C1,C2,C3)的风压系数分布云图可明显发现：随着女儿墙高度的增加，无论是平均风压还是脉动风压，其绝对值都是迅速减小的.当女儿墙的高度达到h0=0.3h时，最小平均风压系数和最大脉动风压系数分别为-0.7和0.17，且屋面风压分布趋于均匀化.此外，C1,C2工况在背风屋角处都出现了正风压区，而当女儿墙达到一定高度时，屋面正风压区已不再存在，正如C3工况屋面的风压分布所示.联系上文对应工况的旋涡结构图，可推断上述现象是女儿墙的存在影响了气流的分离，阻碍了锥形涡尾部与侧面脱体涡的相互作用，抬高了屋面锥形涡的位置，扩大了锥形涡的范围，使得两个锥形涡的间隙变窄所致.

2.3 风压时程特性研究

模拟过程中，在屋面两锥形涡的范围内对称布置了12个监测点(详见图4)，以便观测屋面吸力较大区域的风压随时间变化的特性.同样，将瞬时风压系数定义为:

(9)

式中:P(t)为屋面瞬时风压，其他参数的意义同式(7)和式(8).

本文截取各工况屋面迎风前缘的点P1,Q1或P3,Q3较稳定的6 s时程数据进行分析，具体结果如图9所示.

时间/s

时间/s

时间/s

时间/s

时间/s

时间/s

时间/s

从图中A0和B1工况下点P3,Q3的风压系数时程曲线的变化可显著地观察到屋面锥形涡强度随时间的交替波动现象，如在2.4 s左右时刻(图中竖向箭线处)沿屋角线对称的P3点与Q3点的风压系数差值分别达到了0.4和0.5，可见这种此消彼长波动的现象是十分明显的.同时也可看出沿屋面对角设置的竖向分隔板对风压系数的数值及波动变化影响不大，仅屋面迎风尖角附近的监测点P1，Q1旋涡强度大小的交替波动现象有所加强．这可能是由于竖向分隔板的存在更有利于迎风尖角处气流的分离，从而缩短了形成锥形涡的气流过渡区.由上述分析可知，屋面两个锥形涡彼此的联系并不是直接在屋面上通过涡流的相互作用而建立的.比较A0和B2工况中P3，Q3点的风压系数时程曲线，可发现在尾流区加设竖向分隔板后旋涡强度大小的交替波动现象已基本消失，这说明屋面两锥形涡之间的联系已被切断.而B2工况所设分隔板阻挡了尾流中漩涡的脱落与相互作用，可知屋面锥形涡的强度变化与模型侧面脱体涡的相互作用和脱落有着紧密的联系.

再对比A0工况和设女儿墙的各工况(C1，C2，C3)的风压时程曲线，在大致2 s的时刻，C1工况中点P3，Q3的风压系数的差值约为0.2，而C2,C3工况中都接近为0．可见随着女儿墙高度的增加，沿屋角线对称的监测点的风压系数时程曲线趋于一致，此消彼长的波动现象也逐渐消失.究其原因，联系前文分析可能是随着女儿墙高度的增加，锥形涡范围扩大，使得监测点间的正相关性加大，同时屋面四周封闭的女儿墙阻隔了顶面旋涡与背风区侧面分离涡的相互作用.其他监测点的风压时程规律与上述规律相类似，限于篇幅原因此处不再赘述.

为了更好地了解屋面锥形涡之间的特性，特选取屋面较大风吸力区的监测点P1～P3，Q1～Q3进行脉动风压的相关性分析，具体结果见表1.

表1 监测点脉动风压相关系数

Tab.1 Correlation coefficient of fluctuating wind pressures for monitoring points

监测点P1～Q1P1～P2P1～P3P2～Q2P2～P3P2～Q3P3～Q3工况A0-0.210.240.13-0.740.84-0.72-0.72工况B1-0.590.880.88-0.770.94-0.76-0.75工况B20.240.270.350.500.180.13-0.11工况C10.950.780.57-0.330.74-0.32-0.30工况C20.970.880.820.770.900.670.54工况C30.980.960.900.920.960.880.82

从表中的数值变化可看出：工况B1中屋面上的点P2与Q2，点P3与Q3的脉动风压相关系数与工况A0较为相近，且均为负值.而工况B1中点P1与Q1的脉动风压相关系数的绝对值较工况A0有明显的增加，这是因为工况A0的迎风屋角附近锥形涡还没完全生成，沿屋角线加设分隔板后更有利于锥形涡的形成.这也再次验证了屋面上两个锥形涡之间并没有直接的相互作用或相互的影响很小，同时也说明了锥形涡强度是以屋面对角线轴呈交替变化的.至于工况B2中相应点的脉动风压相关系数则与工况A0，B1有较明显的差异，而这种差异也符合上文关于B2工况所设分隔板切断了屋面两个锥形涡间联系的结论.对比表中工况C1，C2，C3各点之间的脉动风压相关系数的数值大小和正负可发现，随着女儿墙高度的增加，屋面各点之间脉动风压的相关性明显提高.当女儿墙的高度h0=0.3h时(C3工况)，各点之间的脉动风压相关系数大都在0.9以上，这也进一步说明随着女儿墙高度的增加，屋面锥形涡与背风侧面分离涡的作用逐渐减小，屋面两个锥形涡的强度交替增减现象将逐渐消失，而呈现出的是沿屋面对角线的对称同步变化.

3 结 论

本文重点对45°风向角均匀来流作用下建筑物顶面产生的锥形涡的现象进行了大涡模拟(LES)，研究了屋面上锥形涡的相关特性并探讨了屋面常见的女儿墙对锥形涡的影响.得到的有关锥形涡的结果如下：

1)大涡模拟(LES)能很好地捕捉到屋面的锥形涡.

2)基于Q准则的旋涡判别法能较好的识别斜风向下屋面形成的锥形旋涡结构.

3)无女儿墙的平屋盖上锥形涡并不是完全对称的，其强度也不是对称变化的，而是随时间的不断变化，主要表现为两个锥形涡强度的大小交替波动，此消彼长.

4)通过在模型顶面及尾流区加设竖直分隔挡板，发现模型顶面的锥形涡在屋面上并没有直接的相互作用和影响，而是沿流动方向锥形涡逐渐变大，到达背风区后屋面锥形涡与侧面脱体涡相互作用并随尾流的漩涡脱落反馈到屋面锥形涡的强度变化上.

5)女儿墙的存在会抬高屋面锥形涡，扩大屋面锥形涡的作用范围，阻隔顶面旋涡与背风区侧面分离涡的相互作用.当在周边设置较高女儿墙时屋面风压分布趋于均匀，且屋面峰值吸力明显减小.

[1] LI Q S, HU S Y, DAI Y M,etal. Field measurements of extreme pressures on a flat roof of a low-rise building during typhoons [J]. Journal of Wind Engineering and Industrial Aerodynamics, 2012, 111: 14-29.

[2] 李秋胜,胡尚瑜,戴益民,等. 低矮房屋屋面实测峰值风压分析[J]. 湖南大学学报:自然科学版, 2010, 37(6): 11-16.

LI Qiu-sheng, HU Shang-yu, DAI Yi-min,etal. Analysis of the field measured suction peak pressure coefficients on the flat roof of a low- rise building [J]. Journal of Hunan University: Natural Sciences, 2010, 37(6): 11-16. (In Chinese)

[3] MAHAOOD M. Experiments to study turbulence and flow past a low-rise building at oblique incidence [J]. Journal of Wind Engineering and Industrial Aerodynamics, 2011, 99: 560-572.

[4] VAN DE LINDT J W, GRAETTINGER A, GUPTA R,etal. Performance of wood-frame structures during Hurricane Katrina [J]. Journal of Performance of Constructed Facilities, 2007, 21(2): 108-116.

[5] KAWAI H. Structure of conical vortices related with suction fluctuation on a flat roof in oblique smooth and turbulence flow [J]. Journal of Wind Engineering and Industrial Aerodynamics, 1997, 69/71: 579-588.

[6] BANKS D, MERONEY R N, SARKAR P P,etal. Flow visualization of conical vortices on flat roofs with simultaneous surface pressure measurement [J]. Journal of Wind Engineering and Industrial Aerodynamics, 2000, 84: 65-85.

[7] KAWAI H. Local peak pressure and conical vortex on building [J]. Journal of Wind Engineering and Industrial Aerodynamics, 2002, 90: 251-263.

[8] 陈学锐,顾志福,李燕. 锥形涡诱导下建筑物顶面风荷载[J]. 力学学报, 2007, 39(5): 655-660.

CHEN Xue-rui, GU Zhi-fu, LI Yan. Conical vortex induced wind loading on the roof of a building [J]. Journal of Theoretical and Applied Mechanics, 2007, 39(5): 655-660. (In Chinese)

[9] 李鹏年,邹正平,陈学锐,等. 建筑物顶部锥形涡空间演化的数值模拟[J]. 航空动力学报, 2006, 21(5): 891-896.

LI Peng-nian, ZOU Zheng-ping, CHEN Xue-rui,etal. Numerical study of conical vortices space evolvement on a flat roof [J]. Journal of Aerospace Power, 2006, 21(5): 891-896. (In Chinese)

[10]陈青松. 不同高度和坡度的建筑物顶面风荷载及锥形涡研究[D]. 北京：北京大学, 2006: 13-20.

CHEN Qing-song. Study on wind load and conical vortex of building roof with different heights and sloping angles [D]. Beijing: Peking University, 2006: 13-20. (In Chinese)

[11]HUANG S H, LI Q S, XU S L. Numerical evaluation of wind effects on a tall steel building by CFD [J]. Journal of Constructional Steel Research, 2007, 63: 612-627.

[12]卢春玲,李秋胜,黄生洪,等. 大跨度屋盖风荷载的大涡模拟研究[J]. 湖南大学学报:自然科学版, 2010, 37(10): 7-12.

LU Chun-ling, LI Qiu-sheng, HUANG Sheng-hong,etal. Large eddy simulation of wind loads on long-span roofs [J]. Journal of Hunan University: Natural Sciences, 2010, 37(10): 7-12. (In Chinese)

[13]FLUENT 6.3. The user guide of FLUENT 6.3 [M]. USA: FLUENT Inc, 2006: 871-879.

[14]胡子俊,张楠,姚惠之,等. 涡判据在孔腔涡旋流动拓扑结构分析中的应用[J]. 船舶力学, 2012, 16(8): 839-846.

HU Zi-jun, ZHANG Nan, YAO Hui-zhi,etal. Vortex identification in the analysis on the topology structure of vortical flow in cavity [J]. Journal of Ship Mechanics, 2012, 16(8): 839-846. (In Chinese)

[15]JCR H, WRAY A, MOIN P. Eddies, streams, and convergence zones in turbulent flows[R]. Center for turbulence research report CTR-S88, 1988: 193-208.

[16]NISHIMURA H, KAWAI H. Switching phenomenon of conical vortices on the flat roof of a low-rise building [J]. Journal of Wind Engineering, 2010, 35(4): 99-106.

Large Eddy Simulation of the Characteristics of Conical Vortex on a Flat Roof

LI Qiu-sheng1, 2†, LIU Shun1

(1. College of Civil Engineering, Hunan Univ, Changsha, Hunan 410082, China;2. City Univ of Hong Kong, Hong Kong 100013, China)

This paper was concerned with the unsteady characteristics of the conical vortex on a flat roof on the basis of numerical results by large eddy simulation (LES) when an approaching flow attacks at about 45 degrees against the wall. The results indicate that the large eddy simulation (LES) can give fairly good results, compared with the experiment pressure data of the AIJ flat roof model and other wind tunnel tests. Based on the Q criterion, this paper presented the flow visualization of conical vortices on the roof. By analyzing the mean and fluctuating pressure coefficients on the roof, it is observed that the strong suction acts on near leading edge on the flat roof of the building. The strong suction is caused by a pair of conical vortices whose strength is unbalanced in any time. When the strength of one vortex recedes, the vortex on the other side is strengthened. By adding splitter plates and parapets, the LES results show that the switching of conical vortices may have inseparable relations with the interaction of vortices formed in the wake shed from the sidewalls of the building. Additionally, the parapets can mitigate the high corner suctions on the roof and lift the position of vortex core, while the dimension of the conical vortex is greater with the higher parapets.

large eddy simulation; conical vortex; Q criteriton; the switching phenomenon; parapets

2014-10-06

国家自然科学基金重大研究计划重点项目(90815030)；国家自然科学基金资助项目(51178179),National Natural Science Foundation of China(51178179)

李秋胜(1962-)，男，湖南永州人，湖南大学教授

†通讯联系人，E-mail:bcqsli@cityu.edu.hk

1674-2974(2015)11-0072-08

TU247.1；TU973.32

A