Stability condition of finite difference solution for viscoelastic wave equations

1.   Introduction

The finite difference (FD) method plays one of the most important numerical roles for solving wave equations. Stability and numerical dispersion are two inherent and significant problems in finite-different numerical solution. Discussion of the two problems usually aims at elastic medium. As a result, the obtained stability condition and the analytic result of dispersion are also based on elastic assumption. However, the actual Earth is imperfectly elastic medium and its viscoelastic properties will make the stability of finite difference solution worse. When we solve viscoelastic wave equations with FD method, the former stability condition for elastic equations does not guarantee the stability of the solution.

On the research of stability problems, Aki and Richards (1980) gave examples of finite-difference stability algorithm, then the stability condition expressions of the first-order velocity-stress equation for P-wave and SV-wave were given respectively by Virieux (1986). Some researchers, such as Levander (1988), Carcione et al (2002), Moczo et al (2000), de Basabe and Sen (2007) also analyzed the stability of finite difference scheme of elastic wave equations successively. Manning and Margrave (2000) introduced a correction method for instability problem in 2-D elastic finite-difference forward modeling. Long-time stable modeling of elastic wave with free surface was discussed by Hesthoim (2003). The effects of boundary condition on the instability of finite difference solution in elastic medium have been elaborated by Ilan and Loewenthal (1976), while Vidale and Clayton (1986) discussed the stability of finite difference solutions in elastic medium with the free surface boundary. In China, Dong et al (2000) progressed a quantitative analysis on stability of staggered grid finite difference in elastic medium. And Mu and Pei (2005) presented the stability condition of high-order difference scheme of elastic wave equation. For viscoelastic case, Robertsson et al (1994) explored the FD method of viscoelastic medium, and they introduced the highest allowable formation velocity in stable difference scheme under ultra high frequency, but that for low frequency case which involves in exploration seismology and earthquake research is absent. Saenger and Bohlen (2004) explored the rotated staggered-grid finite difference forward modeling methods of viscoelastic anisotropic medium, and they described the stability condition of difference scheme. Rotated staggered-grid difference scheme is superior to conventional rectangular grid in suppressing numerical dispersion, but its stability is weaker than that of conventional grid to some extent.

It is evident that the status quo on stability study mostly focuses on the elastic medium. In the finite difference solution of viscoelastic wave equation, the visco-elasticity would increase the numerical solution's instability. Based on the viscoelastic theory, the stability conditions of the popularly used Kelvin-Voigt model and Maxwell model in wave field modeling with conventional rectangular grid are explored and further relationships between stability conditions and seismic wave velocity, frequency, quality factor, grid size and difference coefficients are investigated.

2.   Derivation of the stability conditions of basic viscoelastic models

Studies of anelastic medium mostly involve two basic models, namely Kelvin-Voigt model and Maxwell model. The former suggests that medium stress under certain strain consists of two parts, strain related and strain rate related. And it can be regarded as a parallel connection of an elastic element (spring) and a Newtonian viscous element (damper). While the latter indicates that medium strain under certain deformation is composed of two parts, namely, stress related and stress rate related. It can be seen as a series connection of an elastic element (spring) and a Newtonian viscous element (damper), it is also known as the elastic viscous body. Both of the stress-strain constitutive relational expressions are respectively given by:

for Kelvin-Voigt model and

for Maxwell model, where σ is the stress, ε is the strain and M₁, M₂ are elastic coefficient and viscous coefficient of the medium, respectively.

With the harmonic stress σ=Pe^iωt and strain ε=Ee^jωt, the complex modulus of the model can be expressed as M=P/E=ke^jδ, the quality factors Q are satisfied, respectively

for Kelvin-Voigt model, and

for Maxwell model, where δ is the argument of complex modulus.

It can be drawn that the former Q is inversely proportional to frequency ω, and the latter one is proportional to frequency ω. As the relationships of Q versus ω for the two models are different, the impacts of quality factor on the stability also vary. According to constitutive relations, the stability conditions of finite difference solutions for the two models are derived in the following.

2.1   Kelvin-Voigt viscoelastic model

In 3D medium, the bulk strain of elastic body is

where u is the displacement vector. Substitute equation (3) into equation (1a), suppose stress is P, we can get

The wave equation can be rewritten as:

We simultaneously calculate second-order partial derivatives with respect to t in equation (4), and first-order partial derivatives to x, y, z in (5) as following:

Then we do space Fourier transform to stress P of equation (6). The expression in wave number domain is

where is the space Fourier transform of P and k is the wave number.

Do difference approximation to the time derivative in equation (7), we get

