COMSOL 6.3 - Transient Gaussian Explosion

This example introduces some important concepts to have in mind when solving transient problems. In particular, it examines the relationship between the frequency content in the sources driving the model, the mesh resolution, and the time step.

Model Definition

An ellipse with sound-hard walls has the interesting property that an acoustic signal emanating from one of the foci refocuses at the other focal point b/c seconds later, where b (SI unit: m) is the major axis length and c (SI unit: m/s) is the speed of sound.

Inspired by Ref. 1 and Ref. 2, this model involves a Gaussian explosion at one focus of an ellipse to illustrate some properties of time-dependent acoustic problems. The major and minor axis lengths are 10 m and 8 m, respectively. The major axis coincides with the x-axis and the foci are located at x = −3 m and x = 3 m. Because of symmetry, the model can be limited to the upper half plane.

Denoting the fluid density by ρ and the speed of sound by c, the acoustic pressure field p(x, t), inside the elliptical chamber is governed by the scalar wave equation

where the point-source term on the right-hand side is given by

where δ(2) is the two-dimensional Dirac delta function. The time dependence of the explosion is determined by the cutoff Gaussian pulse

describing the rate of airflow (SI unit: m2/s) away from the source, located at x = x0. The parameter f0, which is proportional to the pulse bandwidth, is chosen as f0 = 380 Hz. As the following plots show, by taking τ = 1/f0 the pulse very closely approximates a full Gaussian, the effect of the cutoff tails being numerically insignificant.

Figure 1: Normalized Gaussian pulse (left) and its derivative (right).

A particularly interesting property of the Gaussian function is that its Fourier transform, depicted in Figure 3, is equally simple (neglect cutoff effects):

where ω0 = 2π f0. The magnitude of the Fourier transform falls off quickly for increasing angular frequencies ω. Practically all the energy in the signal is contained in the frequency band −2ω0 < ω < 2ω0 with most of it concentrated between −ω0 and ω0.

Therefore, when using a forcing function of this type, it is enough to resolve wavelengths corresponding to the angular frequency ω0, which in turn corresponds to the frequency f0. In this model, the parameter that defined the maximal frequency to resolve is fmax, it is set equal to 600 Hz which is slightly higher than f0 (and will resolve the pulse fully).

In order to well resolve the propagating pulse in space, the rule of thumb is to use a minimum of 5 to 6 elements per wavelength (at the maximal frequency), when using quadratic Lagrange elements. The typical mesh size is given by

Where hmax is a typical maximal mesh-element size, and N is the number of elements per wavelength required to resolve a harmonic wave with some accuracy. The following discussion uses N = 5. The important point is that there is little to gain in prescribing a forcing function that contains frequencies that the mesh cannot resolve.

In addition to controlling the pulse shape, the mesh resolution imposes a restriction on the internal time steps size taken by the solver. In the transient interfaces of the Acoustics Module, the solver is controlled from the section found at the top physics level. The user should enter the in the model. It is also possible to choose a (method). Here it is recommended to use the default option.

Time Stepping Explanation

The logic for the automatic choice made is as follows. The relationship between mesh size and time-step size is closely related to the CFL number (Ref. 3), which is defined as

This nondimensional number can be interpreted as the fraction of an element the wave travels in a single time step. A CFL number around 1 would correspond to the same resolution in space and time if the discretization errors were of the same size; however, that is normally not the case.

By default, COMSOL Multiphysics uses the implicit second-order accurate method generalized-α to solve transient acoustics problems. In space, the default is second-order Lagrange elements. Generalized-α introduces some numerical damping of high frequencies but much less than the BDF method.

The temporal discretization errors for generalized-α are larger than the spatial discretization errors when second-order elements are used in space. The limiting step size, where the errors are of roughly the same size, can be found somewhere at CFL < 0.2. You can get away with a longer time step if the forcing does not make full use of the mesh resolution; that is, if high frequencies are absent from the outset.

When the excitation contains all the frequencies the mesh can resolve, there is no point in using an automatic time-step control which can be provided by the time-dependent solver. The tolerances in the automatic error control are difficult to tune when there is weak but important high-frequency content. Instead, you can use your knowledge of the typical mesh size, speed of sound, and CFL number to calculate and prescribe a fixed time step. This is exactly the default behavior when the method is chosen in the section. The option corresponds to the automatic time-step control but with some tighter controls of the allowed time-steps. This latter option is still not recommended as the manual option typically yields much better results (and is faster).

The internal time step generated by the option and the entered is set by assuming that the user has generated a mesh that properly resolves the same maximal frequency (minimal wavelength). The following step is generated

Assuming that N is between 5 and 6 and the CFL number is roughly 0.1. These values give a good margin of safety.

