卢旦突然开孔时孔口气流动力特性参数的数值模拟重要

卢旦突然开孔时孔口气流动力特性参数的数值
模拟重要
Pleasure Group Office【T985AB-B866SYT-B182C-BS682T-STT18】
第23 卷第10 期工程力学
2006 年10 月Oct. 2006 ENGINEERING MECHANICS 55
文章编号：1000-4750(2006)10-0055-06
突然开孔时孔口气流动力特性参数的数值模拟
卢旦，*楼文娟
(浙江大学建筑工程学院, 杭州310027)
摘要：建筑物突然开孔时瞬态内压响应的规律可以用一个二阶非线性常微分方程来描述，由于开孔形状和孔边物理特征的复杂性，采用风洞试验获取方程中各项系数的方法往往比较困难。

首先利用计算流体动力学的方法，模拟建筑物突然开孔瞬时的流场变化，获得风致内压的Helmholtz 频率和孔口气流的等效线性阻尼。

再由这两个
参数计算得到的内压增益并与风洞试验结果相比较，二者结果吻合较好，表明数值计算能够准确地模拟突然开孔结构孔口处的气体流动状态。

其次，利用孔口气流振荡曲线，采用参数拟合的方法进一步获得了内压传播方程中的各特征参数。

关键词：风洞试验；风致内压；数值模拟；突然开孔；动力特征参数
中图分类号：文献标识码：A
NUMERICAL SIMULATION OF FLOW DYNAMICAL PARAMETERS AT
SUDDEN OPENING
LU Dan , *LOU Wen-juan
(Department of Civil Engineering, Zhejiang University, Hangzhou 310027, China)
Abstract: The transient response of internal pressure following a sudden openning was described by a
second-order, non-liner, ordinary differential equation. Because of the complexity of the geometry and the
opening configuration, it is quite difficult for wind tunnel test to correctly estimate the constant coefficients in the
equation. The transient response of flow pattern at sudden opening is simulated by Computational Fluid Dynamics
method. Both the Helmholtz frequency and equivalent viscous losses are obtained. The gain function obtained by
using these parameters coincides with the result of wind tunnel test fairly accurately. It is shown that the
numerical computation can simulate the flow pattern at sudden opening correctly. Finally the flow dynamical
parameters are derived by fitting the differential equation to the numerical response.
Key words: wind tunnel test; wind-induced internal pressure; numerical simulation; sudden openning;
dynamical parameter
传统上人们对结构抗风的设计往往只考虑建Helmholtz 共振等的动力效应，都将导致建筑物内
筑物的外表面所承受的风荷载。

但在风灾天气下，
部风压增大[1]，使屋盖、幕墙等结构受内、外压共
建筑物的门、窗、幕墙等围护结构特别容易被强风
同作用而更容易遭受风致破坏[2]。

1979 年Holmes[3]
吹开，风从开孔突然涌入。

无论是内压静力效应还第一个提出具有单一房间和单一开孔的建筑物就是由于突然开孔所引起的内压超载和紊流引起的像声学中的Helmholtz 谐振器，室内空气在孔口处———————————————
收稿日期：2004-12-03；修改日期：2005-04-06
基金项目：国家自然科学基金()
作者简介：卢旦(1978)，男，浙江上虞人，博士生，从事结构风工程研究(E-mail)；
*楼文娟(1963)，女，浙江绍兴人，教授，博士，博导，从事结构风工程研究(E-mail)。

56 工程力学
进出振荡的频率称为Helmholtz 频率。

研究内压响方程(1)满足伯努利方程，稳态流中孔口处的气
应的瞬态过程，并由此获得Helmholtz 频率及孔口流将会发生收缩，因此方程中出现孔口收缩系数c。

气流阻尼能够方便地计算出内压响应的增益，为进Vickery 和Bloxham 研究认为在非稳态流中将
一步研究开孔结构的风振响应提供条件。

不出现孔口气体收缩，因此方程中不应该出现c，随着计算机硬件和软件技术的日新月异，采用而代以另外两个系数C I和C L：
计算流体动力学(CFD)对钝体绕流进行数值模拟分析，获取建筑物表面风压及风速流线的技术已经比ρVρqV
L′&&+| &| &+=
2
a C C C C C
e0 a0
C
pi L pi pi pi pe
γA p2(γ 2
A p)
0 a
0 a
(3)
较成熟[4,5]。

