Isotropic-Anisotropic Sample:Elastic Wave Propagation

Isotropic-Anisotropic Sample:
Elastic Wave Propagation

This tutorial shows how to use the Elastic Waves, Time Explicit physics interface to model the propagation of elastic waves in heterogeneous anisotropic linear elastic media. The case where the material properties differ from one subdomain to another requires an special handling. Namely, the Material Discontinuity interior boundary condition (or the Continuity pair condition for assemblies) is to be explicitly imposed upon the interface between such subdomains.

In addition, the model demonstrates how to retrieve the displacements at the points of interest from the solution that contains the structural velocity and the strain.

Model Definition

A two-dimensional space is occupied by two linear elastic materials: an anisotropic material to the left of the axis x = 0 and an isotropic material to the right. The material properties are given in Table 1 (kg/m3 for the density and 1010 Pa for the stiffness). The other material parameters are c33 = c22, c13 = c23 = c12, and c44 = c55 = c66.

Table 1: Material properties.
	ρ	c11	c12	c22	c66
isotropic	7100	16.5	8.85	16.5	3.96
anisotropic	7100	16.5	5.00	6.2	3.96

The elastic waves are excited by a point force applied in the y direction at the point
x0 = −2 cm, y0 = 0 cm. The source distribution in time is given by a Ricker (also known as Mexican hat) wavelet depicted in Figure 1. It has the dominant frequency f0 = 170 kHz, the time delay t0 = 6 μs, and the amplitude of 1013 N.

The elastic 2D space is modeled as a square with the side length of 2L, L = 20 cm and the center located at (0, 0). The subdomains with the anisotropic and the isotropic properties lie to the left and to the right of the line x = 0, respectively. Four listening points are put to record the seismograms (displacement at points) at the coordinates (−10.5 cm, −8 cm), (−3.5 cm, −8 cm), (−1.0 cm, −8 cm), and (10.5 cm, −8 cm). The model geometry layout is shown in Figure 2.

The propagation of the elastic waves in infinite space is modeled by imposing the Absorbing Layers from all four sides of the square. Combined with the Low-Reflecting Boundary condition that is imposed upon the outer boundary of the square, the Absorbing Layers ensure the reduction of the spuriously reflected waves in the physical domain.

Figure 1: A Ricker wavelet that defines the source distribution in time.

The Elastic Waves, Time Explicit interface solves for the structural velocity and the strain. The displacements can be retrieved from the velocities by solving the corresponding ordinary differential equation

In this model, this ODE is only solved at the probe points.

Figure 2: Model geometry layout.

Results and Discussion

Figure 3 shows the quantitative evolution of the waves propagating from the source at t = 30, 60, and 90 μs. One can see the round profiles of the longitudinal and the shear waves in the isotropic material on the right side. At the same time, the quasilongitudinal and the quasishear waves in the anisotropic material on the right side have the elliptic and the conical profiles, respectively. Note that the absorbing layer domains are not shown in Figure 3.

Figure 4 shows the profiles of the vertical displacement (seismograms) retrieved from the computed velocities at the probe points. The results have a good agreement in their shapes with those given in Ref. 1.

Figure 3: Velocity magnitude profiles at 30, 60, and 90 μs.

Figure 4: Seismograms at the probe points.

Notes About the COMSOL Implementation

The force load at the point (x0, y0) is modeled as a Body Load with a Total force applied to the physical domain. In theory, the load applied on the domain level must have proportional to the Dirac delta to be equivalent to a point load. That is,

One of the Dirac delta representations reads

This fact is used to approximate the point source by a domain source that has the form of a product of the temporal part and the spatial part. Here, the former is the Ricker wavelet depicted in Figure 1 and the latter is a Gaussian pulse with a relatively small extent a.

The default approximation order used in the Elastic Wave, Time Explicit interface is quartic. In this case, the maximum mesh element size that resolves the minimum wavelength λmin must not exceed λmin/1.5. The minimum wavelength is defined by the frequency content of the source, which is the Ricker wavelet here. Its frequency domain representation reads

