超大直筒−锥段型钢结构冷却塔内压效应

杜凌云，柯世堂



超大直筒−锥段型钢结构冷却塔内压效应

杜凌云，柯世堂

(南京航空航天大学 土木工程系，江苏 南京，210016)

基于计算流体力学方法(CFD)对某在建亚洲最高(189 m)超大直筒−锥段型钢结构冷却塔进行数值风洞模拟，提出直筒和锥段典型断面处表面绕流和尾迹特性，同时分析直筒和锥段内表面平均风荷载分布特性，并给出相应的体型系数拟合公式。基于有限元方法对超大直筒−锥段型钢结构冷却塔分别进行3种工况(内压作用、外压作用、内外压共同作用)下的静风响应分析，研究内压效应对超大型直筒−锥段钢结构冷却塔的主筒、加强桁架和附属桁架三部分结构风致响应的影响规律。研究结果表明：锥段内表面体型系数在迎风面与直筒段数值相近，其他部位偏离较大；内压效应对结构负向位移、主筒负压极值区的横杆轴力、附属桁架横杆轴力及上部加强桁架支撑轴力影响相对较大。

钢结构冷却塔；直筒−锥段型；数值风洞；平均风荷载；风致响应；内压效应

超大直筒−锥段型钢结构冷却塔作为一种造型新颖、材料独特的新型冷却塔结构形式，其塔筒由不同截面形式与材料等级的钢管组成复杂的桁架体系，与传统的钢筋混凝土冷却塔[1−2]相比，钢结构冷却塔阻尼更小、频率更低[3−4]，对风荷载更加敏感。考虑到塔筒顶部露天敞开且底部由支柱支撑从而形成巨大的风通道[5−6]，其内表面风荷载对冷却塔运营安全性能的影响不容忽略，尤其对于此类超大空间薄壳[7−8]直筒−锥段钢结构冷却塔来说，其内压效应更加显著。针对大型冷却塔的内外风压分布取值问题，冷却塔规范[9]给出了外压系数沿环向分布的三角函数拟合公式，但并未给出内压系数的明确说明；国内外很多学者也针对多种透风率[10]、不同高度与纬度[11]及考虑内部构件[12]等多种工况下的结构内吸力进行了全面系统研究，相关研究成果很好地指导了传统的双曲线型钢筋混凝土冷却塔结构的内压取值问题，然而，对于造型新颖的直筒−锥段型钢结构冷却塔，其内、外表面风荷载分布特性的研究较少，内压效应对风致受力性能影响的定性和定量的研究更少。鉴于此，本文作者以国内某在建塔高189 m的超大直筒−锥段型钢结构冷却塔为例，采用CFD方法[13−14]对其周围流场、表面绕流及尾迹特性进行数值风洞模拟，分析并归纳直筒和锥段内表面平均风荷载分布特性，并给出相应的体型系数三角函数拟合曲线。基于有限元理论，建立直筒−锥段钢结构冷却塔有限元计算模型，进行内压作用、外压作用、内外压共同作用3种工况下的结构静风响应计算，讨论内外表面风荷载对主筒、加强桁架及附属桁架等3部分结构的风致受力性能影响特点。

1 工程简介

该在建超大直筒−锥段型钢结构冷却塔位于中国山西，所处地貌类型为B类，基本风压取为0.45 kPa。塔体采用钢框架支承，由主筒、加强桁架和附属桁架3部分组成，主筒共18层，加强桁架共5层，分别设置在高度为32.5，67.0，107.0，148.0和189.0 m处，主筒和加强桁架均采用Q345钢材；附属桁架共30榀，钢材等级为Q235B。整体塔高189 m，进风口高度为32.5 m，直径为144.5 m，进风口高度以上结构外覆挡风钢板，表面光滑平整。冷却塔主要结构尺寸如表1所示，其整体结构示意图见图1。

表1 大型钢结构冷却塔主要结构尺寸

图1 冷却塔整体结构示意图

2 平均风荷载数值模拟

2.1 网格划分及参数设置

建立直筒−锥段钢结构冷却塔足尺物理模型，并以30%的透风率考虑百叶窗工作开启状态[10]。整个计算域顺风向长24，横风向宽15，高度方向为4(其中为塔筒底部直径，为塔高)。划分网格时将其分为局部和外围2个计算域，局部计算域内含结构物理模型，采用非结构化网格进行划分，外围区域形状规整，采用高质量的结构化网格进行划分，总体网格数量达1 640万，如图2所示。