CFD 技术不但能节约试验成本、缩短H elmholtz 频率的表达式也不包含流动系数c
试验周期，而且与试验相比可以获得非常详尽的资料，有助于对问题发生的机理进行研究。

Sharma
1 γA p
f( =0 a(4) HH kery) ′
Vic
2πρL V
a e0
和Richards[6]采用数值模拟研究了结构突然开孔后其中L e′=L0 +C I A0 是气流的有效长度，C I是惯性内压发生的静脉收缩式振动现象，比较了紊流引起
系数，C L是损失系数。

的粘滞阻尼损失和内外气流相互作用引起的能量
文献[6]建议，对于墙面较薄的开孔采用下式：
损耗。

本文通过CFD 数值模拟，获得突然开孔结构内压瞬态响应方程中的Helmholtz 频率及孔口气流
的等效线性阻尼，计算出内压响应的增益，使之与试验结果相比较，研究利用流体动力学数值计算方ρρ 2
L V&&qV| &| &
a e C C C C
0 +a0 pi pi
γcA p2(γA p)2
pi L
0 a0 a
?
?
?
?
μ?
eff?
?
PL V
?&
r
e0
?
C
pi
γ2
c A p
2
0 0
+C
pi
=C
pe
+
(5)
法来模拟开孔结构孔口处气流的运动状态，并由此
对于墙面较厚、孔道长度大于孔道直径的开获得相关参数的可行性。

孔，建议采用如下方程：
1 内压传播方程理论ρ′ρ 2
L V&&qV| &| &
a C+C C C
e0 a0
γA p2(γ)2
pi L pi pi
A p
0 a0 a
+
迄今为止，各国学者已经在内压响应方程类型上达成共识，即认为可以用二阶常微分方程来描述突然开孔结构的内压响应。

Liu 和Saathoff[7]采用非?
?
?
?
μ?
eff′
?
PL V
?
?&
r
e0
?
C
pi
γA p
2
0 0
+C
pi
=C
pe
(6)
定常等熵伯努利方程对孔口处入射气流进行分析μ是考虑紊流的有效摩擦系数，?r为流速从平均
eff
后得到：
ρL&&+ρ| &| &+=
V qV
C
2
a C C C C
e0 a0
(1)
pi pi pi pi pe γcA p2(γcA p)2
cA p cA p)
0 a0 a
式中，γ、ρa、p a分别是开孔周围空气的比热、密度和压强；A0 开孔面积，V0 是建筑物内部容积；C,C分别为内压力系数和外压力系数，q是参考pi
pe 速度到达零速度时经过的距离，P为孔口周长。

尽管方程(5)、方程(6)已经明确表达了瞬态内压响应的规律，但方程各项系数的取值却受到孔口物理性质的影响。

不同的问题对应有不同的系数，而问题的关键在于开孔处流动状态的区别。

由于受到试验条件的限制，通过风洞模型试验获取这些参数
点风压，即q=1/ 2ρ
2
，v h为参考点风速；c是气
a
v点风压，即q=1/ 2ρ 2 ，v h为参考
点风速；c是气
h
的方法往往会很困难。

而数值模拟获取参数的可行
流的孔口收缩流动系数；L e=L0 +αA0 / π是孔口性和可靠性正是本文所要研究的主题。

处空气柱的有效长度，L0 是开孔的实际深度，一般
指墙体的厚度，系数α依赖于开孔的形状等其他因
2风洞试验
素需由试验测得。

上式中的一至四项可分别被视作风洞试验刚性模型选用优质木材做成，屋面为惯性项、阻尼项、弹性项和系统所受的外力。

在无800mm×500mm 的长方形，底裙高200mm，屋面
阻尼自由振动条件下(即令方程的第二项和第四项板和底裙之间为固支，用胶带纸密封。

测压管为外为零)求解(1)式便得到Helmholtz 频率：径1mm 的铜管，埋在屋面板上下表面中间预先设
1 γcA p
f( =(2)
0 a
HH/ S) ρ
Liu
2πL V
a e0 计的夹层内。