where f0 is the dominant frequency. Although the upper cutoff frequency of the wavelet is 3f0, λmin is chosen to correspond to 2f0 in this model. This is a compromise that on one hand accounts for the greatest part of the source energy, and on the other hand does not result in a too fine mesh that would slow down the simulation and tremendously increase the model file size.

Reference

1. J. de la Puente, M. Käser, M. Dumbser and H. Igel, “An arbitrary high-order discontinuous Galerkin method for elastic waves on unstructured meshes – IV. Anisotropy”, Geophys. J. Int., vol. 169, issue 3, 2007.

Acoustics_Module/Elastic_Waves/isotropic_anisotropic_sample

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select .

Click .

Click

In the tree, select .

Click

Global Definitions

Parameters 1

In the window, under click .

In the window for , locate the section.

Click

Browse to the model’s Application Libraries folder and double-click the file isotropic_anisotropic_sample_parameters.txt.

Create the source time function given by a Ricker wavelet.

Definitions

Analytic 1 (an1)

In the toolbar, click

and choose .

In the window for , type A in the text field.

Locate the section. In the text field, type (1 - 2*pi^2*f0^2*(t - t0)^2)*exp(-pi^2*f0^2*(t - t0)^2).

In the text field, type t.

Locate the section. In the text field, type s.

In the text field, type N.

Locate the section. In the table, enter the following settings:


Argument	Lower limit	Upper limit
t	0	5*t0

Click

The signal should look like the one in Figure 1.

Results

Impulse Frequency Content

In the window for , type Impulse Frequency Content in the text field.

Click to expand the section. In the text area, type FFT of G (Pa).

Locate the section. Select the check box.

In the associated text field, type FFT of the signal (Pa).

Function 1

In the window, expand the node, then click .

In the window for , locate the section.

From the list, choose .

In the toolbar, click

Select the check box.

In the text field, type 1000000.

In the toolbar, click

The Fourier transformation of the signal should look like this. The plot indicates that the mesh should resolve frequency content up to 400 to 500 kHz.

Now, create a two-dimensional Gaussian function that will be used to approximate the Dirac delta at (x0,y0).

Definitions

Analytic 2 (an2)

In the toolbar, click

and choose .

In the window for , type G in the text field.

Locate the section. In the text field, type 1/(pi*dS)*exp(-((x - x0)^2 + (y - y0)^2)/dS).

In the text field, type x, y.

Locate the section. In the text field, type m.

In the text field, type 1.

Locate the section. In the table, enter the following settings:


Argument	Lower limit	Upper limit
x	-0.04	0
y	-0.02	0.02

Click

Geometry 1

Square 1 (sq1)

In the toolbar, click

In the window for , locate the section.

In the text field, type 2*L.

Locate the section. From the list, choose .

Click to expand the section. Select the check box.

In the table, enter the following settings:


Layer name	Thickness (m)
Layer 1	L/6

Select the check box.

Create a line segment that separates the subdomains with isotropic and anisotropic materials.

Line Segment 1 (ls1)

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

In the text field, type -L.

Locate the section. From the list, choose .

In the text field, type L.

Create four probe points for calculating the seismograms.

Point 1 (pt1)

In the toolbar, click

In the window for , locate the section.

In the text field, type -10.5[cm] -3.5[cm] -1[cm] 10.5[cm].

In the text field, type -8[cm] -8[cm] -8[cm] -8[cm].

Click

The geometry should look like the one in Figure 2.

Definitions

In the next steps, create Integration operators that will be used to probe the displacement at the probe points.

Integration 1 (intop1)

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Select Point 9 only.

Integration 2 (intop2)

Right-click and choose .

In the window for , locate the section.

Click

Select Point 10 only.

Integration 3 (intop3)

Right-click and choose .

In the window for , locate the section.

Click

Select Point 11 only.

Integration 4 (intop4)

Right-click and choose .

In the window for , locate the section.