计算区域及边界条件如图3所示。进行CFD数值模拟时选用3D单精度、分离式求解器，空气模型等效为不可压缩流体，其三维不可压缩N−S方程如式(1)和(2)所示。

图3 计算区域与边界条件示意图

对流项离散为二阶迎风差分格式，流场求解采用SIMPLEC算法[15−17]，计算模型采用剪切应力SST−两方程混合模型，在近壁面保留了原始−模型，在远离壁面处应用−模型，其方程和方程如式(3)~(11)所示。

其中：ν为涡黏系数；为涡量；1和2为混合函数。

2.2 绕流特性

图4和图5所示分别为立面和3个典型断面(锥段、直筒段、锥段与直筒交接区)的涡量分布图和压力速度叠合图。从图4和图5可知：

1) 背风区涡旋强度沿塔高先增大后减小，且峰值位于塔筒中部，顺风向涡量增值区域沿塔高变化规律则与之相反；下部锥段背风面的涡量增值区域沿横风向分布相对锥段−直筒段交接区和直筒段较大，且分离点更加靠前。

2) 随着塔筒高度的增加，背风面速度呈现先增大后减小的趋势，极值出现在塔筒中部，但其回流及涡旋相对较弱；冷却塔两侧速度场关于风轴方向基本呈对称分布，不同断面的迎风区流动分离点基本相同，背风区出现不同程度的回流及涡旋。

3) 部分气流透过百叶窗进入塔筒内部，在塔体内表面附着流动、撞击并向上爬升，由于直筒与锥段交接部位半径减小，上升气流受到阻碍，形成完整的三维气流涡旋，造成该部位内表面压力系数显著增大；气流越过交接部位后平稳上升，相应的直筒段压力系数沿高度分布变化微弱。

2.3 外表面平均风压特性

图6给出了CFD数值模拟上部直筒段与下部锥段数值模拟得到的平均风压与火工规范光滑双曲线混凝土冷却塔和建筑荷载规范圆截面构筑物风压的对比曲线。

由图6可知：1) 2种规范所得结果中，圆截面构筑物背风区风压平台较小，负压极值出现滞后且数值较大；2) 下部锥段的负压极值点和分离点对应角度与火工规范中双曲冷却塔的一致，但在背风区负压数值要明显小于规范值；3) 上部直筒段风压系数分布曲线与荷载规范中圆截面构筑物的风压系数分布曲线基本吻合。因此，本文数值模拟具有一定的有效性，可为此类直筒−锥段型冷却塔的风荷载设计取值提供参考。

2.4 内表面平均风压特性

图7(a)所示为塔筒典型断面内表面体型系数沿环向分布曲线。分别将锥段和直筒段2部分的典型断面体型系数沿高度和环向角度分别进行平均，如图7(b)所示。分析图7可知：

1) 由于下部锥段距离百叶窗较近，气流由百叶窗进入冷却塔内部，直接作用至锥段背风面下部，造成其内压急剧减小，并在240°处达到最小值−0.3；

图4 冷却塔立面及典型断面涡量图

图5 冷却塔立面及断面压力与速度流线叠合图

图6 数值模拟平均风压系数与规范值对比

2) 随着高度的增加，锥段背风区体型系数逐渐增大，由于直筒与锥段交接处受风面积突然减小，空气受到挤压，体型系数在交界处断面环向90°~270°范围内显著增大，最大增幅达75.5%；

3) 直筒段内表面体型系数分布曲线沿子午向分布趋势较为平稳，并沿环向关于风轴方向基本呈对称分布；

4)直筒段典型断面体型系数平均值沿环向分布波动较小，在−0.49~−0.51之间；锥段体型系数平均值在−80°~80°范围内与直筒段数值较接近，但在背风面相差较大。

以环向角度为变量，拟合给出上部直筒段和下部锥段及其交接部位内表面体型系数平均值计算公式：

其中：si()为锥段和直筒段及其交接部位内表面体型系数；为环向角度；a，b和c为冷却塔锥段或直筒段或其交接部位的内表面体型系数拟合公式参数，如表2所示。

表2 冷却塔锥段与直筒段及其交接部位内表面体型系数拟合公式参数取值