为了使从测压点处引出的铜管不影响模型内部风场，在屋面板下的模型内设置空心柱
工程力学57
子，屋面的测压管通过空心柱与底部的测压模块相由于开孔前内、外压存在压力差，因此开孔瞬连接。

模型在迎风面墙板正中开孔，开孔面积为间外压的激励类似于一个矩形脉冲。

而从图2 所示A =×。

试验时先将孔口用一挡板密封，
的内压系数时程中可以看出，开孔瞬间内压突然增
在风洞中吹几秒钟后用预先连在挡板上的绳子突大，但由于试验中使用的是湍流流场，实际内压响然拉开，模拟实际建筑门窗产生风致破坏的情形。

应含有脉动风压，从其时间历程中不易看出模型如图1 所示。

Helmholtz 振动；又由于本模型Helmholtz 振动的能
量较小，因此依据其功率谱也难以获得Helmholtz
振动的信息。

3 CFD 数值模拟
采用计算机进行数值模拟不但能对结构进行
准确地三维计算，采用非稳态计算还能实时记录下
流体的压强、流速矢量等参数的变化，逼真地再现
气流的运动变化规律。

而这些都是风洞试验所难以
或无法做到的。

数值模拟对象为前文中图1 所示的风洞模型。

图 1 突然开孔风洞试验刚性模型
计算流域为长×宽×高＝6m×3m×2m，符合风洞Stiff model of the structure with a sudden opening
试验阻塞比小于4%的要求，采用非结构化网格划由于来流在屋檐处分离造成迎风屋面边缘附
近形成高负压区，这种高负压区仅限于边缘区域[8]，
分，孔口和建筑物附近网格加密。

为提高计算效率，
考虑到结构的对称性采用对称建模，在计算过程中所以在屋面上表面的边缘测压点布置较密。

考虑到
采用自适应网格技术。

计算网格模型如图3 所示。

室内气压较均匀，所以在屋面下表面及内墙面的测
点进行均匀布置。

由于模型具有对称性，因此对模
型屋盖上、下表面的测点都采取1/2 单边布置(见图
1)，但在迎风墙面孔口附近重点研究位置则采取整
个墙面对称布置。

本次风洞试验在南京航空航天大学603 研究所
NH-2 低速风洞中进行。

试验段截面为矩形带小切
角，高，宽，长。

在试验段的进风
口设置格栅以形成一定的湍流度。

数据的采集由多
通道电子扫描脉动风压测量系统完成。

图 3 开孔结构模型的计算网格2
Computing grid of the opening model
因为是迎风面开孔，同时又考虑到孔道内气体
流动具有强烈的各向异性，因此湍流模型选用可实
1 现的k ?ε模型(Realizable k ?ε)。

流场为均匀湍流场，参考风速14m/s，湍流度I =20%。

由于脉动风
对平均风的影响是通过湍流度计入运算过程，因此计算结果得到将是开孔瞬时平均风速的变化过程。

0 5 10 15 计算分为两步：首先，关闭(deactivate)结构内部网格进行稳态计算(steady-state solution)，用来模拟开
时间t/s 孔前的流场环境，计算收敛的结果保存作为下一步图 2 内压系数时程非稳态(unsteady-state)计算的初始条件。

第二步，激Time history of internal pressure coefficient 活(activate)结构内部网格使之参与非稳态计算。

根
58 工程力学
据初步估算，结构的Helmholtz 频率约为60Hz，由内孔口附近气流仍是往内回的。

由此可见，出风时
此设定计算时间步长为。

值得指出的是，通不但孔口有效面积较小，而且由于气流的复杂流
过改变时间步长可以控制突然开孔的延续时间，这动，阻尼损耗也较大，整个气流循环过程是一个严一点是风洞试验无法精确做到的，同时计算时间步重的非线性过程。

长的大小只会影响响应峰值的大小(步长越小峰值由于Helmholtz 频率和孔口气流阻尼是开孔结
越高)，但不会改变Helmholtz 频率和阻尼值。

本次构内压响应所固有的特性，它们只和孔口本身的物
计算时间步总数为100 步，为节省计算机资源采取理性质有关，因此这些参数同样可以从进出孔口气自适应迭代，每个时间步最大迭代次数为30 次。

流的瞬态响应中获取。

图5 是数值模拟得到的突然与以往研究风荷载时考虑空气为不可压缩气开孔后孔口气流的振荡曲线。