Click

Select Point 16 only.

Define a for each of the probe points. These probes will use the integration operators defined previously to measure uy, a variable that will be defined through an ODE that integrates the vertical velocity to obtain the vertical displacement.

Global Variable Probe 1 (var1)

In the toolbar, click

and choose .

In the window for , locate the section.

In the text field, type intop1(uy).

Select the check box.

In the associated text field, type y displacement in (-10.5 cm, -8 cm) (m).

Global Variable Probe 2 (var2)

Right-click and choose .

In the window for , locate the section.

In the text field, type intop2(uy).

In the text field, type y displacement in (-3.5 cm, -8 cm) (m).

Global Variable Probe 3 (var3)

Right-click and choose .

In the window for , locate the section.

In the text field, type intop3(uy).

In the text field, type y displacement in (-1 cm, -8 cm) (m).

Global Variable Probe 4 (var4)

Right-click and choose .

In the window for , locate the section.

In the text field, type intop4(uy).

In the text field, type y displacement in (10.5 cm, -8 cm) (m).

Create groups to facilitate the navigation through the model.

Global Variable Probe 1 (var1), Global Variable Probe 2 (var2), Global Variable Probe 3 (var3), Global Variable Probe 4 (var4)

In the window, under , Ctrl-click to select , , , and .

Right-click and choose .

Global Probes

In the window for , type Global Probes in the text field.

Integration 1 (intop1), Integration 2 (intop2), Integration 3 (intop3), Integration 4 (intop4)

In the window, under , Ctrl-click to select , , , and .

Right-click and choose .

Integration Operators for the Probes

In the window for , type Integration Operators for the Probes in the text field.

Probe Points

In the toolbar, click

In the window for , locate the section.

From the list, choose .

Select Points 9–11 and 16 only.

In the text field, type Probe Points.

Now, set up the absorbing layers (sponge layers) used to truncate the computational domain.

Absorbing Layer 1 (ab1)

In the toolbar, click

Select Domains 1–4, 6, 7, and 9–12 only.

Setting up the materials is simpler if the correct material model is selected in the physics. Here we want to use an anisotropic material.

Elastic Waves, Time Explicit (elte)

Elastic Waves, Time Explicit Model 1

In the window, under click .

In the window for , locate the section.

From the list, choose .

Materials

Material Isotropic

In the window, under right-click and choose .

In the window for , type Material Isotropic in the text field.

Select Domains 7–12 only.

Locate the section. In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Elasticity matrix, Voigt notation	{DVo11, DVo12, DVo22, DVo13, DVo23, DVo33, DVo14, DVo24, DVo34, DVo44, DVo15, DVo25, DVo35, DVo45, DVo55, DVo16, DVo26, DVo36, DVo46, DVo56, DVo66} ; DVoij = DVoji	{16.5e10, 8.58e10, 16.5e10, 8.58e10, 8.58e10, 16.5e10, 0, 0, 0, 3.96e10, 0, 0, 0, 0, 3.96e10, 0, 0, 0, 0, 0, 3.96e10}	Pa	Anisotropic, Voigt notation
Density	rho	7100	kg/m³	Basic

Material Anisotropic

Right-click and choose .

In the window for , type Material Anisotropic in the text field.

Select Domains 1–6 only.

Locate the section. In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Elasticity matrix, Voigt notation	{DVo11, DVo12, DVo22, DVo13, DVo23, DVo33, DVo14, DVo24, DVo34, DVo44, DVo15, DVo25, DVo35, DVo45, DVo55, DVo16, DVo26, DVo36, DVo46, DVo56, DVo66} ; DVoij = DVoji	{16.5e10, 5.0e10, 6.2e10, 5.0e10, 5.0e10, 6.2e10, 0, 0, 0, 3.96e10, 0, 0, 0, 0, 3.96e10, 0, 0, 0, 0, 0, 3.96e10}	Pa	Anisotropic, Voigt notation
Density	rho	7100	kg/m³	Basic