时间步计算

configure.dat输入参数：

  • fluid（水）
    • Density (kg/m^3) : 1000
    • Bulk modulus (Pa) : 220000
  • float（浮筒）
    • Density (kg/m^3) : 6000
    • Elastic modulus (Pa): 1.5e12
    • Thickness (m): 0.01
  • gas（浮筒内空气）
    • Density (kg/m^3) : 12.6
    • Initial pressure (Pa): 120000
  • strip（绑带）
    • Density (kg/m^3): 3500
    • Elastic modulus (Pa): 1.5e9
    • Thickness (m): 0.01
  • dt（CFL stable condition）
    • dt = 0.1

按照介质中波速的计算公式，手动计算得到的结果为：

$wave_{fluid}=\sqrt{\frac{Bulk}{\rho}}=\sqrt{\frac{220000}{1000}}=14.8324 m/s=14832.4mm/s$

$wave_{gas}=\sqrt{\frac{Bulk}{\rho}}=\sqrt{\frac{120000}{12.6}}=14.8324 m/s=97.59m/s=97590mm/s$

$wave_{float}=\sqrt{\frac{E}{\rho}}=\sqrt{\frac{1.5e12}{6000}}=1.58114e4 m/s=1.58114e7mm/s$

$wave_{strip}=\sqrt{\frac{E}{\rho}}=\sqrt{\frac{1.5e9}{3500}}=654.64 m/s=654654mm/s$

code中进行了单位转换，实际用于计算的单位为mm，并且在计算浮筒和绑带波速时，杨氏模量乘以了相应的thickness=0.01，因此：

$wave_{fluid}=\sqrt{\frac{Bulk}{\rho}}=\sqrt{\frac{220000}{1000}}=14.8324 m/s=14832.4mm/s$

$wave_{gas}=\sqrt{\frac{Bulk}{\rho}}=\sqrt{\frac{120000}{12.6}}=14.8324 m/s=97.59m/s=97590mm/s$

$wave_{float}=\sqrt{\frac{E}{\rho}}=\sqrt{\frac{1.5e12 \times 0.01}{6000}}=1.58114e3 m/s=1.58114e6mm/s$

$wave_{strip}=\sqrt{\frac{E}{\rho}}=\sqrt{\frac{1.5e9 \times 0.01}{3500}}=65.464 m/s=65465.4mm/s$

这与cmdline输出的wavespeed一致，即：

1
2
3
Minimal element size: 0.0124994 (m)
Time step (s) : 7.90533e-07
Wavespeed: fluid=14832.4; airplane=25000; gas=97590; membrane=25000; float=1.58114e+06; strip=65465.4//此处的Wavespeed的单位是mm/s

现在按cmdline输出的参数来计算code里的timestep：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
double unit_factor = 1.0e3;
Cs *= unit_factor; // Cs = 25 * 1000 = 25000 mm/s

double C_fluid = sqrt(Bulk_fluid/rho_fluid); //14832.4 mm/s
double C_gas = sqrt(Bulk_gas/rho_gas); //97590 mm/s
double C_float = sqrt(E_float/rho_float); //1.58114e6 mm/s
double C_strip = sqrt(E_strip/rho_strip); //65465.4 mm/s

double wavespeed = C_fluid > Cs ? C_fluid : Cs;	// wavespeed = Cs = 25000 mm/s
if ( pFloat != NULL ) {
  wavespeed = wavespeed > C_gas ? wavespeed : C_gas; // wavespeed = C_gas =97590 mm/s
  wavespeed = wavespeed > C_float ? wavespeed : C_float; // wavespeed = C_float =1.58114e6 mm/s
  wavespeed = wavespeed > C_strip ? wavespeed : C_strip; // wavespeed = wavespeed =1.58114e6 mm/s
}

dt *= hc / wavespeed; // dt = 0.1, hc =12.4994 mm

按照上述波速：

$wave_{float}>wave_{gas}>wave_{strip}>wave_{airplane}=wave_{membrane}>wave_{fluid}$

$d_t = dt \times \frac{hc}{wave_{float}} = 0.1 \times \frac{0.0124994 \times 1000}{1.58114e6}=7.90533 \times 10^{-7}$

这与cmdline输出的结果是一致的。

问题

存在的问题是计算浮筒和绑带的波速时为何要将其杨氏模量乘以thickness，不乘以thickness的话它们的波速应该大10倍：$wave_{float}=1.58114e7mm/s$，$wave_{strip}=654654mm/s$，进而时间步会缩小10倍：

$$
d_t = dt \times \frac{hc}{wave_{float}} = 0.1 \times \frac{0.0124994 \times 1000}{1.58114e7}=7.90533 \times 10^{-8}
$$

  • thickness的作用是什么？
  • code中单位换算、时间步计算以及cmdline中输出信息的代码实现比较乱；
  • 不过感觉时间步不是主要的问题，从测试结果来看，浮筒材料和绑带材料的变形还不太正确，可能需要将绑带的材料属性增加进去。