where at time n-1, n, and n+1. Assuming , then

Also, we have relation as

(Sun and Yin, 2007), where v is wave velocity in this medium. From equation (2a), we obtain M₂=M₁/(Qω). Thus

Obviously, the stability condition of recursive equation (9) is the modulus of eigenvalues λ of state transition matrix A must be less than 1, where

Suppose c=a+b, then the eigenvalues λ of the matrix A is

And the stability condition of Kelvin-Voigt model becomes

where Δt is time sampling interval and v′²=[(Q²+1)^1/2+Q]Qv²/[2(Q²+1)]. In the case of 2N-order spatial difference accuracy, assuming that the grid spacing in x, y, z directions are Δx, Δy, Δz and wave-number components are k_x, k_y, k_z, respectively. It is simple to obtain (Mu and Pei, 2005)

where a_l is the difference coefficient of 2N-order accuracy, specific values are shown in Table 1 and they satisfy a_-l=a_l. As a result,

Table  1.  Weighting coefficients of the second derivative with 2N-order accuracy (from Sun, 2007)

 | Show Table

DownLoad: CSV

By substituting equation (12) into equation (11), the expression of the stability condition becomes

Specially, when Δx=Δy=Δz, stability condition is simplified as

Generally, difference scheme is often carried out in 2D case, and the stability condition can be further simplified as

When M₂→0, Q→∞, the above expression deteriorates to , which coincides well with the stability condition in elastic medium.

2.2   Maxwell viscoelastic model

For the Maxwell viscoelastic model, the derivation is similar to Kelvin-Voigt case. Main difference lies in that the Maxwell model is originated from equation (2b), then M₂=M₁Q/ω, that is M₁/M₂=ω/Q, also we have approximation as M₁/ρ.=[(Q²+1)^1/2+Q]·v²/(2Q). After a comparable derivation, the stability conditions are as follows:

where Δt is time sampling interval, and v′²=[(Q²+1)^1/2+Q]·v²/(2Q).

In the allowable range of the wave number, stability expression is given by

Specially, when Δx=Δy=Δz=h, it can be simplified as

For the 2D case, there is

Similarly, when M₂→∞, Q₂→∞, the above equation can be rewritten as , which is the perfectly elastic case.

2.3   Low quality factor effect

Conclusions in equations (13)-(15) and (17)-(18) are derived under the assumption of Q≥5 which is described after equation (9). When quality factor Q is smaller than five, M₁/ρ ≈ v² becomes unacceptable, the above conditions can not guarantee the stability of numerical result. Actually, for arbitrary quality factor Q, the wave velocities satisfy

for Kelvin-Voigt model, and

for Maxwell model, respectively (Sun and Yin, 2007), where v²=M₁/ρ. In such case, it needs to replace v_K-Ⅴ for the velocity v in the conclusions (13)-(15), and substitute v_M for the velocity v in (17)-(18). However, quality factors of most rocks and layers can satisfy the approximate condition of Q≥5, so in actual calculation they can still maintain the stability with above approximate condition. The following numerical examples are all carried out using approximated condition; the results give a proof to this conclusion.

3.   Numerical analysis of stability of basic viscoelastic models

Assume that Δx=Δy=Δz=5 m, the analytic stability conditions for 3D Kelvin-Voigt viscoelastic wave equation with sixth-order spatial accuracy are shown in Figure 1. Figure 1a shows the variation of requested time sampling intervals with the formation velocity in the case of Q=5, 50, 100 and 1 000. For a certain value of Q, as velocity increases, the time sampling interval becomes smaller and smaller, that is to say the stability worsens; when the velocity is constant, with the decrease of the media quality factor Q, the stability worsen, too; and with the increase of the quality factor, the more tendency to elastic case, the better the stability. Difference between the two curves that Q=100 and Q=1 000 is very small, while the interval widens between the two curves of Q=50 and Q=100, and the intervals between the curve of Q=5 and other curves are much remarkable. It is demonstrable that as the media quality factor decreases, the stability condition requested by difference solution will sharply be tougher. Figure 1b shows the relationships between the requested time sampling interval and quality factors under different velocities for a stable scheme. From this figure, the same conclusions as Figure 1a can be drawn. Figure 2 shows the comparison of stabilities under different order spatial difference accuracy. For the same velocity model, the higher precision the spatial difference has, the smaller time sampling interval is needed; that is to say, higher spatial difference precision means worse stability. Under such circumstance, high-density time sampling interval is necessary to keep stability. The conventional method usually adopted is improving difference accuracy to suppress the difference spatial dispersion, but this way of work also elicits higher requirements on time stability, the two are a contradiction.

