非线性晃动问题的 ALE 边界元方法

标, 网格不动, 此时变为 Eu ler 描述。
第1期
岳宝增等: 非线性晃动问题的AL E 边界元方法
3
3 AL E 网格速度的确定与网格更新 在 AL E 描述下, 网格可以按任意合适的方式进行设计, 但为了能够达到跟踪边界尤其
是流体晃动时的自由面, 其上的网格点必须永远在边界上, 这样网格点的速度设计就受到限 制即速度的约束条件。 假设边界面 (或线) 的方程为:
由 (5) 式, 自由液面上网格结点移动的速度为:
uδ1 = 0
(6)
uδ2 =
n n
1 2
u
1
+
u2
(7)
在求解过程中, 每一时间步网格需要更新一次, 对于 tn+ 1时刻, 网格结点的坐标可按下式计 算:
∫ x
n+ i
1
=
x
n i
+
tn+ 1
w id t
(8)
tn
4 数学模型
本文以流体的不可压缩, 无粘, 非定常无旋理论为基础, 物理模型如图 2
2
宇航学报
第 19 卷
面上的动力边界条件 (B ernou lli’s) 及自由液面的位置随着时间的变化而移动, 因此, 求解晃 动问题是一个既困难又麻烦的课题。有限元法被应用于求解晃动问题获得了极大的成功。用 有限元法求解液体晃动问题时, 运动的流体与有限元网格之间的关系即动力学描述是极为 重要的。由于拉格朗日方法和欧拉方法这两种经典的描述方法各有长短, 这就使人们把这两 种方法结合起来使用, 得到了AL E 有限元方法[1, 2]。边界元方法以有限元方法为基础, 突出 优点是可将求解空间降低一维, 80 年代以来用以求解带有自由液面的晃动问题[3]取得了极 大进展。 本文推导了 AL E 边界元计算格式及非线性边界元方程的预报——多次校正迭代 格式, 对时间变量采用N ewm a rk 方法进行离散。 在跟踪自由液面时, 提出了一种简便合理 的网格更新方法, 通过设计满足一定条件的网络移动速度使边界上的网格点保持在边界上, 并使单元形状不发生畸变。
55 5n
T
=
5< 5n
,
1
5< 5n
, …,
2
5< 5n
n1
, n1 为自由液面 S 1 上的结点个数, A
为 n 阶方阵, B
为 n×n1 阶矩阵。 对于动力边界条件 (10) , 为了消除解的不稳定性[3, 5], 引入误差修正项 D ,
用下式代替 (10) 式:
5D 5t
=
d d
<
t
-
1 2
55 5n
,
其中插值函数
N
T
=
1 2
(1-
2s l
)
,
1 2
(1+
2s ) l
, 5 T = [ <j , <j+ 1 ],
55 5n
T
=
5< 5n
,
j
5< 5n
j+ 1
, l 为单元长度。 根据边
界元方法[4 ] (9) 式可离散为:
A5 +
B
55 5n
=
0
(12)
而: 5 T = [ <1, <2, …, <n ], n 为 S = S 1+ S 2 上结点个数
流体粘度阻尼效应的所谓 R ayleigh 阻尼项。 Λ 一 般采用实验数据。Αx (t) 为一水平加速度激励。S 1 为
自由液面,
S
2
为贮箱刚性壁。
由:
D D
<t =
55<t +
5< 5x i
5x 5t
i
,
可推出:
D D
<t =
55<t +
5< 5x
2
+
5< 5y
2
, 代入动力学边界
条 件 可 得:
简计为df
dt
v i=
5x i 5t
: 定义为质点速度。
X
w
i=
5x i 5t
xδ: 定义为网格速度。
5xδi 5t
: 定义为质点在参考坐标系中的速度。
X
可以看到:
vi =
5x i (xδ, 5t
t)
=
5x i 5t
xδ +
5x i 5xδi 5xδi 5t
(2)
由 (1) 式得:
f=
Df Dt
=
5f 5t
F (x , t) = 0
(4)
由于流体质点和网格点始终在边界上即永远满足方程 (4) 式, 所以有:
从而有关系: (v- w ) · 则可得
DF Dt
=
0;
dF dt
=
0
F (x , t) = 0, 而边界面方程函数的梯度方向即边界面的法线方向 nο,
v n= w n
(5)
在满足上式的条件下, 自由液面上的网格点的移动速度可按如下及种方式设计。
主题词 液体晃动 边界元法 AL E 边界元方法
THE AL E BO UNDARY EL EM ENT M ETHOD S IN SOL V ING NO NL INEAR SHO SH ING PRO BL EM S
Yue B aozeng W ang Zhao lin Kuang J in lu
1 3
5<j 5n
+
1 5<j+ 1 6 5n
T=
-
1 6
5<j 5n
-
1 5<j+ 1 3 5n
1 6
5<j 5n
+
1 5<j+ 1 3 5n
3
5<j 5n
+
5<j+ 1 5n
H=
1 12
5<j 5n
+
5<j+ 1 5n
5<j 5n
+
5<j+ 1 5n
5<j 5n
+
3
5<j + 5n
1
单元集成后可写成如下形式:
xδ +
5f 5x
ci
i
=
df dt
+
c
f= 0
(3)
其中: ci= u i- w i, 定义为 AL E 描述下的对流速度。 (1) 当 xδ= X , w = v 时参考坐标取为物质
坐标, 网格跟随物质质点, 变为L ag range 描述。 (2) 当 xδ= x , w = 0 时, 参考坐标取为空间坐
c5α+
c15 +
G
55 5n
=
R
(15)
时间变量采用N ewm a rk 方法离散, 化简后可得如下形式的方程:
(c +
c1Χ∃ t) 5 n+ 1 +
Χ∃
tG
1 2
5< 5s
2
-
tgΒ
5< 5n
5< 5s
+
Λ< +
Αx x +
gy =
0
(11)
其中:
d d
<t =
5< 5t
,即x
x
保持不变, 竖直方向跟随自由液面的时间导数, Β 表示自由液面和水
平面的夹角。
5 系统方程的数值离散
对方程 (9) 进行边界元离散, 在每一单元上设:
<=
N T5,
55n<=
NT
f=
D f (x , t) Dt
=
5f 5t
x+
5f 5x i 5x i 5t
=5xδj 5xδi 5t
( i =
X
1, 2, …N Α)
(1)
其中 表示固定 x 时求导, N Α 是时间维数。 x
5f 5t
: 定义为当地导数。
x
5f 5t
xδ:
定义为网格导数,
d<
dt
+
5< 5x
cx
+
5< 5y cy -
1 2
5< 5x
2
+
5< 2 5Y
+ Λ<+ Αx x + gy = 0, 由关系式:
5< 5x
=
55n<sinΒ +
55<s co sΒ
5< 5y
=
55n<co sΒ +
55<s sinΒ
经过化简可得:
图 2 物理模型
d< dt -
1 2
5< 2 5n +
2 AL E 描述方法 物质点 t 时刻在空间坐标中的位置为: x = x (X , t) , 物质点 t 时刻在参考坐标系中的位
置为: x = x (xδ, t) , 一个函数 f 一般表示为空间和时间的函数, 同时在参考坐标及随动坐标 中也可表示为各自坐标和时间的函数, 从而有关系式:
f (x , t) = f (X , t) = f (xδ, t) x , X , xδ 分别为空间坐标, 物质坐标, 网络参考坐标。 而
(D ep artm en t of Eng ineering M echan ics, T singhua U n iverity, Beijing·100084·Ch ina)