将内表面体型系数拟合曲线与模拟曲线进行对比，如图8所示，其中上部直筒段示意图中增加了与传统双曲冷却塔内表面体型系数[15]的对比曲线。由图8可见：拟合效果相对较好，直筒段内表面体型系数相对传统双曲冷却塔较小，沿环向波动性较弱，但沿环向和高度方向变化均不明显，近似为均匀分布，由此可认为本文数值模拟结果可信。

图8 锥段与直筒段及其交接处内表面体型系数拟合曲线与模拟曲线对比

3 动力特性分析

3.1 有限元建模

采用大型通用软件ANSYS建立直筒−锥段型钢结构冷却塔有限元模型，各部分构件均离散为空间梁单元，主筒底部斜杆与地面固接，内部加强桁架外围节点与主筒对应节点共用，附属桁架上部和下部分别与主筒和地面铰接。

钢结构冷却塔表面蒙皮不提供结构刚度，但考虑到表面蒙皮的质量会对后续模态和静风响应分析造成一定影响，因此，本文通过增大桁架构件密度方法考虑蒙皮的附加质量，将模型中钢材密度乘以1.15作为钢结构冷却塔的整体等效密度[18]。

3.2 动力特性分析

图9所示为直筒−锥段型冷却塔结构前100阶自振频率分布图。由图9可得：1) 结构基频仅为0.78 Hz，且前100阶自振频率分布密集，均在3.75 Hz以内，整体结构柔度较大；2) 频率随阶数的增长出现明显的平缓段与陡直段，且平缓段的结构振型在一定阶数范围内变化微弱，多呈现为附属桁架的环向扭转，但陡直段的振型沿阶数变化显著。图10所示为相应典型阶振型示意图。由图10可知：频率沿振型阶数变化斜率较大的区域主筒变形剧烈，在38阶左右出现明显的呼吸振型。

图9 结构前100阶自振频率图

图10 典型阶结构振型示意图

4 静风响应分析

4.1 参数选取与工况设置

文献[19]对本工程直筒−锥段型钢结构冷却塔进行了动态风效应和风振系数研究，基于研究结论结构风振系数取为1.9，此外，考虑到未来二期工程冷却塔的影响，干扰系数取1.1。外表面平均风荷载采用文献[18]中给出的拟合公式，内表面平均风荷载采用本文直筒段和锥段相应拟合公式进行计算。计算工况分为3类：工况1为仅内压作用；工况2为仅外压作用；工况3为内、外压共同作用。

基于静风响应计算结果，提取主筒、加强桁架和附属桁架典型节点的径向位移及构件内力，对比研究内压效应对超大直筒−锥段型钢结构冷却塔风致受力性能的影响。

4.2 主筒响应

图11所示为不同工况下主筒所有节点位移响应沿高度和环向角度变化的等势线图，其中，背离圆心为正，指向圆心为负。由图11可知：

1) 在结构仅施加内表面风荷载时，主筒节点径向位移关于180°呈明显的轴对称分布，迎风面内吸力较强，出现径向位移极大负值；由于背风区内吸力相对较弱，受周围塔体向内的挤压作用，使得变形沿径向朝外，径向位移由负转为正。

图11 3种工况下节点径向位移等势线图

2) 仅外表面风荷载作用下的塔筒径向位移亦沿风轴方向对称分布，位移峰值位于正迎风面及±(70°~100°)范围内的直筒段与锥段交接部位；位移次正极值沿着高度的增大分布角度范围逐渐变窄，且向背风区偏离；

3) 当内、外表面风荷载共同作用时外压起主导作用，内压效应改变了交接处的局部位移分布，但次极值的变窄现象并不明显。

分别选取主筒正迎风面与负压极值区的横杆和斜杆轴力进行不同工况下的对比分析，结果如图12所示。从图12可知：

图12 正迎风面与负压极值区横杆和斜杆轴力对比

1) 内表面风荷载使正迎风面和负压极值区的横杆均产生挤压作用，而对负压极值区斜杆的作用则相反，内压效应使正迎风面的斜杆轴力在外围有附属桁架的主筒下部产生拉力，在中部和上部为压力；

2) 内压效应减小了负压极值区的横杆轴力，但增大了其他区域杆件的响应，仅在塔筒顶部区域减小了杆件轴力。

