非线性波浪变形计算的三维边界元方法

第4期
孙大鹏 、李玉成 : 非线性波浪变形计算的三维边界元方法
399
单纯地加大阻尼消波段长度 B 以谋求消波效果的做法无疑是不经济的 。为此 , 在 Sponge layer 阻尼消波段的后边界 S 2 上 , 设置透过边界 , 即波势函数满足 Sommerfeld 放射条件 。即 Sponge
摘要 : 采用 0 - 1 混合型边界元剖分计算域边界面 , 以提高边界元方法模拟三维波场波浪变形问 题的数值计算精度 。采用 Sponge layer 阻尼消波和 Sommerfeld 放射条件相匹配的边界处理方式 , 以 消除计算域内的反射波 。借助时域内的波面位置追踪 , 实现了三阶 Stokes 波在数值波动水槽中的 波浪变形计算 , 其结果吻合于理论值 。 关 键 词 : Laplace 方程 ; 非线性波浪 ; 三维边界元 ; 0 - 1 混合元 中图分类号 : TV 13912 文献标识码 : A 文章编号 :100126791 (2002) 042397206
s
x
1
5μ
( on S 2 )
( 9)
式中 C 为波速 ; x1 为 Sponge layer 阻尼消波层的前端位置坐标 ; x 2 为 Sponge layer 阻尼消波层 的后端位置坐标 ; <s 为阻尼消波段内的水面波势函数 。
Байду номын сангаас
2 基本方程和边界条件
2 11 0 - 1 混合元
对于有势问题的边界元数值方法 , 考察非光滑边界面上具有相同势函数值 ( 或势函数法向 导数值) 作为已知值的两个相邻单元 , 经过数值求解后两单元势函数法向导数值 ( 或势函数值 ) 会略有差别 , 究其原因是因为计算域内存在着如图 3 所示的法向导数不连续的角点[ 3 ] 。由于线 性波浪理论下的波浪变形计算是定常场问题 , 计算仅进行一次便告结束 。因此采用常数元剖分 边界所表现出来的横向 ( 波峰线上) 差别对计算结果精度的影响不很明显 。但是在做时域内的非 线性波浪计算时 , 上述误差的积累 , 会造成形如图 4 所示的 “横向振动” 。鉴于势函数在整个 边界面上的连续性和势函数法向导数在个别角点上的非连续性 , 这里构造一种单元 , 即对单元 上的势函数法向导数采用常数元 , 称之 “0 ”型元 ; 而对于单元上的势函数采用线性元 , 称之 “1”型元 。而用这样的 0 - 1 混合元剖分图 1 的三维波场边界面 , 时域内的波浪数值计算得以 顺利实施 。
若在四节点的四边形单元内基函数取作 : Ψ1 = 1 ( 1 - ξ ) (1 - η ) ; Ψ2 = 1 ( 1 + ξ ) (1 - η ) 4 4 Ψ3 = 1 ( 1 + ξ ) (1 + η ) ; Ψ4 = 1 ( 1 - ξ ) (1 + η ) 4 4
( 10)
为了合理地解释和定量地描述起因于有限振幅波动的波浪传播非线性现象 , 数值模式应该 考虑自由表面水质点的运动速度以及波面斜率 。对于计算域内的有势波动来讲 , 无论是微幅波 还是有限振幅波 , 其速度势函数都应满足基于水体运动质量守恒原理而导出的 Laplace 方程 , 而借助 Green 公式转换后 , 两者都应具有同样形式的积分表达式[1 ] , 所不同的是有限振幅波在 自由水面上 , 应满足以伯努利方程给出的关于波势函数和波面高度的非线性边界条件 , 从而决 定了非线性波浪变形数值计算是一个时域内具有初始值的边值问题 。在计算技术方面 : ( 1) 在 时间上 , 每个时步都要首先确定自由表面的波动位置 , 进而剖分三维波场边界面 。为了模拟非 线性波的传播特征并得到稳定的计算结果 , 数值计算的计算时间长度不宜太短而时间步长又不 应太大 ; ( 2) 在空间上 , 要设定足够长的传播距离作为波浪变形的计算区域 , 同样出于计算稳 定方面的考虑 , 较之线性模式非线性数值计算时边界应该剖分得更细 ; 而在数值模式的建立方 面 , 时域内的波面追踪方法 、水面动边界的剖分方式以及计算域内边界条件的引入和处理是三 维非线性波浪数学模型所要解决的关键问题 。 针对非线性波浪变形数值模式的上述特征 , 本文采用 0 - 1 混合元剖分波场边界面和
第 13 卷 第 4 期 2002 年 7 月
水
科
学
进
展
ADVANCES IN WATER SCIENCE
Vol113 ,No14 Jul. ,2002
非线性波浪变形计算的三维边界元方法
孙大鹏 , 李玉成
( 大连理工大学海岸和近海工程国家重点实验室 , 辽宁 大连 116024)
t = t0
= η( x , y , t 0 )
( 12)
则此时自由表面上的未知量仅为势函数法向导数5 </ 5 n 。在进行下一时步 t 0 +Δt 的数值 计算之前 , 作为计算的初始值 , 利用台劳级数给出 t 0 +Δt 时刻波面位置和波面势函数 :
X ( t 0 + Δt ) = X ( t ) + Δ X Y ( t0 +Δt ) = Y ( t ) +Δ Y Z ( t 0 + Δt ) = Z ( t ) + Δ Z ( 13)
T ij
<j ( i = 1 , 2 , …, N ; j = 1 , 2 , …, m )
( 11)
式中 N 为边界面上的节点总数 ; m 为单元数 ; <j 为 j 单元的势函数法向导数 。
hij = qij1 = qij3 =
T
κr
1
j
D ( x , y , z) dξ dη; ) D (ξ,η
qij = { qij1 , qij2 , qij3 , qij4 } ;
T
κ κ
Ψj1 5 1 5n r j Ψj3 5 1 5n r j
=
q = D ( x , y , z) dξ dη; ij2 ) D (ξ,η D ( x , y , z) dξ dη; ) D (ξ,η qij4 =
5
1
r
ds
( 4)
设定三维波场如图 1 , 计算域边界面由波浪入 射边界 S 0 、自由表面 S F 及固定边界 S B 构成 。 在入射边界上 , 给出运动学边界条件 , 即波动 水质点的法向流速 : 5 <( x , y , z , t ) = q ( t ) ( x , y , z ) ∈ S 0 ( 5) 5n 在海底及固体边界上 , 采用全反射边界条件 , 即
<( x , y , z , t0 +Δt ) = <( x , y , z , t0) +Δ < ΔX 、 ΔY、 Δ η及Δ < 为Δt 时段上波面水质点位置坐标和势函数的增量 。忽略式 ( 13) 中 式中 2 Δt 以上的高阶项 , 即 η Δ X = d XΔt ; ΔY = d Y Δt ; Δ η = d Δt ; Δ< = d < Δt dt dt dt dt 而在自由表面上 :
图2 Sponge layer 阻尼消波
Fig12 Arrange of the Sponge layer for wave damping
图1 三维波场的边界面构成
采用 Sponge layer 阻尼消波方式 , 即在反射边界前 设置一距离为 B 的阻尼消波段 , 如图 2 , 以消除计
算域内的反射波 。而以往对 Sponge layer 阻尼消波的应用经验表明[ 2 ] , 对于波长大于消波段长 度 B 的波浪 , Sponge layer 的消波效率并不高 。而且 , 从计算机容量及占用机时的两方面考虑 ,
5 <( x , y , z , t ) Fig11 Definition sketch of wave field = 0 ( x , y , z ) ∈ S B ( 6) 5n 因自由表面是非固定边界 , 须同时给出波面水质点的运动边界条件和动力边界条件 , 并认 为水面上压力为常值 , 即 dη 5 < = ( 7) dt 5 z z =η( x , y , t) 5< 1 5< 2 5< 2 5< 2 + + + 5t 2 5x 5y 5z ( 8) + gη = 0 | z =η( x , y , t) η( x , y , t ) 为相对于静水面的波动位置高度 。 式中
利用上式的基函数将任意的四边形单元转换为 (ξ,η) 坐标下的正方形基本单元 , 并用 j1 、
400
水
科
学
进
展
第 13 卷
j2 、j3 、j4 表示 j 单元的四个节点 , 采用 0 - 1 混合元剖分边界 , 则式 ( 2) 的积分可离散成 ci <i +
6
j
m
qij <j =
T
6h
j
m
T 即在每个单元上计算 qij 形成单元矩阵系数 , 合成后得到式 ( 11) 线性方程组的总体系数矩阵 。
2 13 时域内的波面追踪
在具有初始值的自由表面 S F 上 , 若已知 t 0 时步波面位置高度和波面势函数 : η( x , y , t ) | t = t0 = η( x , y , t 0 ) <( x , y , z , t ) |
ϕ ∫
t
( 2)
上式中的 CM 是仅与边界形状和边界剖分方式有关的系数 , 在完成边界剖分后 , 可以首先 对于一个等势场问题 : < ( x , y , z ) = const ( 3) 5< = 0 5n 求解式 ( 2) 的积分方程 , 得到 :
cM = 1 12 边界条件
κ5n
Γ
t
398
水
科
学
进
展