Results and Discussion

Using the properties of air for the medium and selecting a mesh density based on the parameters. The model runs for 0.035 s so that you can study the refocusing at the right-hand focus point at roughly 0.0315 s.

Figure 2: The refocusing occurs at roughly 0.0315 s. An animation gives a better feeling for the process.

For the selected combination of mesh size, pulse shape, and time step, the solution can be shown to be both smooth and accurate. Selecting a smaller value for N leads to oscillations if the CFL number is small enough (choosing a too low frequency to resolve), while selecting a higher CFL number (and consequentially a larger time step) leads to an inaccurate solution.

Figure 3 shows the Fourier transform of the source signal. Plotting the frequency content of the source signal is a good way for selecting the maximal frequency to resolve in the model.

Figure 3: Fourier transform of the source signal.

References

1. B. Yue and M.N. Guddati, “Dispersion-reducing Finite Elements for Transient Acoustics”, J. Acoust. Soc. Am., vol. 118, no. 4, pp. 2132–2141, 2005.

2. H.-O. Kreiss, N.A. Peterson, and J. Yström, “Difference Approximations for the second order wave equation,” SIAM J. of Num. Analys., vol. 40, 1940–1967, 2002.

3. R. Courant, K.O. Friedrichs, and H. Lewy, “On the Partial Difference Equations of Mathematical Physics,” IBM Journal, vol. 11, pp. 215–234, 1956.

Acoustics_Module/Tutorials,_Pressure_Acoustics/gaussian_explosion

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

Geometry 1

Ellipse 1 (e1)

In the toolbar, click

In the window for , locate the section.

In the text field, type 5.

In the text field, type 4.

In the text field, type 180.

Point 1 (pt1)

In the toolbar, click

In the window for , locate the section.

In the text field, type -3.

Click

Click the

button in the toolbar.

The completed geometry should look like that in the figure below.

Global Definitions

Parameters 1

In the window, under click .

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Value	Description
c_air	343[m/s]	343 m/s	Speed of sound in air
f_max	600[Hz]	600 Hz	Maximum frequency to resolve
N	5	5	Elements per wavelength
h_max	c_air/f_max/N	0.11433 m	Typical element size
A	4[m^2/s]	4 m²/s	Point-source amplitude
f0	380[Hz]	380 Hz	Source frequency bandwidth
t0	1/f0	0.0026316 s	Source pulse half width

Add Material

In the toolbar, click

to open the window.

Go to the window.

In the tree, select > .

Click the button in the window toolbar.

In the toolbar, click

to close the window.

Pressure Acoustics, Transient (actd)

Enter the maximal frequency to be resolved in the model. This value can be assessed from a Fourier analysis of the sources as depicted in Figure 3.

Under the section enter the and set it to f_max.

Point Source 1

In the toolbar, click

and choose .

Select Point 2 only.

In the window for , locate the section.

From the list, choose .

In the A text field, type A.

In the f0 text field, type f0.

In the tp text field, type t0.

Mesh

Proceed and generate the mesh using the functionality. The frequency controlling the maximum element size is per default taken . Set the desired in the study step. 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, use the default option, which gives 5 elements per wavelength.

Study 1

Step 1: Time Dependent

In the window, under click .

In the window for , locate the section.

In the text field, type range(0,0.5e-3,35e-3).

This setting gives you a solution output at every 0.5 ms from t = 0 to t = 35 ms. This should not be confused with the time steps actually taken by the solver. These are set automatically based on the settings entered above. To see where they end up, do the following:

Solution 1 (sol1)

In the toolbar, click

In the window, expand the node, then click .

In the window for , click to expand the section.

The solver automatically picks up the settings entered at the physics level (in the Transient Solver Settings section) when the default solver is generated or when solve is clicked the first time. It is in general not necessary to edit these settings.

Experienced users maybe find it useful to edit these settings.

Mesh 1

In the window, under right-click and choose .

Study 1

In the toolbar, click

Results

Acoustic Pressure (actd)

The default plot shows the pressure at the final time. To get a more attractive plot, you can add a height.

Height Expression 1

In the window, expand the node.

Right-click and choose .

You can select different times to look at the wave using the parent Acoustic Pressure (actd) node. At t = 0.0315 s, you are close to the moment when the waves refocus.

Acoustic Pressure (actd)

In the window, under click .

In the window for , locate the section.

From the list, choose .

In the toolbar, click

Click the

button in the toolbar.

The resulting plot is found in Figure 2. It is illustrative to animate transient problems in general and wave propagation in particular; you can do this by right-clicking the Export node and adding an Animation feature.

It is also possible to plot the field on all sides of the point source by defining a mirror dataset.

Mirror 2D 1

In the toolbar, click