PDF
Acoustic Reflections off a Water-Sediment Interface
Introduction
This validation model computes the reflection coefficient of acoustic waves off a water-sediment interface. Homogeneous waves are incident from a fluid (water) domain and are reflected and transmitted at a water-sediment interface. The sediment domain is modeled using Biot’s theory with the use of the Poroelastic Waves interface. It includes the detailed interaction between the pressure waves in the saturating fluid and the elastic waves in the porous matrix (the sediment). The model results are compared with those obtained by Stoll and Kan (Ref. 1). A visual comparison of the results with those presented in Ref. 1 shows good agreement.
Figure 1: Sketch of the water-sediment system. The computational unit cell has the width W and the height H with perfectly matched layers at the top and bottom. λ is the acoustic wavelength in water. An acoustic wave with wave vector k is incident at an angle θ0 on the water-sediment interface.
Knowledge about the plane wave reflection and refraction coefficients are important when interpreting acoustic data taken from marine seismology. The kind of model studied here can provide insight into the different phenomena at play. The sediment is modeled using Biot’s theory with the addition of the porous matrix skeleton having viscoelastic properties. This is achieved by simply allowing the bulk modulus and the shear modulus to have a small imaginary component (a constant attenuation factor). Inside the porous material, three types of waves can propagate: a fast and a slow pressure wave, and a shear wave. The interplay between these waves and the incoming homogeneous pressure wave in the fluid results in nontrivial behavior of the reflection coefficient. See Ref. 1 for further discussion.
Model Definition
A sketch of the modeled system is shown in Figure 1. An acoustic wave is incident on the water-sediment interface at an angle θ0, defining the wave vector
where ω is the angular frequency 2πf. The wave vector is used to define an incident plane wave of the form
where x = (x, y). The wave vector is defined as a variable in the model and used to define the direction of the plane wave in the Background Pressure Field feature, as well as the wave vector for the Floquet periodicity in the model. Note that in the nomenclature of the Acoustics Module the incident field is the background acoustic field.
The fluid domain is modeled by classical pressure acoustics (lossless) solving the Helmholtz equation. The sediment domain is modeled by Biot’s theory solving for both the pressure and the displacement field of the porous matrix. The coupled multiphysics problem is set up using the Acoustic-Poroelastic Waves Interaction multiphysics interface.
The problem is 2D and is modeled as a periodic structure using Floquet periodic conditions on both the acoustic and the porous domains. To ensure that the problem is well defined numerically, the domain size is set to depend on the wavelength in water (thus the frequency). The width of the computational unit cell W is equal to the wavelength and the height H is equal to two times the wavelength. The computational domain is truncated, in both the fluid and in the porous domain, using perfectly matched layers (PMLs). The thickness is set equal to one third of the wavelength.
The material parameters necessary for defining the properties of the porous domain are taken from Ref. 1 and some are calculated using these values. The porous material is here assumed to be sand. The most important parameters and their definitions are given in Table 1 below.
Table 1: Porous Material Parameters for Sand.
Symbol
Value
Description
 µF
1·10-3 Pa·s
Fluid viscosity (water)
 ρF
1000 kg/m3
Fluid density (water)
 εP
0.47
Porosity
 κP
1·10-6 cm2
Permeability
 Kb
4.36·107 Pa
Bulk modulus of frame (drained porous matrix)
 Ks
3.6·1010 Pa
Bulk modulus of solid grains (constituting the frame)
 αB
1−Kb/Ks = 0.999
Biot-Willis coefficient: Because it is very close to 1 the frame is said to be limp.
 logdec
0.15
Logarithmic decrement factor
 Kbc
(1+i·logdec/π)·Kb
Complex bulk modulus of frame
 G
2.61·107 Pa
Shear modulus of frame (drained porous matrix)
 Gc
(1+i·logdec/π)·G
Complex shear modulus of frame
 ρs
2650 kg/m3
Solid density (of material constituting the frame)
 ρd
(1-εP)·ρs
= 1404.5 kg/m3
Drained density (of the porous matrix)
 τ
1.25
Tortuosity (chosen, not specified in Ref. 1)
 a
4·10-3 cm
Pore size parameter (chosen to have fc = 100 Hz)
 fr
