Solution of the dispersion relation

A variety of numerical methods may be used to solve (dispersionrelation), including for example Crout's reduction method (Crout, 1941). However, since the system is relatively small ($3\times 3$) and since our purposes include gaining some physical insight into the processes involved, it will prove instructive to do some more analysis on the problem prior to the ultimate numerical calculations.

First, note that our analysis will be considerably simplified by an obvious rearrangement of the determinant (dispersionrelation) so that

(- v^2) = 0,   where is the identity matrix and

= ^-1^-1^-1,   having complex matrix elements g_ij, for i,j = 1,2,3. The matrix was given previously in (variablechange) and is clearly invertible for finite values of the two porosities. Matrix is also always invertible for realistic choices of material parameters.

Using the properties of determinants, it is not hard to show that the $3\times 3$ determinant (rearrangement) can be reduced to a $2\times 2$ determinant in either of two convenient forms:

&& (g_11-v^2)(g_22-v^2) - g_12g_21 & -1 g_12g_23g_31 + g_13g_32g_21 - (g_22-v^2)g_13g_31 - (g_11-v^2)g_23g_32 & g_33 - v^2 &=& (g_11-v^2)(g_33-v^2) - g_13g_31 & -1 g_12g_23g_31 + g_13g_32g_21 - (g_33-v^2)g_12g_21 - (g_11-v^2)g_23g_32 & g_22 - v^2 &=& 0.   It is useful to write the determinant in each of these two ways in order to make connections with single-porosity models. The first version in (newdet) has as its upper left-hand element an expression that corresponds precisely to the determinant for the $2\times 2$system when the only porosity present is the matrix (or storage) porosity. Similarly, the second version in (newdet) has as its element in the upper left-hand position a term that corresponds to the determinant for the $2\times 2$ system when the only porosity present is the fracture/crack porosity. In each case the remaining terms determine the effect of coupling to the second type of porosity. Depending on other system parameters such as permeabilities and porosities, either of these limits may be useful to consider, and may also provide good starting points for the iterative (Newton-Raphson) procedure that we will use to solve the determinant equation for the complex velocity v.

If we express the determinant (rearrangement) as a polynomial D(x) where x = v², then

D(x) = -x^3 &+& (g_11+g_22+g_33)x^2                  &-& (g_22g_33+g_33g_11+g_11g_22          &-& g_23g_32 - g_13g_31 - g_12g_21)x          &+& (g_11g_22g_33 - g_23g_32g_11 - g_13g_31g_22 &-& g_12g_21g_33 + g_12g_23g_31 + g_13g_32g_21) = 0.   The Newton-Raphson method (Hildebrand, 1956) for solving this equation for x is an iteration process starting from some initial choice x⁽⁰⁾ and computing

x^(i)=x^(i-1)-D(x^(i-1))/D'(x^(i-1)),    for    i = 1,2, ..., N_c   (where D' is the first derivative of D with respect to x) until some convergence criterion has been met at i = N_c. The fact that the coefficients and the solution of the problem are complex adds no special complication to this procedure. However, since the polynomial is complex, it is important that good starting values x⁽⁰⁾ be obtained for at least two of the three roots of (polynomial), as a search procedure in the complex plane would be considerably more difficult to implement than is pure Newton-Raphson iteration.

From (newdet), we see that two choices of starting values for the complex parameter x are given by (for example)

x_1,2^(0) = 12(g_11+g_22 (g_11-g_22)^2 + 4g_12g_21).   Once the Newton-Raphson iteration has converged from both of these starting values to their final values of x₁ and x₂, then the third solution of the dispersion relation is obtained directly by recalling that

D(x) _n=0^3 D_n(-x)^n = -(x-x_1)(x-x_2)(x-x_3).   Using the linear independence of the terms in powers of x in (factoring), we have three equations showing that

x_3 = D_2 - x_1 - x_2 = D_1-x_1 x_2x_1+x_2 = D_0x_1 x_2.   Any one of these three formulas may be used to calculate x₃ directly, or combinations of them (such as the geometric mean of any two) may be used to reduce the errors that might be introduced by premature termination of the Newton-Raphson iteration process for x₁ and x₂. An alternative method (but not the one we used in the examples of the present paper) is to use the formulas of (p3) in a different way to arrive at a quadratic formula for x₂ and x₃, once that x₁ is known from the iteration procedure. The resulting formulas are given (for example) by

x_2,3 = 12[(D_2 - x_1) (D_2 - x_1)^2 - 4D_0/x_1].   In this approach it is clearly preferable to solve for the fast wave solution as x₁, and then use the formulas in (p2p3) to obtain the two slow wave solutions since the fast wave solution will virtually always be well-behaved.

For each value of x that solves the dispersion equation, we can then compute the wave velocity and inverse of the quality factor Q using the definition

1x 1v() 1v_r(1 + i2Q),   where the choice of square root is determined by the requirement of physical realizability as discussed in the following section. This definition of 1/Q is accurate when attenuation losses are low, but should be carefully interpreted for high loss situations (e.g., low frequency diffusive modes). In particular, it is often stated that, when losses are high, (vandQ) needs to be modified so that if

1v() 1v_r(1 + iv_r),   where $\alpha$ is the attenuation coefficient (having units of inverse length), then [as Hamilton (1972) and Bourbié et al. (1987) show]

1Q 2v_r- ^2 v_r^2/.   However, for diffusive modes the imaginary and real parts of the velocity are of comparable size, and therefore it is possible that the formula (highlossQ) will become singular at low frequencies for such modes. This complication can and does happen in practice. To avoid this complication, we use the definition (vandQ) in all cases, and then interpret those situations in which $1/Q \to 2$ as an indication that the mode under consideration is actually diffusive rather than propagatory.