Figure  1.  The relationship between stability, quality factor and velocity of Kelvin-Voigt model where velocity is in unit of km/s.

DownLoad: Full-Size Img PowerPoint

Figure  2.  The relationship between stability and spatial accuracy when Q equals to 100.

DownLoad: Full-Size Img PowerPoint

To be convenient for model trail, 2D case is considered. Grid size and spatial difference accuracy are the same as above mentioned. With Assumption of wave velocity is v=3 km/s and frequency is 35 Hz, the calculated stability condition are Δt≤0.422 ms, Δt≤0.914 4 ms and Δt≤0.958 7 ms when Q equals to 5, 100 and ∞, respectively. We conduct the finite difference numerical test according to the results, and the trends with the number of iterations are shown in Figure 3, where abscissa is the cycle number n of iterations and vertical coordinate is the numerical results. For Q =5 (Figure 3a), it is stable when Δt=0.3 ms and Δt=0.4 ms, but unstable as Δt=0.45 ms. For Q=100 (Figure 3b), it is stable when Δt=0.5 ms and Δt=0.9 ms, but unstable when Δt= 0.92 ms. All of the above verify the validity of the results.

Figure  3.  The convergence of Kelvin-Voigt viscoelastic model by finite difference for Q=5 (a) and Q=100 (b) where n denotes the cycle number of calculation.

DownLoad: Full-Size Img PowerPoint

Figure 4 shows the 3-D analytic solution of stability for Maxwell model with sixth-order spatial differential accuracy. The basic characteristics are same as that of Kelvin-Voigt model. With the increase of velocity and the improvement of spatial difference accuracy, the stability tends to be weaker and the time sampling interval should be smaller. However, the influence of quality factor on stability is less significant than that in the former model. Considering a sixth-order scheme, supposing the ratio of grid size to velocity is 0.003, when quality factor changes from 5 to 1 000, the changes of requested time sampling interval is only less than 2%. The influence is so small that can be neglected. Figure 5 shows the relationships between stability of two models and frequency of seismic wave. On condition that the other parameters are same, with the increase of frequency, the stability of difference solution of Kelvin-Voigt model gets better, while that of Maxwell model becomes worse.

Figure  4.  The variation of stability with quality factor, spatial accuracy and velocity (in unit of km/s) of Maxwell model.

DownLoad: Full-Size Img PowerPoint

Figure  5.  The relationship between stability and frequency for Kelvin-Voigt model (a) and Maxwell model (b).

DownLoad: Full-Size Img PowerPoint

The above numerical results are all based on the approximate conditions (13)-(15) and (17)-(18). From equation (19) we can see that the real velocity of KelvinVoigt model is greater than approximate velocity, while that of Maxwell model is smaller. To illustrate the error scale caused by the velocity of approximation, Figure 6 shows the changes of relative errors between real velocities and their approximations with quality factor for the two viscoelastic models. With the increase of quality factor, the relative errors tend to be zero, and the velocities of viscoelastic models approach to the velocity of elastic media. The relative errors of the two models are less than 1.48% and 0.49% when Q≥5, respectively. When Q is too small, equations (13)-(15) can not guarantee the stability of difference scheme because the real velocity of Kelvin-Voigt model is much larger than the approximate value, while the equations (17)-(18) can guarantee the stability due to a smaller velocity of Maxwell model than the approximate value.

Figure  6.  The relative errors between accurate velocities and their approximation of two viscoelastic models.

DownLoad: Full-Size Img PowerPoint

4.   Conclusions

From some theory analysis, stability conditions of Kelvin-Voigt model and Maxwell model with varied difference accuracy are analyzed. Combined with analysis on numerical tests, the following conclusions can be drawn.

1) The greater the quality factor Q and the lower the velocity, the better the stability of difference scheme of the two models. Quality factor has less influence on stability of Maxwell model than that of Kelvin-Voigt model.

2) On condition that value of Q is fixed, with the increase of the frequency, the stability of difference scheme of Kelvin-Voigt model gets better, while that of Maxwell model becomes worse with higher frequency.

3) The stability conditions obtained from this paper are suitable for quality factor Q≥5; when Q is too small, the stability of difference scheme of Kelvin-Voigt model deteriorates dramatically, for which accurate velocity are needed to define the stability conditions.

From the analysis we can see that the stability of difference scheme of Maxwell model is superior to that of Kelvin-Voigt model to some extent. But Kelvin-Voigt model is popularly used in solving viscoelastic problems at present because there are certain gaps between inherent Q~ω relationship of Maxwell model and the current observational results.