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 .

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 1e13*(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

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

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/sqrt(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

Explicit 1

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 1 (mat1)

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 2 (mat2)

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