SPECFEM手册的一些翻译（Chapter 4）

scouttttt

已于 2024-09-23 22:34:24 修改

阅读量1.1k

点赞数 9

文章标签： github

于 2024-09-23 20:06:33 首次发布

本文链接：https://blog.csdn.net/scouttttt/article/details/142461692

版权

要运行求解器，在主工作目录下输入以下命令：

bin/xspecfem2D
（如果你使用了并行支持的编译，可以使用 mpirun 或等效命令）。

运行后，求解器将在 OUTPUT_FILES/ 目录中输出地震记录和不同时间步长的波前快照。

可视化输出结果：

要查看波场的 Postscript 文件，输入：
```
gs OUTPUT_FILES/vect.ps*
```
在这些文件中，波场用小箭头表示，流体/固体的匹配界面用粗粉线表示，吸收边界用粗绿线表示。
要查看彩色波场快照，输入：
gimp OUTPUT_FILES/image.gif*
这将显示你选定输出的波场的像素化图像（或压力场，取决于你在 DATA/Par_file 中的选择）。

在运行求解器时，请注意以下几点：

提供的 DATA/Par_file 文件可以正常运行，测试代码时无需修改。
地震记录 OUTPUT_FILES/*.sem* 是简单的 ASCII 文件，包含两列数据：第一列是时间，第二列是振幅。你可以使用任何工具来可视化这些数据，例如 "gnuplot"。如果你希望以 Seismic Unix 格式输出二进制地震记录，可以使用参数 SU_FORMAT。在这种情况下，所有地震记录将写入一个带有 .bin 扩展名的文件。要查看头文件信息，可以使用以下命令（取决于你的 Seismic Unix 安装）：
```
surange < Uz_file_single.bin
suoldtonew < Uz_file_single.bin | surange
```
要查看地震记录的波形图，请将 surange 替换为 suxwigb。
如果 DATA/Par_file 中的 MODEL 标志设置为 default，则速度和密度模型通过 nbmodels 和 nbregions 来确定。否则，nbmodels 的值将被忽略，速度和密度模型由用户提供的文件或子程序确定。
使用 Intel 的 ifort 编译时，使用 -assume byterecl 选项来创建显示波场的二进制 PNM 图像。
目录 utils/ 中包含一些有用的脚本和 Fortran 程序。
你可以在 EX2DDIR 中找到用于计算简单介质的解析解的 Fortran 代码，许多文章中的基准测试都使用了该代码作为参考。此代码在 Berg et al. [1994] 中有描述。

关于 `DATA/Par_file` 参数的说明

DATA/Par_file 文件包含详细的注释，几乎是自解释的。生成网格时，请参考相应的说明。以下是一些影响求解器的参数的更详细说明：

1. ATTENUATION_VISCOELASTIC 或 ATTENUATION_VISCOACOUSTIC

关于衰减（粘弹性和粘声学性），在 Par_file 中需要选择标准线性固体（N_SLS）的数量，以模拟恒定的 Q 品质因数。通常使用 N_SLS = 3 是安全的。如果（仅在你确定知道自己在做什么时）想降低模拟成本，可以尝试减少这个值。
图 4.2 提供了一些可供参考的值（再次提醒，只有在你确定知道自己在做什么时才使用这些值）。该表格由 Zhinan Xie 创建，通过比较真实恒定 Q 与基于 N_SLS 标准线性固体近似的结果得出。比较依据是 Kristeková et al. [2009] 提出的时频失配和拟合优度标准。该表格针对一个无量纲参数绘制，表示传播距离。

2. USE_TRICK_FOR_BETTER_PRESSURE

这个选项仅适用于所有接收器记录压力且位于声学元素中的情况。此选项通过使用一种技巧来提高流体（声学）元素中压力地震记录的准确性：
- 使用源的二阶导数作为源时间函数，而不是直接使用源本身，然后记录 potential_acoustic() 作为压力地震记录，而不是 potential_dot_dot_acoustic()。
- 从数学上讲，这两者是等价的，但在数值上，使用 potential_acoustic() 更加准确，因为在显式 Newmark 时间方案中，加速度是零阶精度，而位移是二阶精度。因此，在流体元素中，potential_dot_dot_acoustic() 仅是零阶精度，而 potential_acoustic() 是二阶精度，因而含有显著更少的数值噪声。

这两个参数对模拟的影响较大，特别是在涉及粘弹性或压力地震记录的情况下，建议根据需要谨慎调整这些参数。

关于 `DATA/SOURCE` 参数的说明

DATA/ 目录中的 SOURCE 文件应按如下方式编辑：

source_surf：将此标志设置为 .true.，以强制源位于模型的表面；否则源将位于介质内部。
xs 和 zs：源的位置，xs 为 x 方向的位置（单位：米），zs 为 z 方向的位置（单位：米）。
source_type：设置源类型：
- 1 表示弹性力或声学压力源。
- 2 表示矩张量源。
- 对于包含自由表面反射和转换波的平面波：
  - P 波 = 1
  - S 波 = 2
  - 瑞利波 = 3
- 对于不包含自由表面反射和转换波的平面波（仅入射波）：
  - P 波 = 4
  - S 波 = 5
- 入射平面波由 DATA/Par_file 中的 initialfield 参数启用。
time_function_type：选择源时间函数类型：
- 1 使用 Ricker 波（高斯函数的二阶导数）。
- 2 使用高斯函数的一阶导数。
- 3 使用高斯函数。
- 4 使用 Dirac 函数。
- 5 使用 Heaviside 源时间函数。
Ricker 波形的标准定义为：
Ricker(t)=(1−2at2)e−at2\text{Ricker}(t) = (1 - 2at^2)e^{-at^2}Ricker(t)=(1−2at2)e−at2
其中 a=π2f02a = \pi^2 f_0^2a=π2f02，该波形的傅里叶变换如图 4.3 所示。
f0：设置源的主频率。对于使用 Heaviside 源时间函数的点源模拟（time_function_type = 5），建议将频率参数 f0 设置为较高值，这相当于模拟一个阶跃源时间函数，即一个 δ 函数的矩率函数。

源的半持续时间由 1/f0 给出。如果使用高斯源时间函数（time_function_type = 3），源时间函数的半宽度等于其半持续时间。通常我们在模拟中设置半持续时间为 0，并在后处理时卷积合成地震图，以方便使用多种源时间函数。
t0：对于单一源，建议将时间移位参数 t0 设置为 0.0。时间移位可以在后处理时应用，但对于多个源，t0 可以设置为非零值。
anglesource：源的角度（仅适用于力源）；对于平面波，这表示入射角。对于矩张量源，该参数不适用。
Mxx, Mzz, Mxz：矩张量分量（仅适用于矩张量源，source_type = 2）。请注意，在 SPECFEM2D 和 SPECFEM3D 中，矩张量分量的单位不同：
- 在 SPECFEM3D 中，单位是 dyne*cm。
- 在 SPECFEM2D 中，单位是 N*m。
要将 Strike / Dip / Rake 转换为 CMTSOLUTION 矩张量格式，可以使用 SPECFEM3D_GLOBE 提供的两个小型 C 程序：
- strike_dip_rake_to_CMTSOLUTION.c
- CMTSOLUTION_to_AkiRichards.c
但是在 2D 平面应变 P-SV 模型中，这实际上是一个线源，而不是 3D 点源。有关更多详细信息，请参考 Pilant [1979] 的书中的第 7.3 节“二维点源”，以及 Helmberger 和 Vidale [1988] 的研究。
factor：放大因子。

模拟的零时间点

模拟的零时间对应于三角形/高斯函数的中心，或地震事件的质心时间。模拟的起始时间为：
t=−1.2×half duration+t0t = -1.2 \times \text{half duration} + t0t=−1.2×half duration+t0
其中，因子 1.2 用于确保模拟开始时，矩率函数接近于零。对于 Heaviside 函数，这个因子为 2.0。半持续时间由 1/f0 给出。如果你希望手动设置起始时间，可以在 constants.h 文件中将参数 USER_T0 设置为一个正的非零值，在这种情况下，模拟将在 -USER_T0 的时间点开始。

4.1 如何运行弹性波模拟

对于各向同性的弹性材料，有两种选择：

P-SV 波：
- 如果要计算在 x-z 平面传播的 P-SV 波（纵波和剪切波），请在 Par_file 中将 p_sv 设置为 .true.。
SH 波：
- 如果要计算在 x-z 平面传播并具有 y 方向运动分量的 SH 波（膜波），请将 p_sv 设置为 .false.。

此功能仅对弹性材料实现，并且可以计算敏感性核。有关膜表面波的详细信息，请参阅 Tape 等人 [2007] 的研究。

此外，还提供了一个名为 SEM_save_dir.py 的可选 Python 脚本。该脚本允许自动保存给定模拟的所有参数和结果，非常实用。

4.2 如何运行轴对称波模拟

在 SPECFEM2D 中可以进行轴对称模拟。对于这些模拟，计算的二维域物理上是轴对称三维域的子午面。作为入门，建议阅读我们发表的文章 [Bottero et al., 2016]。要将几何设置为轴对称，请在 Par_file 文件中将 AXISYM 标志设置为 .true.：

AXISYM = .true.

此时模型的左边界变为对称轴，计算出的波场物理上是通过将二维波场围绕左边界旋转获得的三维波场。

关于源的注意事项：

在轴对称几何中，整个模型（包括源）相对于对称轴是对称的。因此，如果源不在对称轴上，它物理上会具有圆形形状。这在一些应用中（如无损检测）是可行且有意义的，但在大多数情况下是不需要的。需要牢记这一点。在声学介质中，液体中的爆炸自然是轴对称的，生成的波场具有正确的三维形状。然而，如果源位于弹性固体中，其三维辐射图案将是轴对称的。

快速入门：

在 EXAMPLES/axisymmetric_case_AXISYM_option 目录下提供了一个简单的示例，建议阅读其中的 README 文件。该示例展示了如何使用 AXISYM 选项，并通过使用半解析代码 OASES（Schmidt [2004]）进行验证。在该示例中，研究的域为位于粘弹性介质上方的一层水，源是水中的爆炸，并且该域由 PML（完美匹配层）边界包围。

关于外部网格生成器：

在轴对称几何中可以使用外部网格生成器。EXAMPLES/paper_axisymmetry_example 目录提供了一个 Cubit/Trelis 网格生成器的示例（详见 http://www.csimsoft.com/trelis）。建议参考该示例并阅读前一章了解更多细节。

与平面应变几何的唯一不同之处在于，SPECFEM2D 需要一个额外的文件来定义轴元素。这个文件的路径必须在 Par_file 中指定：

axial_elements_file = /path/to/the/axial_elements_file

轴元素文件的结构如下：

48
1 2 8456 8457
2 2 8457 8458
3 2 8458 8459
4 2 8459 8460
...

第一行包含轴元素的数量，接下来的行包含四列：元素 ID、轴元素描述的节点数（始终为 2）、第一个节点 ID、第二个节点 ID。请注意，轴元素必须包括与轴接触的可能的 PML（上/下）元素。为简化起见，假设与对称轴接触的网格元素是通过完整边接触的，而不是通过单个点接触的，这意味着左侧最外层的网格元素必须是结构化的，其他部分可以是非结构化的。

关于分辨率的注意事项：

在轴对称模拟中，轴元素使用不同的求积法，这使得每个波长所需的点数比平面应变模拟稍大（约 25%）。

关于小错误的注意事项：

代码中仍然存在一个小错误，虽然输出的信号在整个域中是正确的，但包含源的元素中存在异常。这个小错误尚未解决，但不影响代码的正常使用。

学习示例代码：

在 utils/small_SEM_solver_in_Fortran_without_MPI_to_learn 目录中提供了一个简化的演示代码，用于学习光谱元方法在平面应变和轴对称几何中的工作原理。如果你感兴趣，可以查看该代码。在该目录中输入：

./make_Fortran_2D_axisymmetric.csh

进行编译，然后输入 ./xspecfem2D 运行。上述讨论的小错误不在该代码中。

4.3 如何运行各向异性波模拟

在 SPECFEM2D 中，我们遵循 Carcione 等人 [1988a] 的方法，使用经典的简化 Voigt 符号来表示对称张量 [Helbig, 1994, Carcione, 2007]。各向异性弹性固体的本构关系可以用广义胡克定律表达，公式为：

$\sigma_{ij} = c_{ijkl} \varepsilon_{kl}, \quad i, j, k = 1, \dots , 3$

矩阵符号表示：

为了表达横向各向同性介质的应力-应变关系，我们引入了在文献中常用的简化矩阵符号。通过该约定，弹性张量的成对下标由单个数字替换，具体如下：

$(11) \to 1, \quad (22) \to 2, \quad (33) \to 3,$

$\quad (23) = (32) \to 4, \quad (31) = (13) \to 5, \quad (12) = (21) \to 6$

因此，在最一般的二维情况下，应力-应变关系的表示如下：

$\sigma_{xx} = c_{11} \frac{\partial u_x}{\partial x} + c_{13} \frac{\partial u_z}{\partial z} + c_{15} \left( \frac{\partial u_z}{\partial x} + \frac{\partial u_x}{\partial z} \right)$

$\sigma_{zz} = c_{13} \frac{\partial u_x}{\partial x} + c_{33} \frac{\partial u_z}{\partial z} + c_{35} \left( \frac{\partial u_z}{\partial x} + \frac{\partial u_x}{\partial z} \right)$

$\sigma_{xz} = c_{15} \frac{\partial u_x}{\partial x} + c_{35} \frac{\partial u_z}{\partial z} + c_{55} \left( \frac{\partial u_z}{\partial x} + \frac{\partial u_x}{\partial z} \right)$