Page 37 - 201903

第 38 卷 第 3 期 黄武琼等： 共形时域有限差分技术在戏场八字墙中的应用 319

其中，∆t 表示离散时间间隔；∆h 表示离散空间间 n+ 1 2 ( 1 )
+ u y i, j + × l x (i, j + 1)
隔；n表示离散时间刻度；i和j 分别表示在直角坐标 2
n+ 1 ( 1 ) ]
系x和y 方向的位置。对于曲线边界，需采用共形网 − u y 2 i, j − 2 × l x (i, j) , (10)
格技术。在二维直角坐标系下，曲边经过的共形网 其中，l x 和 l y 分别为质点速度对应的棱边在墙体外
格如图 1 所示，深灰色的区域 1 为墙体，深灰色以外 的长度，A为元胞在墙体外的面积。当l x = l y = ∆h，
的区域 2 为墙体外侧，图中斜线阴影网格为曲边附 A = ∆h 时，式(10) 就退化为常规直角网格的递推
2
近的共形网格，需利用共形网格技术来计算。设声 公式 (6)，可见，共形网格的更新方程可以用于所有
压位于直角网格的中心，不管该中心在区域 2 或区 网格的更新方程。为了便于编程，引入相对长度与
域1。共形网格中，质点速度递推方程与常规网格的 相对面积
递推方程一样，只需处理声压的更新方程 [13] 。
l = l x /∆h, (11)
′
x
y ′
l = l y /∆h, (12)
y
Dh
l x ↼i֒j⇁↽ ′ 2 (13)
Dh A = A/∆h .
ӝ۫
l y ↼i֒j↽ P↼i֒j↽ 共形网格的更新方程最终化为
ӝ۫ l y ↼i⇁֒j↽ p n+1 (i, j)
l x ↼i֒j↽ n ∆tρ 0 c 2 [ n+ 1 2 ( 1 ) ′
y
O x = p (i, j) − u x i+ , j l (i + 1, j)
∆h · A (i, j) 2
′
图 1 直角坐标系下的共形网格 ( )
n+ 1 2 1 ′
y
Fig. 1 Conformal grids in rectangular coordinates − u x i − , j × l (i, j)
2
(
n+ 1 2 1 ) ′
对式(3)两边进行二重积分，得 + u y i, j + 2 l (i, j + 1)
x
∫∫ ( ) ]
∂p n+ 1 1
′
dxdy − u y 2 i, j − × l (i, j) . (14)
x
∂t 2
D
∫∫ ( ) 共形网格技术的一个关键点在于相对长度和
= − ρ 0 c 2 ∂u x + ∂u y dxdy. (7)
∂x ∂y 相对面积的求取。若曲线边界是由规则的方程决
D
由格林公式，式(7)右边可化为 定，可由方程与网格的交点算出相对长度，以下部
∫∫ ( ) 分讨论的椭圆房间与八字墙都是规则方程，可以容
− ρ 0 c 2 ∂u x + ∂u y dxdy
∂x ∂y 易计算得到相对长度。共形网格的相对面积一般
D
I 直接取规则网格面积或者规则网格面积的一半，显
= − ρ 0 c 2 u x dy − u y dx. (8)
L 然这样处理太粗糙，本文利用双加权平均法求得相
综合式(7)与式(8)，可以得到 对面积 [14] 。共形网格处理受到变形网格特性的制
∫∫
∂p 约，共形面积过小会使程序不稳定，因此，本文利用
dxdy
∂t SCFDTD 方法 [8] ，将相对面积小于 1/6 的网格近似
D
I
= − ρ 0 c 2 u x dy − u y dx. (9) 为常规网格面积的 1/6，以此来改善计算的稳定性
L
问题。
设墙体边界为刚性边界，那么墙体上的质点速
度为零，仅需考虑网格中墙体外的质点速度的贡献， 2 精度验证
于是由式(9)可得共形网格的FDTD递推方程：
为了验证共形网格技术的精确度，模拟一个
p n+1 (i, j)
椭圆房间内的声场变化。根据费马原理，波总是沿
∆tρ 0 c 2 [ n+ 1 ( 1 )
n 2 着 “波程” 为极小值的路径传播，从椭圆一焦点发
= p (i, j) − u x i + , j l y (i + 1, j)
A(i, j) 2
射出的声线，经过界面的各种反射，最终会聚集在
n+ 1 ( 1 )
− u x 2 i − , j × l y (i, j) 另一焦点 [15] 。以下利用共形网格技术和传统的阶
2