体不同的是，此次问题必须把空气视为可压缩气体
来研究。

即空气的密度和压强必须满足以下公式：
30
原始曲线
趋势曲线
衰减曲线p p
20
=
a
=const
(7)
ρρ
a
10
为提高非稳态计算的精度，对可压缩气体方程的求
解采用2 阶couple implicit 格式算法。

-10
-20
时间t/s 室外孔口
室内
图 5 数值模拟进出孔口气流的振荡曲线
Computational oscillating curve of flow following
sudden opening
结构开孔前孔口风速的平均值为零，在开孔之
(a) 进风后的时间，由于局部风场环境发生了变化，
风速也发生了变动。

图5 可以看到，风速中明显包
含了Helmholtz 振动。

且气流在开孔后时间内
是处于进风状态，在之后的时间内全都处于出风状
态，直至趋向于零风速。

为计算Helmholtz 频率和
孔口阻尼，还需对图中原始曲线进行处理，首先采室外孔口室内
用HHT(Hilbert／Huang transform)方法[9]获得原始
振荡曲线的趋势线(图5 中虚线表示)。

该方法从本
质上讲是对一个信号进行平稳化处理，其结果是将
信号中不同尺度的波动或趋势逐级分解开来，产生
(b) 出风
一系列具有不同特征尺度的数据序列，每一个序列图 4 孔口附近风速矢量图
称为一个本征模函数(Intrinsic Mode Function，IMF) Velocity vectors at the opening
分量。

对经EMD 分解得到的各IMF 分量再进行线图4(a)为开孔后时刻开孔位置处气流从
性、平稳处理即Hilbert 变换，使得变换后得到的结模型外部流入模型内部的矢量图，图4(b)为开孔后
果能够反应真实的物理过程。

将原始数据减去所有时刻开孔位置处气流从模型内部流出模型内
IMF 分量的总和即为原始型号的趋势曲线。

部的矢量图。

从图4 可见，进风和出风均有孔口收
原始曲线减去由HHT 方法得到的趋势曲线后，缩现象，并且进风时在室内的孔口附近、出风时在
室外的孔口附近均伴有涡流产生。

此外，出风时室
可以得到趋势为零的自由衰减曲线(经EMD 分解和
H ilbert 变换后得到的所有IMF 分量总和)。

至此，
工程力学59
便可由常用的计算单自由度振动系统频率和阻尼(10)得出内压增益。

图6 为分别为由CFD 计算和风
的方法计算Helmholtz 频率和孔口阻尼。

其中，孔洞试验得到的内压增益比较。

试验和数值模拟获得
口气流振荡的等效线性阻尼及其他参数阻尼的确定可以分别通过等效线性方程和非线性常微分方的Helmholtz 和等效线性阻尼系数基本一致，分别为f HH=43 Hz，c eq=35 H z。

程参数拟合得到。

从图6 中可以看出，CFD 计算和风洞试验到的4 内压增益及参数确定内压增益吻合得非常好。

说明利用CFD 数值模拟
能够十分准确地获得开孔结构的Helmholtz 振动和增益函数| H(ω) |是反映振动系统动力特性的
孔口气流的阻尼。

最重要的函数之一，定义为响应与激励的傅立叶变为了进一步准确获得孔口气流的非线性动力换之比的模。

它在振动理论中占有很重要的地位。

特性，计算孔口气流运动的动力特征参数，也即方本文讨论的内压增益是指内压时程(内压系数时程) 程(5)或方程(6)中的常系数值。

可以利用图5 中的衰
与孔口位置处外压时程(外压系数时程)的傅立叶变减曲线，采用多变量寻优拟合的方法确定这些常系换之比的模。

选择内压增益作为数值模拟和风洞试数。

根据前文介绍，孔口气流的动力特性只和孔口验的比较指标，可以有效分析数值模拟的准确性。

本身的物理性质有关，因此气流的动力特性参数同对于风洞试验只需分别测出结构外表面开孔样也可以从进出孔口气流的瞬态响应中获取。

首部位附近的风压系数和结构内部风压系数的时程先，利用图5 中的衰减曲线分别计算出进出口风速数据便可由下式获得内压增益。