4.3 加强桁架响应

加强桁架位于主筒内部，在环向起着约束主筒径向变形的作用，其自身的风致受力性能至关重要。图13所示为不同工况下5层加强桁架的节点径向位移对比曲线。由图13分析可知：

图13 3种工况下加强桁架节点径向位移对比曲线

1) 仅内压作用时，加强桁架节点径向位移沿环向均关于风轴对称，呈现先增大后减小的趋势，位移负极值−2.06 mm位于第2层加强桁架(即锥段与直筒交接部位)环向48°和300°位置处，正极值1.00 mm位于顶层加强桁架背风点。

2) 工况2和3下加强桁架节点径向位移分布趋势一致，仅第2层加强桁架径向位移与其余层偏离较大，其内压效应使得第2层加强桁架的位移负值有所 增大。

图14和15所示为第1，3和5层加强桁架腹杆轴力和支撑轴力在不同工况下的对比散点图。由图14和15可知：1) 仅内压作用下，加强桁架的腹杆和支撑轴力不随环向角度变化，腹杆和第1层支撑轴力极小，而第3和5层支撑轴力沿环向角度在−7.5 kN和+6.8 kN附近间隔跳跃；2) 腹杆轴力在工况2和3下沿环向角度呈螺旋分布，数值基本吻合，内压效应与耦合效应十分微弱；3) 第1层支撑轴力在工况2和3下关于风轴对称分布，分布趋势与外表面体型系数曲线相似；但负压效应可以显著减小第3层和第5层支撑轴力。

图14 3种工况下腹杆轴力对比示意图

图15 3种工况下支撑轴力对比示意图

4.4 附属桁架响应

图16所示为3种工况下附属桁架径向位移沿环向和子午向变化云图。从图16可知：

1) 内表面风荷载使得附属桁架径向位移正、负值间隔12°交替分布，分别在背风区和迎风面达到正、负极值；

2) 外表面风荷载作用下附属桁架径向位移正极值达到65 mm位于高50 m环向75°处，负向峰值为达到−95 mm出现在环向0°处；

3) 内压效应显著增大了附属桁架的负向位移，但减小了节点的正向位移。

图16 3种工况下附属桁架节点径向位移对比云图

图17所示为附属桁架迎风面和负压极值区的腹杆、横杆、内部斜杆及外部斜杆轴力对比。从图17可知：

图17 迎风面和负压极值区杆件轴力对比

1) 内表面风荷载在一定程度上增大了附属桁架腹杆迎风面的轴力，减小了负压极值区的轴力，外表面风荷载响应为主；

2) 内压效应在第1~7层横杆轴力中起主导作用，但在与主筒交接区域内压效应减弱，外压是导致轴力突增的主要因素；

3) 内、外部斜杆迎风面的轴力远大于负压极值区，内压效应使得内部和外部斜杆分别承受更大的拉力和压力；内压效应略微减弱了负压极值区的内部斜杆轴力，但在外部斜杆中占比较大，且对其受力不利。

5 结论

1) 冷却塔背风区涡旋强度均沿塔高方向先增大后减小，峰值出现在塔筒中部，而顺风向涡量增值区域沿塔高变化规律则与之相反，回流和涡旋在塔筒中部相对较弱；透过百叶窗的气流附着于塔体内表面向上流动，爬升至直筒与锥段交接部位受到阻碍形成完整的气流涡旋，越过该部位后气流平稳上升。

2) 直筒段内表面体型系数沿环向在−0.49~−0.51范围内发生较小波动；锥段内表面体型系数分布曲线在迎风面与直筒段内压数值较为接近，但受百叶窗和展宽平台的气流干扰，其负压区内表面体型系数与直筒段内压数值偏离较大；最终拟合给出了内表面体型系数沿环向分布的三角函数计算公式以供设计参考。

3) 自振频率低且分布密集，基频仅为0.78 Hz，频率沿振型阶数呈现明显的平缓段与陡直段，平缓段的结构振型多表现为附属桁架的环向扭转，而陡直段的振型以主筒剧烈变形为主，整体结构的呼吸振型出现在第36阶。