µF/(2π·ρF·a2)
= 99.5 Hz
Reference frequency: frequency at which the viscous boundary layer thickness (viscous penetration depth) is equal to the pore size parameter a. This value is automatically calculated by COMSOL
Note: Further information about the physics solved by the Poroelastic Waves interface, about Biot’s theory, and about the meaning of the material parameters is given in the documentation. Press Ctrl+F1 and browse to the Acoustics Module and the User’s Guide. The relevant theory is found in Elastic Waves Interfaces>Theory for the Poroelastic Waves Interfaces.
Results and Discussion
The reflection coefficient R and the absorption coefficient α are plotted for different combinations of the angle of incidence and the frequency (the model parameters). The reflection coefficient is the complex ratio of the reflected (scattered) pressure to the incident (background) pressure, and is given by
(1)
where the subscript re stands for reflected and tot for total. The absorption coefficient α gives a measure of the relative absorbed energy. It is in general defined through the incident power Pin and the reflected power Pre, through
(2)
Note that for plane waves (when no higher oder diffraction modes exist) the absorption coefficient is given by
(3)
The absolute value of the reflection coefficient R is shown in Figure 2 and Figure 3. Note that average of the reflection coefficient across the interface is used. This is because the model extracts the effective homogenized properties of the water-sediment interface. In the first figure, R is depicted as function of the angle of incidence θ0 and in the second figure as function of the frequency f. The absorption coefficient α, corresponding to the reflection coefficient in Figure 2, is shown in Figure 4. An interesting feature that can be deduced from Figure 3 is that the reflection at the interface acts as a filter with respect to a broadband signal; this is especially true at high angles of incidence.
Figure 2: Absolute value of the reflection coefficient R as function of the angle of incidence for the driving frequencies 10 Hz, 100 Hz, 1 kHz, and 10 kHz.
Figure 3: Absolute value of the reflection coefficient R as function of frequency for the angles of incidence 30°, 40°, 50°, 60°, 70°, and 80°.
Figure 4: Absorption coefficient α as function of the angle of incidence for the driving frequencies 10 Hz, 100 Hz, 1 kHz, and 10 kHz.
As can be seen from Figure 2, the reflection coefficient has significantly different behavior for 3 of the driving frequencies at an incidence angle of, for example, 70°. The total pressure is shown in Figure 5 and the displacement field is shown in Figure 6. Note the different length scales as the computational domain (geometry) scales with the wavelength.
From the graphs in Figure 2 and Figure 4, it is evident that the behavior shifts for the different frequencies for increasing angles of incidence. At 90°, the reflection coefficient is of course 1 and the absorption 0.
Figure 5: Total pressure distribution for an incidence angle of 70° for 10 Hz (top left), 100 Hz (top right), 1 kHz (bottom left), and 10 kHz (bottom right). Note that the computational domain size (geometry) scales with the wavelength.
Figure 6: Displacement amplitude with deformation of the porous matrix for an incidence angle of 70° for 10 Hz (top left), 100 Hz (top right), 1 kHz (bottom left), and 10 kHz (bottom right). Note that the computational domain size (geometry) scales with the wavelength.
Notes About the COMSOL Implementation
•
In order for the model to be numerically well behaved for all frequencies studied, the geometry of the periodic unit cell is set to depend on the wavelength lambda0. This means that the geometry depends on the frequency. To model this, a parametric sweep which updates the geometry is added. The sweep parameter is the frequency f0, which is then entered as the study frequency in the Frequency Domain study step.
•
The PMLs do not automatically compensate for the oblique angle of incidence of the waves. It is recommended to use a PML stretching factor of 1/cos(theta) in situations with plane waves at an oblique incidence.
•
Rational PML stretching is unsuitable for evanescent wave components such as the fast wave above the critical angle of total internal reflection, so it is recommended to use polynomial stretching in the sediment domains.
•
The periodic plots shown in Figure 5 and Figure 6 are set up using a Periodic 2D dataset. Using the Floquet-Bloch periodicity option and entering the wave vector components k0x and k0y with this option, the correct phase is automatically added to the underlying solution variables.
•
From the material parameters in Table 1, it is seen that the Biot-Willis coefficient αB is close to 1. This means that the porous matrix is in a so-called limp configuration. At the other extreme, when αB is close to the value of porosity εP, the configuration is named rigid (the skeleton is not moving). When a porous material is in one of these states it is possible to use an equivalent fluid model from the Poroacoustics domain feature (found in the Pressure Acoustics, Frequency Domain interface) to model the porous material. Use, for example, the Johnson-Champoux-Allard (limp frame option) for the case studied here. The advantage of using this model is that it is much less computationally expensive than the full Biot model. An equivalent fluid model generates or mimics the losses in the porous domain and therefore predicts the reflection and absorption coefficient. However, it does not account for the details of the movement of the skeleton and the propagation of the three different types of waves. Moreover, as soon as a porous material is neither limp nor rigid, the full Biot model of the Poroelastic Waves interface is necessary.
Reference
1. R.D. Stoll and T.-K. Kan, “Reflections of Acoustic Waves at a Water-Sediment Interface”, J. Acoust. Soc. Am., vol. 70, no. 1, pp. 149–156, 1981.
Application Library path: Acoustics_Module/Underwater_Acoustics/reflections_water_sediment
Modeling Instructions
Start by adding the Acoustic-Poroelastic Waves Interaction multiphysics interface. The interface consists of the Pressure Acoustics, Frequency Domain interface, the Poroelastic Waves interface, and the necessary Multiphysics coupling.
From the File menu, choose New.
New
In the New window, click  Model Wizard.
Model Wizard
1
In the Model Wizard window, click  2D.
2
In the Select Physics tree, select Acoustics>Acoustic-Structure Interaction>Acoustic-Poroelastic Waves Interaction.
3
Click Add.
4
Click  Study.
5
In the Select Study tree, select General Studies>Frequency Domain.
6
Click  Done.
Geometry 1
First, load the material parameters from a file. The list contains definitions of material parameters (see Table 1).
Global Definitions
Parameters - Material
1
In the Model Builder window, under Global Definitions click Parameters 1.
2
In the Settings window for Parameters, type Parameters - Material in the Label text field.
3
Locate the Parameters section. Click  Load from File.
4
Browse to the model’s Application Libraries folder and double-click the file reflections_water_sediment_material_parameters.txt.
Parameters - Model
1
In the Home toolbar, click  Parameters and choose Add>Parameters.
This second parameter group loaded from a file contains geometry lengths, and the two model parameters (frequency f0 and angle of incidence theta0).
2
In the Settings window for Parameters, type Parameters - Model in the Label text field.
3
Locate the Parameters section. Click  Load from File.
4
Browse to the model’s Application Libraries folder and double-click the file reflections_water_sediment_model_parameters.txt.
Build the geometry, a sketch is given in Figure 1.
Geometry 1
Rectangle 1 (r1)
1
In the Geometry toolbar, click  Rectangle.
2
In the Settings window for Rectangle, locate the Size and Shape section.
3
In the Width text field, type W.
4
In the Height text field, type H.
5
Click to expand the Layers section. Clear the Layers on bottom check box.
6
Select the Layers on top check box.
7
In the table, enter the following settings:
Layer name
Thickness (m)
Layer 1
Hpml
8
Click  Build Selected.
Rectangle 2 (r2)
1
In the Geometry toolbar, click  Rectangle.
2
In the Settings window for Rectangle, locate the Size and Shape section.
3
In the Width text field, type W.
4
In the Height text field, type H.
5
Locate the Position section. In the y text field, type -H.
6
Locate the Layers section. In the table, enter the following settings:
Layer name
Thickness (m)
Layer 1
Hpml
7
Click  Build All Objects.
8
Click the  Zoom Extents button in the Graphics toolbar.
Under the Definitions node add variables as well as three nonlocal couplings. The couplings are used to pick up variable values. One samples a point in the water domain and the others evaluates the average/integral over the water-sediment interface.
Definitions
Variables 1
1
In the Model Builder window, under Component 1 (comp1) right-click Definitions and choose Variables.
2
In the Settings window for Variables, locate the Variables section.
3
Click  Load from File.
4
Browse to the model’s Application Libraries folder and double-click the file reflections_water_sediment_variables.txt.
Integration 1 (intop1)
1
In the Definitions toolbar, click  Nonlocal Couplings and choose Integration.
2
In the Settings window for Integration, locate the Source Selection section.
3
From the Geometric entity level list, choose Point.
4
Select Point 4 only.
Average 1 (aveop1)
1
In the Definitions toolbar, click  Nonlocal Couplings and choose Average.
2
In the Settings window for Average, locate the Source Selection section.
3
From the Geometric entity level list, choose Boundary.
4
Select Boundary 6 only.
Integration 2 (intop2)
1
In the Definitions toolbar, click  Nonlocal Couplings and choose Integration.
2
In the Settings window for Integration, locate the Source Selection section.
3
From the Geometric entity level list, choose Boundary.
4
Select Boundary 6 only.
Now, create selections to be used when setting up the physics.
Porous Domain
1
In the Definitions toolbar, click  Explicit.
2
In the Settings window for Explicit, type Porous Domain in the Label text field.
3
Select Domains 1 and 2 only.
Water Domain
1
In the Definitions toolbar, click  Explicit.
2
In the Settings window for Explicit, type Water Domain in the Label text field.
3
Select Domains 3 and 4 only.
Proceed to set up the physics and boundary conditions of the problem. The condition at the water-sediment interface is defined under the Multiphysics node.
The properties of the porous material are edited after adding the porous material.
Poroelastic Waves (pelw)
1
In the Model Builder window, under Component 1 (comp1) click Poroelastic Waves (pelw).
2
In the Settings window for Poroelastic Waves, locate the Domain Selection section.
3
From the Selection list, choose Porous Domain.
4
Locate the Sound Pressure Level Settings section. From the Reference pressure for the sound pressure level list, choose Use reference pressure for water.
Periodic Condition 1
1
In the Physics toolbar, click  Boundaries and choose Periodic Condition.
2
Select Boundaries 1, 3, 10, and 11 only.
3
In the Settings window for Periodic Condition, locate the Periodicity Settings section.
4
From the Type of periodicity list, choose Floquet periodicity.
5
Specify the kF vector as
k0x
X
k0y
Y
Pressure Acoustics, Frequency Domain (acpr)
1
In the Model Builder window, under Component 1 (comp1) click Pressure Acoustics, Frequency Domain (acpr).
2
In the Settings window for Pressure Acoustics, Frequency Domain, locate the Domain Selection section.
3
From the Selection list, choose Water Domain.
4
Locate the Sound Pressure Level Settings section. From the Reference pressure for the sound pressure level list, choose Use reference pressure for water.
5
In the Physics toolbar, click  Boundaries and choose Periodic Condition.
Periodic Condition 1
1
Select Boundaries 5, 7, 12, and 13 only.
2
In the Settings window for Periodic Condition, locate the Periodicity Settings section.
3
From the Type of periodicity list, choose Floquet periodicity.
4
Specify the kF vector as
k0x
x
k0y
y
Background Pressure Field 1
1
In the Physics toolbar, click  Domains and choose Background Pressure Field.
2
In the Settings window for Background Pressure Field, locate the Domain Selection section.
3
From the Selection list, choose Water Domain.
4
Locate the Background Pressure Field section. In the p0 text field, type 1.
5
From the c list, choose From material.
6
Specify the ek vector as
kx_e
x
ky_e
y
7
Select the Calculate background and scattered field intensity check box.
8
From the ρ list, choose From material.
Materials
Water
1
In the Model Builder window, under Component 1 (comp1) right-click Materials and choose Blank Material, twice, to create two new blank materials.
2
In the Settings window for Material, type Water in the Label text field.
Sediment
1
In the Model Builder window, under Component 1 (comp1)>Materials click Material 2 (mat2).
2
In the Settings window for Material, type Sediment in the Label text field.
Now, go to the Poroelastic Material node and specify the porous matrix material parameters to be taken from the Sediment material that you just created. When done, return to the material, where it is now evident which parameters to define.
Poroelastic Waves (pelw)
Poroelastic Material 1
1
In the Model Builder window, under Component 1 (comp1)>Poroelastic Waves (pelw) click Poroelastic Material 1.
2
In the Settings window for Poroelastic Material, locate the Porous Matrix Properties section.
3
From the Porous elastic material list, choose Sediment (mat2).
4
Locate the Fluid Properties section. From the Viscosity model list, choose Biot’s high frequency range.
5
From the Specify list, choose Characteristic pore size.
6
In the a text field, type a.
Materials
Water (mat1)
1
In the Model Builder window, under Component 1 (comp1)>Materials click Water (mat1).
2
In the Settings window for Material, locate the Material Contents section.
3
In the table, enter the following settings:
Property
Variable
Value
Unit
Property group
Density
rho
rhoF
kg/m³
Basic
Speed of sound
c
c0
m/s
Basic
Dynamic viscosity
mu
muF
Pa·s
Basic
Compressibility of fluid
chif
1/Kf
1/Pa
Basic
Sediment (mat2)
1
In the Model Builder window, click Sediment (mat2).
2
In the Settings window for Material, locate the Material Contents section.
3
In the table, enter the following settings:
Property
Variable
Value
Unit
Property group
Bulk modulus
K
Kbc
N/m²
Bulk modulus and shear modulus
Shear modulus
G
Gc
N/m²
Bulk modulus and shear modulus
Density
rho
rhod
kg/m³
Basic
Porosity
epsilon
epsilonP
1
Basic
Permeability
kappa_iso ; kappaii = kappa_iso, kappaij = 0
kappaP
m²
Basic
Biot-Willis coefficient
alphaB
alphaB0
1
Poroelastic material
Tortuosity factor
tau
tau0
1
Poroacoustics model
Multiphysics
Acoustic-Porous Boundary 1 (apb1)
Now is a good time to inspect the Multiphysics node and look at the multiphysics coupling. Click on Acoustic-Porous Boundary 1. You can see that it is active on the interface between the water domain and the porous sediment.
The last step before building the mesh is to set up and define the perfectly matched layers (PMLs). Define one for the poroelastic domain and another for the water domain.
Definitions
Perfectly Matched Layer 1 (pml1)
1
In the Definitions toolbar, click  Perfectly Matched Layer.
2
Select Domain 1 only.
Switch the Physics to Poroelastic Waves. The polynomial stretching is generally most appropriate when there is a mix of different wave types in the model, as it is the case for the Poroelastic Waves domain.
3
In the Settings window for Perfectly Matched Layer, locate the Scaling section.
4
From the Physics list, choose Poroelastic Waves (pelw).
The PMLs do not compensate for the oblique angle of incidence of the waves. Typically, the PML stretching factor should be 1/cos(theta) in situations with plane waves at oblique incidence.
5
In the PML scaling factor text field, type if(theta0==90[deg],1e6,1/cos(theta0)).
Perfectly Matched Layer 2 (pml2)
1
In the Definitions toolbar, click  Perfectly Matched Layer.
2
Select Domain 4 only.
3
In the Settings window for Perfectly Matched Layer, locate the Scaling section.
4
From the Coordinate stretching type list, choose Rational.
The PMLs do not compensate for the oblique angle of incidence of the waves. Typically, the PML stretching factor should be 1/cos(theta) in situations with plane waves at oblique incidence.
5
In the PML scaling factor text field, type if(theta0==90[deg],1e6,1/cos(theta0)).
Keep the default setting in the Physics list, that is Pressure Acoustics, Frequency Domain (acpr).
Build the mesh by following the instructions below. When modeling systems with different materials it is important to have a mesh that solves for the slowest wave traveling through each of the materials. In this particular example, where the background field could travel in any direction, it is important to have enough elements in the PML domain in both directions. It is also important for models with periodic conditions to have the same mesh at source and destination. Ideally, a mesh sensitivity analysis should be made for models of this kind. In this model, the mesh is set up manually. Proceed by directly adding a mapped mesh in order to reduce the computational cost of the analysis.
Mesh 1
Mapped 1
In the Mesh toolbar, click  Mapped.
Size
1
In the Model Builder window, click Size.
2
In the Settings window for Size, locate the Element Size section.
3
Click the Custom button.
In general, 5 to 6 second-order elements per wavelength are needed to resolve the waves. For more details, see Meshing (Resolving the Waves) in the Acoustics Module User’s Guide. In this model we use 5 elements per wavelength. Use the wavelength in water in the water domain.
4
Locate the Element Size Parameters section. In the Maximum element size text field, type lam0/5.
5
In the Minimum element size text field, type lam0/10.
Mapped 1
1
In the Model Builder window, click Mapped 1.
2
In the Settings window for Mapped, locate the Domain Selection section.
3
From the Geometric entity level list, choose Domain.
4
Click in the Graphics window and then press Ctrl+A to select all domains.
Size 1
1
Right-click Mapped 1 and choose Size.
Use the wavelength of the slowest wave in the poroelastic waves domains.
2
In the Settings window for Size, locate the Geometric Entity Selection section.
3
From the Selection list, choose Porous Domain.
4
Locate the Element Size section. Click the Custom button.
5
Locate the Element Size Parameters section.
6
Select the Maximum element size check box. In the associated text field, type cs_poro/f0/5.
7
Select the Minimum element size check box. In the associated text field, type cs_poro/f0/10.
Add a distribution to make sure that the PML domains have enough number of elements to resolve all the waves.
Distribution 1
1
Right-click Mapped 1 and choose Distribution.
2
Select Boundaries 1 and 7 only.
3
In the Settings window for Distribution, locate the Distribution section.
4
In the Number of elements text field, type 30.
Proceed and add a single boundary layer mesh at the interface between the water and the sediment. This layer will ensure a good evaluation of the gradient involved in the intensity computation.
Boundary Layers 1
1
In the Mesh toolbar, click  Boundary Layers.
2
In the Settings window for Boundary Layers, locate the Domain Selection section.
3
From the Geometric entity level list, choose Domain.
4
From the Selection list, choose Water Domain.
5
Click to expand the Transition section. Clear the Smooth transition to interior mesh check box.
Boundary Layer Properties
1
In the Model Builder window, click Boundary Layer Properties.
2
Select Boundary 6 only.
3
In the Settings window for Boundary Layer Properties, locate the Layers section.
4
In the Number of layers text field, type 1.
5
Click  Build All.
The finalized mesh should look like that in the figure below.
Study 1 - Detailed Angular Sweep
1
In the Model Builder window, click Study 1.
2
In the Settings window for Study, type Study 1 - Detailed Angular Sweep in the Label text field.
3
Locate the Study Settings section. Clear the Generate default plots check box.
Step 1: Frequency Domain
1
In the Model Builder window, under Study 1 - Detailed Angular Sweep click Step 1: Frequency Domain.
2
In the Settings window for Frequency Domain, locate the Study Settings section.
3
In the Frequencies text field, type f0.
4
Click to expand the Study Extensions section. Select the Auxiliary sweep check box.
5
From the Sweep type list, choose All combinations.
6
Click  Add.
7
In the table, enter the following settings:
Parameter name
Parameter value list
Parameter unit
theta0 (Angle of incidence)
range(0,10[deg],50[deg]) range(51[deg],1[deg],89[deg]) 89.5[deg]
deg
Parametric Sweep
1
In the Study toolbar, click  Parametric Sweep.
2
In the Settings window for Parametric Sweep, locate the Study Settings section.
3
Click  Add.
4
In the table, enter the following settings:
Parameter name
Parameter value list
Parameter unit
f0 (Driving frequency)
10 100 1000 10000
Hz
5
In the Study toolbar, click  Compute.
Results
Add a selection to the dataset to only plot results in the physical domains (not in the PMLs). You can study the behavior of the PMLs by removing this selection again or creating a separate dataset without the selection.
1
In the Model Builder window, expand the Results node.
Study 1 - Detailed Angular Sweep/Parametric Solutions 1 (sol2)
In the Model Builder window, expand the Results>Datasets node, then click Study 1 - Detailed Angular Sweep/Parametric Solutions 1 (sol2).
Selection
1
In the Results toolbar, click  Attributes and choose Selection.
2
In the Settings window for Selection, locate the Geometric Entity Selection section.
3
From the Geometric entity level list, choose Domain.
4
Select Domains 2 and 3 only.
Now, also add an array dataset in order to easily plot the result in the periodic structure.
Array 2D 1
1
In the Results toolbar, click  More Datasets and choose Array 2D.
2
In the Settings window for Array 2D, locate the Data section.
3
From the Dataset list, choose Study 1 - Detailed Angular Sweep/Parametric Solutions 1 (sol2).
4
Locate the Array Size section. In the X size text field, type 4.
5
Click to expand the Advanced section. Select the Floquet-Bloch periodicity check box.
6
Find the Wave vector subsection. In the X text field, type k0x.
7
In the Y text field, type k0y.
Now, create four different 2D plots representing the incident pressure (background pressure), the reflected pressure (scattered pressure), the total pressure, and the displacement in the porous matrix.
Incident Pressure
1
In the Results toolbar, click  2D Plot Group.
2
In the Settings window for 2D Plot Group, type Incident Pressure in the Label text field.
3
Click to expand the Title section. From the Title type list, choose Manual.
4
In the Title text area, type Incident Pressure.
5
In the Parameter indicator text field, type f0 = eval(freq) Hz, theta0 = eval(theta0, deg) deg.
Surface 1
1
Right-click Incident Pressure and choose Surface.
2
In the Settings window for Surface, locate the Expression section.
3
In the Expression text field, type acpr.p_b.
4
Locate the Coloring and Style section. Click  Change Color Table.
5
In the Color Table dialog box, select Wave>Wave in the tree.
6
Click OK.
7
In the Settings window for Surface, locate the Coloring and Style section.
8
From the Scale list, choose Linear symmetric.
Arrow Surface 1
1
In the Model Builder window, right-click Incident Pressure and choose Arrow Surface.
2
In the Settings window for Arrow Surface, locate the Expression section.
3
In the X-component text field, type k0x.
4
In the Y-component text field, type k0y.
5
Locate the Arrow Positioning section. Find the X grid points subsection. In the Points text field, type 5.
6
Locate the Coloring and Style section. From the Color list, choose Black.
Incident Pressure
1
In the Model Builder window, click Incident Pressure.
2
In the Settings window for 2D Plot Group, locate the Data section.
3
From the Dataset list, choose Study 1 - Detailed Angular Sweep/Parametric Solutions 1 (sol2).
4
From the Parameter value (freq (Hz),theta0 (deg)) list, choose 4: freq=10000 Hz, theta0=30 deg.
Select any of the desired combinations of the frequency f0 and the angle of incidence theta0 and plot to have a look at the results. You may need to use the Zoom Extent tool as the geometry depends on the frequency.
5
In the Incident Pressure toolbar, click  Plot.
6
Click the  Zoom Extents button in the Graphics toolbar.
Reflected Pressure
1
In the Home toolbar, click  Add Plot Group and choose 2D Plot Group.
2
In the Settings window for 2D Plot Group, type Reflected Pressure in the Label text field.
3
Locate the Data section. From the Dataset list, choose Study 1 - Detailed Angular Sweep/Parametric Solutions 1 (sol2).
4
From the Parameter value (freq (Hz),theta0 (deg)) list, choose 4: freq=10000 Hz, theta0=30 deg.
5
Locate the Title section. From the Title type list, choose Manual.
6
In the Title text area, type Reflected Pressure.
7
In the Parameter indicator text field, type f0 = eval(freq) Hz, theta0 = eval(theta0, deg) deg.
Surface 1
1
Right-click Reflected Pressure and choose Surface.
2
In the Settings window for Surface, locate the Expression section.
3
In the Expression text field, type acpr.p_s.
4
Locate the Coloring and Style section. Click  Change Color Table.
5
In the Color Table dialog box, select Wave>Wave in the tree.
6
Click OK.
7
In the Settings window for Surface, locate the Coloring and Style section.
8
From the Scale list, choose Linear symmetric.
9
In the Reflected Pressure toolbar, click  Plot.
Surface 2
1
In the Model Builder window, right-click Reflected Pressure and choose Surface.
2
In the Settings window for Surface, locate the Expression section.
3
In the Expression text field, type pelw.p_s.
4
Click to expand the Inherit Style section. From the Plot list, choose Surface 1.
5
In the Reflected Pressure toolbar, click  Plot.
Note that this gives the scattered field in the water domain while it is equal to the total acoustic field in the porous domain.
In the next figure, plot the total pressure (water and porous domains). Use the periodic dataset to display the periodic nature of the modeled system.
The Floquet-periodic nature of the solution is recovered with the selected Floquet-Bloch periodicity option in the Array 2D dataset. With this option the correct phase is added to the underlying solution variables.
Total Pressure
1
In the Home toolbar, click  Add Plot Group and choose 2D Plot Group.
2
In the Settings window for 2D Plot Group, type Total Pressure in the Label text field.
3
Locate the Data section. From the Dataset list, choose Array 2D 1.
4
From the Parameter value (freq (Hz),theta0 (deg)) list, choose 4: freq=10000 Hz, theta0=30 deg.
5
Locate the Title section. From the Title type list, choose Manual.
6
In the Title text area, type Total Pressure.
7
In the Parameter indicator text field, type f0 = eval(freq) Hz, theta0 = eval(theta0, deg) deg.
Surface 1
1
Right-click Total Pressure and choose Surface.
2
In the Settings window for Surface, click Replace Expression in the upper-right corner of the Expression section. From the menu, choose Component 1 (comp1)>Multiphysics>apb1.p_t - Total acoustic pressure - Pa.
The multiphysics pressure variable is defined as acpr.p_t in the water domain and pelw.p_t in the sediment domain.
3
Locate the Coloring and Style section. Click  Change Color Table.
4
In the Color Table dialog box, select Wave>Wave in the tree.
5
Click OK.
6
In the Settings window for Surface, locate the Coloring and Style section.
7
From the Scale list, choose Linear symmetric.
8
In the Total Pressure toolbar, click  Plot.
9
Click the  Zoom Extents button in the Graphics toolbar.
Now, plot the displacement of the porous matrix. Again use the periodic dataset to represent the periodic nature of the modeled system, just as for the total pressure.
Displacement
1
In the Home toolbar, click  Add Plot Group and choose 2D Plot Group.
2
In the Settings window for 2D Plot Group, type Displacement in the Label text field.
3
Locate the Title section. From the Title type list, choose Manual.
4
In the Title text area, type Displacement.
5
In the Parameter indicator text field, type f0 = eval(freq) Hz, theta0 = eval(theta0, deg) deg.
6
Locate the Data section. From the Dataset list, choose Array 2D 1.
7
From the Parameter value (freq (Hz),theta0 (deg)) list, choose 4: freq=10000 Hz, theta0=30 deg.
Surface 1
1
Right-click Displacement and choose Surface.
2
In the Settings window for Surface, locate the Expression section.
3
In the Expression text field, type pelw.disp.
Deformation 1
1
Right-click Surface 1 and choose Deformation.
2
In the Settings window for Deformation, locate the Expression section.
3
In the x-component text field, type u.
4
In the y-component text field, type v.
5
Locate the Scale section.
6
Select the Scale factor check box. In the associated text field, type 2e9.
7
In the Displacement toolbar, click  Plot.
8
Click the  Zoom Extents button in the Graphics toolbar.
The total pressure and the displacement are depicted in Figure 5 and Figure 6 for several combinations of the driving frequency and the angle of incidence.
Wave Speed
1
In the Home toolbar, click  Add Plot Group and choose 2D Plot Group.
2
In the Settings window for 2D Plot Group, type Wave Speed in the Label text field.
Surface 1
1
Right-click Wave Speed and choose Surface.
2
In the Settings window for Surface, click Replace Expression in the upper-right corner of the Expression section. From the menu, choose Component 1 (comp1)>Pressure Acoustics, Frequency Domain>Material properties>acpr.c - Speed of sound - m/s.
3
Locate the Coloring and Style section. From the Color table transformation list, choose Reverse.
4
In the Wave Speed toolbar, click  Plot.
Surface 2
1
In the Model Builder window, right-click Wave Speed and choose Surface.
2
In the Settings window for Surface, locate the Expression section.
3
In the Expression text field, type pelw.cp_fast.
4
Click to expand the Title section. From the Title type list, choose None.
5
Locate the Inherit Style section. From the Plot list, choose Surface 1.
Filter 1
1
Right-click Surface 2 and choose Filter.
2
In the Settings window for Filter, locate the Element Selection section.
3
In the Logical expression for inclusion text field, type (y>-H/3).
4
In the Wave Speed toolbar, click  Plot.
Surface 3
1
In the Model Builder window, under Results>Wave Speed right-click Surface 2 and choose Duplicate.
2
In the Model Builder window, click Surface 3.
3
In the Settings window for Surface, locate the Expression section.
4
In the Expression text field, type pelw.cp_slow.
Filter 1
1
In the Model Builder window, click Filter 1.
2
In the Settings window for Filter, locate the Element Selection section.
3
In the Logical expression for inclusion text field, type (y<-H/3)*(y>-2*H/3).
Surface 4
1
In the Model Builder window, under Results>Wave Speed right-click Surface 2 and choose Duplicate.
2
In the Model Builder window, click Surface 4.
3
In the Settings window for Surface, locate the Expression section.
4
In the Expression text field, type pelw.cs_poro.
Filter 1
1
In the Model Builder window, click Filter 1.
2
In the Settings window for Filter, locate the Element Selection section.
3
In the Logical expression for inclusion text field, type (y<-2*H/3).
Annotation 1
1
In the Model Builder window, right-click Wave Speed and choose Annotation.
2
In the Settings window for Annotation, locate the Annotation section.
3
In the Text text field, type Speed of sound = eval(acpr.c) m/s.
4
Locate the Position section. In the X text field, type W.
5
In the Y text field, type H/2.
Annotation 2
1
Right-click Wave Speed and choose Annotation.
2
In the Settings window for Annotation, locate the Annotation section.
3
In the Text text field, type Speed of fast pressure wave = eval(real(pelw.cp_fast)) m/s.
4
Locate the Position section. In the X text field, type W.
5
In the Y text field, type -H/6.
Annotation 3
1
Right-click Wave Speed and choose Annotation.
2
In the Settings window for Annotation, locate the Annotation section.
3
In the Text text field, type Speed of slow pressure wave = eval(real(pelw.cp_slow)) m/s.
4
Locate the Position section. In the X text field, type W.
5
In the Y text field, type -3*H/6.
Annotation 4
1
Right-click Wave Speed and choose Annotation.
2
In the Settings window for Annotation, locate the Annotation section.
3
In the Text text field, type Speed of shear-wave = eval(real(pelw.cs_poro)) m/s.
4
Locate the Position section. In the X text field, type W.
5
In the Y text field, type -5*H/6.
6
In the Wave Speed toolbar, click  Plot.
Create Figure 2 from Ref. 1 by plotting the absolute value of the reflection coefficient as function of the angle of incidence for four values of the frequency f0.
Reflection Coefficient vs. Angle of Incidence
1
In the Home toolbar, click  Add Plot Group and choose 1D Plot Group.
2
In the Settings window for 1D Plot Group, type Reflection Coefficient vs. Angle of Incidence in the Label text field.
3
Locate the Data section. From the Dataset list, choose Study 1 - Detailed Angular Sweep/Parametric Solutions 1 (sol2).
4
Click to expand the Title section. From the Title type list, choose Label.
5
Locate the Plot Settings section.
6
Select the y-axis label check box. In the associated text field, type Reflection coefficient: |R|.
7
Locate the Axis section. Select the Manual axis limits check box.
8
In the x minimum text field, type -2.
9
In the x maximum text field, type 92.
10
In the y minimum text field, type -0.05.
11
In the y maximum text field, type 1.05.
12
Locate the Legend section. From the Position list, choose Upper left.
Global 1
1
Right-click Reflection Coefficient vs. Angle of Incidence and choose Global.
2
In the Settings window for Global, locate the y-Axis Data section.
3
In the table, enter the following settings:
Expression
Unit
Description
absR
1
Absolute and averaged value of the reflection coefficient
4
Locate the x-Axis Data section. From the Axis source data list, choose Inner solutions.
5
From the Parameter list, choose Expression.
6
In the Expression text field, type theta0.
7
From the Unit list, choose °.
8
Click to expand the Legends section. Find the Include subsection. Clear the Description check box.
9
In the Reflection Coefficient vs. Angle of Incidence toolbar, click  Plot.
The plot should look like the one in Figure 2.
Next, create a second study to sweep over the frequency for a given number of angles of incidence.
Add Study
1
In the Home toolbar, click  Add Study to open the Add Study window.
2
Go to the Add Study window.
3
Find the Studies subsection. In the Select Study tree, select General Studies>Frequency Domain.
4
Click Add Study in the window toolbar.
5
In the Home toolbar, click  Add Study to close the Add Study window.
Study 2 - Detailed Frequency Sweep
1
In the Model Builder window, click Study 2.
2
In the Settings window for Study, type Study 2 - Detailed Frequency Sweep in the Label text field.
3
Locate the Study Settings section. Clear the Generate default plots check box.
Step 1: Frequency Domain
1
In the Model Builder window, under Study 2 - Detailed Frequency Sweep click Step 1: Frequency Domain.
2
In the Settings window for Frequency Domain, locate the Study Settings section.
3
In the Frequencies text field, type f0.
4
Locate the Study Extensions section. Select the Auxiliary sweep check box.
5
From the Sweep type list, choose All combinations.
6
Click  Add.
7
In the table, enter the following settings:
Parameter name
Parameter value list
Parameter unit
theta0 (Angle of incidence)
30[deg] 40[deg] 50[deg] 60[deg] 70[deg] 80[deg]
deg
Parametric Sweep
1
In the Study toolbar, click  Parametric Sweep.
2
In the Settings window for Parametric Sweep, locate the Study Settings section.
3
Click  Add.
4
In the table, enter the following settings:
Parameter name
Parameter value list
Parameter unit
f0 (Driving frequency)
10^{range(1,4/14,5)}
Hz
5
In the Study toolbar, click  Compute.
Results
Create Figure 3 from Ref. 1 by plotting the absolute value of the reflection coefficient as function of the frequency f0 for six values of the angle of incidence.
Reflection Coefficient vs. Driving Frequency
1
In the Home toolbar, click  Add Plot Group and choose 1D Plot Group.
2
In the Settings window for 1D Plot Group, type Reflection Coefficient vs. Driving Frequency in the Label text field.
3
Locate the Data section. From the Dataset list, choose None.
4
Locate the Title section. From the Title type list, choose Label.
5
Locate the Axis section. Select the Manual axis limits check box.
6
In the x minimum text field, type 8.5.
7
In the x maximum text field, type 117171.
8
In the y minimum text field, type -0.05.
9
In the y maximum text field, type 1.05.
10
Select the x-axis log scale check box.
11
In the Reflection Coefficient vs. Driving Frequency toolbar, click  Plot.
12
Locate the Legend section. From the Position list, choose Lower left.
Global 1
1
Right-click Reflection Coefficient vs. Driving Frequency and choose Global.
2
In the Settings window for Global, locate the Data section.
3
From the Dataset list, choose Study 2 - Detailed Frequency Sweep/Parametric Solutions 2 (sol8).
4
From the Parameter selection (freq, theta0) list, choose Manual.
5
In the Parameter indices (1-6) text field, type 1.
6
Locate the y-Axis Data section. In the table, enter the following settings:
Expression
Unit
Description
absR
1
Absolute and averaged value of the reflection coefficient
7
Locate the x-Axis Data section. From the Axis source data list, choose Outer solutions.
8
From the Parameter list, choose Expression.
9
In the Expression text field, type f0.
10
Locate the Legends section. From the Legends list, choose Manual.
11
In the table, enter the following settings:
Legends
theta0 = 30[deg]
Global 2
1
Right-click Global 1 and choose Duplicate.
2
In the Settings window for Global, locate the Data section.
3
In the Parameter indices (1-6) text field, type 2.
4
Locate the Legends section. In the table, enter the following settings:
Legends
theta0 = 40[deg]
Global 3
1
Right-click Global 2 and choose Duplicate.
2
In the Settings window for Global, locate the Data section.
3
In the Parameter indices (1-6) text field, type 3.
4
Locate the Legends section. In the table, enter the following settings:
Legends
theta0 = 50[deg]
Global 4
1
Right-click Global 3 and choose Duplicate.
2
In the Settings window for Global, locate the Data section.
3
In the Parameter indices (1-6) text field, type 4.
4
Locate the Legends section. In the table, enter the following settings:
Legends
theta0 = 60[deg]
Global 5
1
Right-click Global 4 and choose Duplicate.
2
In the Settings window for Global, locate the Data section.
3
In the Parameter indices (1-6) text field, type 5.
4
Locate the Legends section. In the table, enter the following settings:
Legends
theta0 = 70[deg]
Global 6
1
Right-click Global 5 and choose Duplicate.
2
In the Settings window for Global, locate the Data section.
3
In the Parameter indices (1-6) text field, type 6.
4
Locate the Legends section. In the table, enter the following settings:
Legends
theta0 = 80[deg]
5
In the Reflection Coefficient vs. Driving Frequency toolbar, click  Plot.
The plot should look like the one in Figure 3.
Finally, plot the absorption coefficient as function of the frequency f0 for six values of the angle of incidence. Simply duplicate the plot with the absorption coefficient, change absR to alpha in the plots, and rename titles where necessary.
Absorption Coefficient vs. Angle of Incidence
1
In the Model Builder window, right-click Reflection Coefficient vs. Angle of Incidence and choose Duplicate.
2
In the Settings window for 1D Plot Group, type Absorption Coefficient vs. Angle of Incidence in the Label text field.
3
Locate the Plot Settings section. In the y-axis label text field, type Absorption coefficient: \alpha.
4
Locate the Legend section. From the Position list, choose Lower left.
Global 1
1
In the Model Builder window, expand the Absorption Coefficient vs. Angle of Incidence node, then click Global 1.
2
In the Settings window for Global, locate the y-Axis Data section.
3
In the table, enter the following settings:
Expression
Unit
Description
alpha
1
Absorption coefficient
4
In the Absorption Coefficient vs. Angle of Incidence toolbar, click  Plot.
The plot should look like the one in Figure 4.