|
fft(C)
Hω=(8)
( ) |
pi
fft(C)
pe
的一阶和二阶导数，再对非线性方程(5)或方程(6)
的常系数进行参数拟合。

由于系数中ρa、V0 、
γ、
p、A0 、q和P这几个参数都可以由模型的几何
a
8
试验数据尺寸和风场环境方便地获得，因此方程中剩余的L e/c、C L、μeff/(?r?c) 等系数便可随之获得。

本
CFD 模拟
文模型采用方程(5)，拟合所得常系数和对应孔口动
6 4 力特性参数为：
ρL V L e
a e
0 =1 ×?5 ?=
.3 10
γcA p c
0 a
(11)
2
ρ 2
qVρ 2
C 2 9 10 1 25
L C
a0 . ?=.
=×?7
L
2( A p)γ 2
0 a
(12)
0 20 40 60 80 100
频率f/Hz
图 6 CFD 计算与风洞试验的增益函数比较?
?
?
?
μ?
eff
?
PL V
?
??5
r
e0
?
=×10 ?
(13)
γ 2 2
c A p
0 a
μ
eff
=0 ?-2 ?s-
1
.12kg m
?r?c
Comparison of gain functions obtained from
measurements and CFD
对于CFD 计算，首先将内压传播方程线性化[10] C&&+&+ω 2 =ω 2 (9)
pi c C C C
eq pi HH pi HH pe
式中，c eq称之为等效线性阻尼系数。

再将方程(9) 经拉普拉斯变换得到内压增益的计算公式。

根据孔边材料的物理性质或相关材料试验，例如取μeff等于层流粘性系数×10-5kg·m-1·s-1)，
?r mm 。

则由式(13) 可计算得到本模型的=
c，再由式(11)可得L e==L0 +。

=
5 结论
?1/ 2
??
2
??
ωH2 H c
ω2
??ω?
??????
eq
|
H( ) |=+
ω(10)
????
?ωω?
2 2
????
??
HH HH
因此，只要CFD 模型能准确计算出Helmholtz
频率ωHH=2πf HH和等效线性阻尼系数，便能由式
突然开孔结构的瞬态内压传递函数已被公认为是一个二阶非线性微分方程，方程各项系数随孔口物理性质的不同而不同。

若能利用数值模拟获得上述系数将大大节省研究成本。

本文利用CFD 数值计算模拟了突然开孔结构的孔口气流变化，并与
60 工程力学
on high buildings [J]. Journal of Xi’an Jiaotong
University, 2001, 35(5): 471~474. (in Chinese)
试验结果相比较得出以下结论：
(1) CFD 数值模拟能对突然开孔结构的孔口气
[5] 卢旦, 楼文娟, 陈勇, 唐锦春. 双幕墙建筑通风性能的流变化进行准确地模拟。

数值模拟研究[J]. 浙江大学学报, 2005, 39(1): 46~50.
Lu Dan,
Lou
Wenjuan,
Chen
Yong,
Tang
Jinchun.
Study
(2) 采用内压增益作为比较指标，可以直观地
on
ventilation
in
double-
skin
facades
by
numerical
检验数值模拟结果的正确性。

simulation
[J].
Journal
of
Zhejiang
University,
2005,
(3) 模拟得到的响应曲线可以对瞬态内压传递
39(1):
46~50. (in
Chinese) 方程中的各项系数进行寻优拟合，从而获得突然开
孔时孔口气流动力特性参数。