4) 在内压、外压及内外压共同作用下，主筒、加强桁架及附属桁架三部分节点径向位移均关于风轴对称分布，内压效应对各部分结构的负向位移影响相对较大；外表面风荷载通常为杆件轴力的主要控制因素，但内表面风荷载对主筒负压极值区横杆轴力、附属桁架横杆轴力及上部加强桁架支撑轴力的影响较大。

[1] NOH H C. Nonlinear behavior and ultimate load bearing capacity of reinforced concrete natural draught cooling tower shell[J]. Engineering Structures, 2006, 28(3): 399−410.

[2] KE Shitang, GE Yaojun, ZHAO Lin, et al. Wind-induced responses characteristics on super-large cooling towers[J]. Journal of Central South University of Technology, 2013, 20(11): 3216−3227.

[3] 柯世堂, 余文林. 高度200 m特大型冷却塔二维风振系数取值方法及分布规律研究[J]. 建筑结构学报, 2017, 38(10): 78−87. KE Shitang, YU Wenlin. Research for value obtained methods and distribution of two-dimensional wind vibration coefficient of super large cooling towers with 200 m height level[J]. Journal of Building Structures, 2017, 38(10): 78−87.

[4] 赵林, 展艳艳, 王志男, 等. 大型冷却塔结构风致气动和气弹效应综合研究[J]. 实验流体力学, 2017, 31(3): 1−15. ZHAO Lin, ZHAN Yanyan, WANG Zhinan, et al. Comprehensive investigation into wind-induced aerodynamic and aeroelastic effects of large cooling towers[J]. Journal of Experiments in Fluid Mechanics, 2017, 31(3): 1−15.

[5] ZHAO Lin, CHEN Xu, GE Yaojun. Investigations of adverse wind loads on a large cooling tower for the six-tower combination[J]. Applied Thermal Engineering, 2016, 105: 988−999.

[6] MAHDI I, KHOSROW B. Natural draft steel hyperbolic cooling towers: optimization and performance evaluation[J]. The Structural Design of Tall and Special Buildings, 2014, 23: 713−720.

[7] HLAWICZKA S, KORSZUN K, FUDALA J. Acidity of vapor plume from cooling tower mixed with flue gases emitted from coal-fired power plant[J]. Science of the Total Environment, 2016, 554/555: 253−258.

[8] 柯世堂, 夏逸鸣, 王法武, 等. 考虑气弹效应的超大型冷却塔脉动风压非高斯特性研究[J]. 中南大学学报(自然科学版), 2013, 44(8): 3302−3309. KE Shitang, XIA Yiming, WANG Fawu, et al. Non-Gaussian features on fluctuating wind pressure for super large cooling towers with aero-elastic effect[J]. Journal of Central South University (Science and Technology), 2013, 44(8): 3302−3309.

[9] DL/T 5339—2006, 火力发电厂水工设计规范[S]. DL/T 5339—2006, Code for hydraulic design of fossil fuel power plants[S].

[10] KE Shitang, LIANG Jun, ZHAO Lin, et al. Influence of ventilation rate on the aerodynamic interference for two IDCTs by CFD[J]. Wind and Structures, 2015, 20(3): 449−468.

[11] 沈国辉, 张陈胜, 孙炳楠, 等. 大型双曲冷却塔内表面风荷载的数值模拟[J]. 哈尔滨工业大学学报, 2011, 43(4): 104−108. SHEN Guohui, ZHANG Chensheng, SUN Bingnan, et al. Numerical simulation of wind load on inner surface of large hyperbolic cooling tower[J]. Journal of Harbin Institute of Technology, 2011, 43(4): 104−108.

[12] 董国朝, 张建仁, 蔡春声, 等. 考虑内部构件影响的超大型冷却塔内压系数研究[J]. 工程力学, 2016, 33(4): 77−83. DONG Guochao, ZHANG Jianren, CAI Chunsheng, et al. Study on internal surface pressure coefficient of super-large cooling tower with different internal main components[J]. Engineering Mechanics, 2016, 33(4): 77−83.

[13] WITTEK U, GROTE K. Substitute wind concept for elastic stability of cooling tower shells[C]// Materials & Member Behavior, Proceedings of the Sessions at Structures Congress’87. New York: ASCE, 1987, 500−513.

