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 incident wave and wave vector are defined as variables in the model.
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. 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 to the incident waves (at the interface) and is given by
where the subscript sc stands for scattered and tot for total. The absorption coefficient gives a measure of the relative absorbed energy and is given by
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, 72°. 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 72° 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 72° 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 model parameters from a file. The list contains definitions of material parameters (see Table 1), geometry lengths, and the two model parameters (frequency f0 and angle of incidence theta0).
Global Definitions
Parameters 1
1
In the Model Builder window, under Global Definitions click Parameters 1.
2
In the Settings window for Parameters, locate the Parameters section.
3
Click  Load from File.
4
Browse to the model’s Application Libraries folder and double-click the file reflections_water_sediment_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 two nonlocal couplings. The two couplings act as probes picking up variable values. One samples a point in the water domain and the other evaluates the average 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.
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
5
In the Model Builder window, click Pressure Acoustics, Frequency Domain (acpr).
6
In the Settings window for Pressure Acoustics, Frequency Domain, locate the Typical Wave Speed for Perfectly Matched Layers section.
7
In the cref text field, type real(acpr.c_c).
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. From the Pressure field type list, choose User defined.
5
In the pb text field, type Pin.
Materials
Material 1 (mat1)
In the Model Builder window, under Component 1 (comp1) right-click Materials and choose Blank Material, twice, to create two new blank materials.
Sediment
1
In the Model Builder window, under Component 1 (comp1)>Materials click Material 2 (mat2).
2
In the Settings window for Material, type Water in the Label text field.
3
In the Label text field, type Sediment.
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
Material 1 (mat1)
1
In the Model Builder window, under Component 1 (comp1)>Materials click Material 1 (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 order to reduce the computation cost of the analysis, use a mapped mesh.
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.
Specify the wavelength of 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.
Specify the wavelength of the slowest wave traveling through 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. Select the Maximum element size check box.
6
In the associated text field, type cs_poro/f0/5.
7
Select the Minimum element size check box.
8
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.
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],90[deg])
rad
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
Clear the Parameter indicator text field.
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. From the Color table list, choose Wave.
5
Select the Symmetrize color range check box.
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 (rad)) list, choose 4: freq=10000 Hz, theta0=0.5236 rad.
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 (rad)) list, choose 4: freq=10000 Hz, theta0=0.5236 rad.
5
Locate the Title section. From the Title type list, choose Manual.
6
In the Title text area, type Reflected Pressure.
7
Clear the Parameter indicator text field.
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. From the Color table list, choose Wave.
5
Select the Symmetrize color range check box.
6
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 (rad)) list, choose 4: freq=10000 Hz, theta0=0.5236 rad.
5
Locate the Title section. From the Title type list, choose Manual.
6
In the Title text area, type Total Pressure.
7
Clear the Parameter indicator text field.
Surface 1
1
Right-click Total Pressure and choose Surface.
2
In the Settings window for Surface, locate the Expression section.
3
In the Expression text field, type if(dom==3,acpr.p_t,pelw.p_t).
The if() statement plots acpr.p_t if in the water domain (domain 3) else it plots pelw.p_t (the sediment domain).
4
Locate the Coloring and Style section. From the Color table list, choose Wave.
5
Select the Symmetrize color range check box.
6
In the Total Pressure toolbar, click  Plot.
7
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
Clear the Parameter indicator text field.
6
Locate the Data section. From the Dataset list, choose Array 2D 1.
7
From the Parameter value (freq (Hz),theta0 (rad)) list, choose 4: freq=10000 Hz, theta0=0.5236 rad.
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. Select the Scale factor check box.
6
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, locate the Expression section.
3
From the menu, choose Component 1 (comp1)>Pressure Acoustics, Frequency Domain>Material properties>acpr.c - Speed of sound - m/s.
4
In the Expression text field, type acpr.c.
5
Locate the Coloring and Style section. Select the Reverse color table check box.
6
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
Click to expand 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
Select the Allow evaluation of expressions check box.
5
Locate the Position section. In the X text field, type W.
6
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
Select the Allow evaluation of expressions check box.
5
Locate the Position section. In the X text field, type W.
6
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
Select the Allow evaluation of expressions check box.
5
Locate the Position section. In the X text field, type W.
6
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
Select the Allow evaluation of expressions check box.
5
Locate the Position section. In the X text field, type W.
6
In the Y text field, type -5*H/6.
7
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 None.
4
Click to expand 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 -2.
7
In the x maximum text field, type 92.
8
In the y minimum text field, type -0.05.
9
In the y maximum text field, type 1.05.
10
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 Data section.
3
From the Dataset list, choose Study 1 - Detailed Angular Sweep/Parametric Solutions 1 (sol2).
4
From the Parameter selection (f0) list, choose From list.
5
In the Parameter values list, select 10.
6
Locate the y-Axis Data section. In the table, enter the following settings:
Expression
Unit
Description
absR
1
Absolute value of the reflection coefficient
7
Locate the x-Axis Data section. From the Parameter list, choose Expression.
8
In the Expression text field, type theta0.
9
From the Unit list, choose °.
10
Click to expand the Legends section. From the Legends list, choose Manual.
11
In the table, enter the following settings:
Legends
f0 = 10 Hz
12
In the Reflection Coefficient vs. Angle of Incidence toolbar, click  Plot.
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 values list, select 100.
4
Locate the Legends section. In the table, enter the following settings:
Legends
f0 = 100 Hz
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 values list, select 1000.
4
Locate the Legends section. In the table, enter the following settings:
Legends
f0 = 1 kHz
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 values list, select 10000.
4
Locate the Legends section. In the table, enter the following settings:
Legends
f0 = 10 kHz
5
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]
rad
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
In the Model Builder window, click Reflection Coefficient vs. Driving Frequency.
13
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 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 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
Global 2
1
In the Model Builder window, click Global 2.
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
Global 3
1
In the Model Builder window, click Global 3.
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
Global 4
1
In the Model Builder window, click Global 4.
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.