[6] Sharma R N, Richards P J. Computational modeling of
the transient response of building internal pressure to a
sudden opening [J]. Journal of Wind Engineering and
参考文献：
Industrial Aerodynamics, 1997, 72: 149~161.
[7] Liu H, Saathoff P J. Building internal pressure: sudden
[1] Liu H, Saathoff P J. Internal pressure and building safety
[J]. J. Struct. Div. ASCE, 1982, 57(3): 2223~2234.
[2] 布占宇. 大跨度柔性屋面结构突然开孔时的风振响应
研究[D]. 杭州: 浙江大学, 2003.
Bu Zhanyu. Study on wind-induced dynamic response of
long-span flexible roof structures in the case of sudden
opening [D]. Hangzhou: Zhejiang University, 2003. (in
C hinese)
[3] Holmes J D. Mean and fluctuating internal pressures
induced by wind [C]. Proc. 5th Int. Conf. on Wind
Engineering. Fort Collins, USA, 1979: 435~450.
[4] 周莉, 席光. 高层建筑群风场的数值分析[J]. 西安交
通大学学报, 2001, 35(5): 471~474.
Zhou Li, Xi Guang. Numerical analysis of the wind field
chang [J]. J. Eng. Mech. Div. ASCE, 1981, 107(Em2):
309~321.
[8] Uematsu Y,Yamada M, Sasaki A. Wind-induced dynamic
response and resultant load estimation for a flat long-span
roof [J]. Journal of Wind Engineering and Industrial
Aerodynamics, 1996, 65: 155~166.
[9] Huang N E, Shen Z, Long S R. A new view of nonlinear
water waves: the Hilbert spectrum [J]. Ann Rev Fluid
Mech, 1999, 31: 3417~3457.
[10] Sharma R N, Richards P J. Computational modeling in
the prediction of building internal pressure gain function
[J]. Journal of Wind Engineering and Industrial
Aerodynamics, 1997, 67 & 68: 815~825.
(上接第67 页)
[11] Wang Xu, Shen Yapeng. Inclusions of arbitrary shape in
magnetoelectroelastic composite materials [J]. I. J.
Engineering Science, 2003, 41(1): 85~102.
[12] Liu Jinxi, Liu Xianglin, Zhao Yongbin. Green's functions
for anisotropic magnetoelectroelastic solids with an
elliptical cavity or a crack [J]. I. J. Engineering Science,
2001, 39(12): 1405~1418.
[13] Tian Wenye, Gabbert U. Multiple crack interaction
problem in magnetoelectroelastic solids [J]. European
Journal of Mechanics A/Solids, 2004, 23: 599~614.
[14] 孙建亮, 周振功, 王彪. 功能梯度压电压磁材料中断
裂问题分析[J]. 力学学报, 2005, 37(1): 9~14.
Sun Jianliang, Zhou Zhengong, Wang Biao. A permeable
c rack in functionally grade
d piezoelectric/piezomagnetic
m aterials [J]. ACTA Mechanica Sinica, 2005, 37(1):
9~14. (in Chinese)
[15] George R. Buchanan. Layered versus multiphase
magneto-electro-elastic composites [J]. Composites: Part
B, 2004, 35: 413~420.
[16] Garcia Lage R, Mota Soares C M. Layerwise partial
mixed finite element analysis of magneto-electro-elastic
plates [J]. Computers and Structures, 2004, 82:
1293~1301.
[17] Ding Haojiang, Jiang Aimin. A boundary integral
formulation and solution for 2D problem in
magneto-electro-elastic media [J]. Computers and
Structures, 2004, 82: 1599~1607. [18] 丁皓江, 江爱民. 横观各向同性磁电弹性体的基本解
和边界积分方程[J]. 中国科学(E 辑), 2003, 33(9):
845~855.
Ding Haojiang, Jiang Aimin. The fundamental solutions
for transversely isotropic magnetoelectroelastic media
and boundary integral formulation [J]. Science in China
(Series E), 2003, 33(9): 845~855. (in Chinese)
[19] 孙焕纯, 张立洲. 无奇异边界元法[M]. 大连: 大连理
工大学出版社, 1999.
Sun Huanchun, Zhang Lizhou. Nonsingularity boundary
element methods [M]. Dalian: Dalian University of
Technology Press; 1999. (in Chinese)
[20] Sun Huanchun, Yao Weian, Virtual boundary
element-linear complementary equations for solving the
elastic obstacle problems of thin plate [J]. Finite
Elements in Analysis and Design, 1997, 27 (2): 153~161.
[21]姚伟岸, 王辉. 压电材料平面问题的虚边界元——等额
配点法[J]. 计算力学学报, 2005, 22(1): 42~46.
Yao Weian, Wang Hui. Virtual boundary element-
equivalent collocation method for plane piezoelectric
materials [J]. Chinese Journal of Computational
Mechanics, 2005, 22(1): 42~46. (in Chinese)
[22] Yao Weian, Wang Hui. Virtual boundary element
integral method for 2-D piezoelectric media [J]. Finite
Elements in Analysis and Design, 2005, 41: 875~891.。