[14] 邹云峰, 何旭辉, 谭立新, 等. 特大型冷却塔单塔内表面风荷载三维效应及其设计取值[J]. 湖南大学学报(自然科学版), 2015, 42(1): 24−30. ZOU Yunfeng, HE Xuhui, TAN Lixin, et al. Three-dimensional effect and design values of inter surface wind loading for single super large cooling tower[J]. Journal of Hunan University (Natural Sciences), 2015, 42(1): 24−30.

[15] 周旋, 牛华伟, 陈政清, 等. 大型冷却塔风荷载干扰系数的取值方法[J]. 中南大学学报(自然科学版), 2014, 45(10): 3637−3644. ZHOU Xuan, NIU Huawei, CHEN Zhengqing, et al. Discussion on different methodologies for calculating interference factor of large cooling towers[J]. Journal of Central South University (Science and Technology), 2014, 45(10): 3637−3644.

[16] KE Shitang, GE Yaojun, ZHAO Lin, et al. Wind-induced vibration characteristics and parametric analysis of large hyperbolic cooling towers with different feature sizes[J]. Structural Engineering and Mechanics, 2015, 54(5): 891−908.

[17] 李鹏飞, 赵林, 葛耀君, 等. 超大型冷却塔风荷载特性风洞试验研究[J]. 工程力学, 2008, 25(6): 60−67.LI Pengfei, ZHAO Lin, GE Yaojun, et al. Wind tunnel investigation on wind load characteristics for super large cooling towers[J]. Engineering Mechanics, 2008, 25(6): 60−67.

[18] 柯世堂, 杜凌云, 刘东华, 等. 直筒−锥段型钢结构冷却塔平均风荷载及静风响应分析[J]. 振动与冲击, 2017, 36(7): 149−155. KE Shitang, DU Lingyun, LIU Donghua, et al. Analysis of average wind pressure and response under static wind for cylinder-conic section steel cooling tower[J]. Journal of Vibration and Shock, 2017, 36(7): 149−155.

[19] 柯世堂. 钢结构冷却塔数值风洞模拟试验及风振系数研究报告[R]. 南京: 南京航空航天大学风工程试验研究中心, 2015: 135−137. KE Shitang. Numerical wind tunnel simulation test and wind vibration coefficient for the steel cooling tower[R]. Nanjing: Nanjing University of Aeronautics and Astronautics Wind Engineering Experiment Research Center, 2015: 135−137.

(编辑 赵俊)

Effect of internal pressures for super large cylinder-conic section steel cooling towers

DU Lingyun, KE Shitang

(College of civil engineering, Nanjing University of aeronautics and astronautics, Nanjing 210016, China)

Using the computational fluid dynamics (CFD) method the numerical wind tunnel was conducted to simulate a 189 m-high cylindrical-conical steel cooling tower, Asia’s highest cooling tower that is still under construction. The features of the airflow around the typical cross-sections and its wake were extracted, and the distribution of mean wind loads along the internal surfaces of the cylindrical and conical sections was obtained. The functions for estimating the internal surfaces shape factors of the cylindrical and conical segments were obtained by fitting to the simulated data. Furthermore, finite element method was used to analyze the static wind-induced response of the cylindrical-conical steel cooling tower under internal pressure, external pressure, or both internal and external pressure. The effect patterns of internal pressures on the wind-induced responses of the main tube, stiffening trusses and auxiliary trusses of the tower were derived from the analysis results. The results show that the inner shape coefficient on the windward surface of the cone section is close to that of the straight section of similar value, but in other parts the deviation is great. The internal pressure effect has greater influence on the negative displacement of the structure, the axial force of the horizontal bar in the negative external pressure area of the main cylinder, the axial force of the auxiliary truss, and the axial force of the upper reinforcement truss.

steel cooling tower; cylinder−conic section; numerical wind tunnel; average wind pressure; wind-induced response; effect of internal pressure

TU279.7

A

1672−7207(2018)03−0684−12

2017−05−22；

2017−08−02

江苏省优秀青年基金资助项目(BK20160083)；国家自然科学基金资助项目(51208254)；中国博士后科学基金资助项目(2013M530255，1202006B) (Project(BK20160083) supported by the Excellent Youth Foundation Program of Jiangsu Province; Project(51208254) supported by the National Natural Science Foundation of China; Projects(2013M530255, 1202006B) supported by the Postdoctoral Science Foundation of China)

柯世堂，博士，副教授，从事结构工程与风洞试验研究；E-mail: keshitang@163.com