#
Wave impedance matrices for cylindrically anisotropic

radially inhomogeneous elastic solids

###### Abstract

Impedance matrices are obtained for radially inhomogeneous structures using the Stroh-like system of six first order differential equations for the time harmonic displacement-traction 6-vector. Particular attention is paid to the newly identified solid-cylinder impedance matrix appropriate to cylinders with material at , and its limiting value at that point, the solid-cylinder impedance matrix . We show that is a fundamental material property depending only on the elastic moduli and the azimuthal order , that is Hermitian and is negative semi-definite. Explicit solutions for are presented for monoclinic and higher material symmetry, and the special cases of and are treated in detail. Two methods are proposed for finding , one based on the Frobenius series solution and the other using a differential Riccati equation with as initial value. The radiation impedance matrix is defined and shown to be non-Hermitian. These impedance matrices enable concise and efficient formulations of dispersion equations for wave guides, and solutions of scattering and related wave problems in cylinders.

## Introduction

Impedance provides a useful tool for solving dynamic problems in acoustics and elasticity. A single scalar impedance is usually sufficient in acoustics, whereas a matrix of impedance elements is required to handle the vector nature of elastic wave motion, particularly in the presence of anisotropy. The use of impedance matrices can offer new insight because their properties are intimately related to the fundamental physics of the problem, as, for instance, the Hermitian property of the impedance matrix which is directly linked to energy considerations. A classical example is Lothe & Barnett’s surface impedance matrix [1, 2] which proved to be crucial for understanding surface waves in anisotropic homogeneous half-spaces, with the result that it provides perhaps the simplest method for finding the Rayleigh wave speed. Biryukov [3, 4] has developed a general impedance approach for surface waves in inhomogeneous half-spaces based on the differential Riccati equation, see also [5]. Direct use of the impedance rather than the full displacement-traction wave field provides an efficient and stable procedure for computing high-frequency dispersion spectra. Several numerical schemes for guided waves and scattering in multilayered structures have been developed on this basis, e.g. [6, 7, 8]. These involve the 33 impedance matrix often called the ‘surface impedance’, although it actually differs from the 33 surface impedance of a half-space. It is useful to further distinguish the familiar 33 (‘conditional’) impedance from a 66 (‘two-point’) matrix more closely related to the matricant of the system equations. The nature of these impedances have been analyzed using the Stroh framework for homogeneous and functionally graded plates [9, 10]. It is noteworthy that both impedances are Hermitian under appropriate physical assumptions; however, their Hermiticity implies a somewhat different energy-flux property than the Hermiticity of the half-space impedance. A bibliography on the impedance matrices for piezoelectric media may be found in [4, 11].

The above review concerns rectangularly anisotropic materials and planar structures. The objective of this paper is to provide an equally comprehensive impedance formalism for time-harmonic modes of th azimuthal order in radially inhomogeneous cylindrically anisotropic materials of infinite axial extent and various circular configurations. An important element in this task is the Stroh-like state-vector formalism developed for such materials in [12]. The results of [12], which are based on the matricant in a Peano-series form (particularly the definition of the ‘two-point’ impedances similar to the case of planar structures) are however only relevant to a cylindrical annulus with no material around the central point The intrinsic singularity of elastodynamic solutions at the origin of the cylindrical coordinate system, which rules out the Peano series, is an essential distinguishing feature as compared to the Cartesian setup. The problem can be readily handled in (transversely) isotropic homogeneous media with explicit Bessel solutions; however, it becomes considerably more intricate for cylindrically anisotropic and for radially inhomogeneous solid cylinders. The main analytical tool in this case is the Frobenius series solution. The milestone results on its application to homogeneous and layered cylinders of various classes of cylindrical anisotropy include [13, 14, 15, 16, 17, 18] (see also the review [19]). The state-vector formalism based on the Frobenius solution for the general case of unrestricted cylindrical anisotropy and arbitrary radial variation of material properties [20] is of crucial importance to the present study. Another vital ingredient is the differential matrix Riccati equation for an impedance [4]. To the best of the authors’ knowledge, this equation has only recently been used for the first time in elasticity of cylinders by Destrade et al. [21] who numerically solved it for an elastostatic problem in tubes.

The presence of the special point distinguishes the solid-cylinder case from its Cartesian counterpart in many ways. Apart from the usual radiation condition at infinity, a similar kind of condition has to be applied at . The Riccati equation simultaneously determines the central-impedance at in a consistent manner while requiring it as the initial value for obtaining the solid-cylinder impedance. No other auxiliary (boundary) condition applies at which is (again) unlike the surface or conditional impedance for, say, a traction-free plane . These observations point to the fundamental role of the impedance formalism in cylindrically anisotropic elastodynamics and actually call for a new type of the impedance matrix appropriate for solid cylinders. The concept, properties and calculation of the solid-cylinder impedance are among the main results of this paper.

The outline is as follows. Background material on the matricant, impedance matrices and Riccati equations is presented in §1 in a general context not specific to cylindrical configurations. In §2 the governing equations for cylindrically anisotropic elastic solids are reviewed and the first order differential system for the displacement-traction vector is described. Some examples of the use of impedance matrices are discussed in §3, and in the process the solid cylinder and the radiation impedance matrices are introduced. Methods for determining the solid cylinder impedance are developed in §4. This section provides a detailed description of the Frobenius solution and its properties, and also discusses the Riccati solution. Both methods involve the crucial central impedance matrix, to which §5 is devoted, where explicit solutions are presented and general attributes delineated, including the important Hermitian property. The radiation impedance matrix is analyzed in §6. Explicit examples are presented in §7 for the central impedance matrix in different types of anisotropy, and the solid cylinder impedance is explicitly presented for transverse isotropy. Numerical results illustrating the Riccati equation solution method are also discussed. Concluding remarks can be found in §8.

## 1 The matricant, impedance matrices and Riccati equations

For the moment the development is independent of the physical dimension and the underlying coordinates. Consider a system of linear ordinary differential equations

(1) |

The dimensional vectors , and the submatrices , all possess uni-dimensional spatial dependence on , which may be a Cartesian or radial coordinate. The system matrix displays an important algebraic symmetry which is a consequence of a general flux continuity condition. The derivative of the scalar quantity , where superscript ‘+’ means the adjoint (complex conjugate transpose) and has block structure with zero submatrices on the diagonal and off-diagonal identity matrices, can be identified with the divergence of the flux vector (to be defined more specifically later). Thus, , and hence, (1) implies the connection between flux continuity and symmetry of the system matrix [10]

(2) |

The vanishing of assumes certain physical restrictions which will be described when the elasticity problem is considered in Section 2.

The matricant is a function of two coordinates defined as the solution of the initial value problem

(3) |

The matricant may be represented formally as a Volterra or multiplicative integral evaluated by means of a Peano series [22], alternatively it may be expanded in a Frobenius series [23]. Let be a set of partial solutions, that is, a complete set of independent solutions of the homogeneous system (1), then , where is the integral matrix (a first-rank tensor) . The propagator nature of the matricant is apparent from the property , and in particular , while the symmetry (2) implies

(4) |

Hence,

(5) |

that is, is -unitary [22].

In solving problems one is often not interested in the individual fields and , but rather in their relationship to one another, and perhaps only at one or two positions such as boundary values of . Accordingly we introduce the conditional impedance matrix defined such that

(6) |

The conditional nature of this impedance arises from an auxiliary condition at another coordinate [9, 10], and may be understood from an equivalent definition of the matricant

(7) |

Now suppose is the conditional impedance at , then

and the conditional impedance at is therefore

(8) |

In practice, is often associated with boundary conditions on the level surface . For instance, ‘zero traction’ and ‘rigid boundary’ conditions are specified by vanishing and , respectively, with conditional impedances

(9) |

While it is possible to define the conditional impedance in terms of solutions of the linear system (1), the same system leads through a process of elimination to a quadratically nonlinear equation for the matrix : the differential Riccati equation [4]

(10) |

In this context the auxiliary impedance serves as an initial condition at which once specified uniquely determines at other positions. The symmetry (2) renders eq. (10) self-adjoint in the sense that if is a solution then so is , which does not imply their equality. It does however imply that the differential Riccati equation (10) produces an Hermitian impedance, , as long as the initial condition is Hermitian, . We will also find useful the algebraic Riccati equation associated with eq. (10),

(11) |

the solution of which determines limiting values of the impedance, e.g. as , and can serve as the initial value for the differential equation (10).

We also introduce a ‘two-point’ impedance distinguished from the conditional impedance by its explicit dependence upon two arguments, and defined such that it relates the constituent parts of the -vector at and according to [9, 10]

(12) |

Comparing (6) and (12) one might be tempted to surmise that the two-point impedance is composed simply of block diagonal elements and identified as and , respectively, where is the conditional impedance, and with zero off-diagonal blocks ( and ). But the two-point impedance is more fundamental and thereby richer, as one can see by comparing (7) and (12), implying

where , . The identity (4) then implies the important properties that the two-point impedance is Hermitian, and that the matricant determinant is of unit magnitude, i.e.,

(14) |

It follows directly from and Jacobi’s formula that the phase satisfies the differential equation with initial condition . The matricant is therefore unimodular if vanishes. Further properties of the impedance may be deduced by swapping the ‘running’ and ’reference’ points and in (12) (i.e. in both , and ), implying the reciprocal form

whence follows an obvious relation

(15) |

The two-point impedance therefore has the structure

(16) |

As an alternative to eq. (8) the conditional impedance at may be expressed in terms of the impedance at by using the two-point impedance,

(17) |

where . Note that is Hermitian if is.

In the same way that the matricant satisfies an ordinary differential equation in , viz. , it is possible to express the dependence of on in differential form. Differentiating (12) with respect to and using (1) to eliminate the traction vectors yields an equation for the two-point impedance,

(18) |

The self-adjoint property of this equation is obvious because the two-point impedance is itself self-adjoint (Hermitian). Direct integration of the differential system (18) subject to initial conditions at is problematic because of the fact that all submatrices of are of the form as , and hence undefined. Differential equations with well defined (finite) initial value conditions can be obtained for the block matrices by simple manipulation of eq. (18), but we do not discuss this further here. It is interesting to note, however, that inspection of the block structure of (18) shows that the equation for decouples from the other submatrices and it is the same as the differential Riccati equation (10) for the conditional impedance (under the interchange ). Furthermore, since becomes unbounded as , eq. (18) implies that the submatrix is the conditional impedance with the auxiliary condition of rigid (infinite) impedance at , an observation that is verified by eqs. (9) and (1).

## 2 Cylindrically anisotropic elastic solids

### 2.1 Equations in cylindrical coordinates

The dynamic equilibrium equations for a linearly elastic material when expressed in cylindrical coordinates are [24]

(19) |

Here is the mass density, the displacement, and the traction vectors , , are defined by the orthonormal basis vectors of the cylindrical coordinates according to , where is the stress, and a comma denotes partial differentiation. With the same basis vectors, and assuming the summation convention on repeated indices, the elements of stress are where is the strain, are elements of the fourth order (anisotropic) elastic stiffness tensor, and denotes transpose. The traction vectors are [12]

where, in the notation of [1], the matrix has components for arbitrary vectors and . The explicit form of the various matrices is apparent with the use of Voigt’s notation

### 2.2 Cylindrically anisotropic materials

The concept of cylindrical anisotropy, which apparently originated with Jean Claude Saint-Venant, and has been elaborated by Lekhnitskii [25], demands the angular independence of material constants in the cylindrical coordinates, but admits their dependence on and as well on . We consider materials with no axial dependence whose density and the elasticity tensor may depend upon , i.e. and . We seek solutions in the form of time-harmonic cylindrical waves as

(21) |

where is the circumferential number.

The dependence of the displacement and traction on the single spatial coordinate allows the elastodynamic equations to be reduced to the canonical form of eq. (1) [12]:

(22) |

where is a vector

(23) |

and the system matrix is defined by

The individual matrices are

with the matrices

The constituent matrices are

where

The matrices and are negative definite and positive semi-definite, respectively, for real-valued and positive definite elastic moduli. Note that the th order modal solution is a function of the radial coordinate, but it is also an implicit function of the frequency and the axial wavenumber , which dependence is here kept tacit. In the same manner, the dependence of upon , and is understood.

The superscript is omitted henceforth, with the exception of the specific cases , , as required.

### 2.3 Cylindrical elasticity in the general context

The results of the previous subsection, particularly eqs. (22) and (23), show that the cylindrically anisotropic system of azimuthal order is a special case of the formulation of Section 1 generally with , and . The physical restrictions required for the Hermiticity condition (2) are real-valued , and material constants (more precisely, Hermitian elastic moduli suffice [11]). Under these conditions the matrix displays the symmetry

(24) |

The 66 matricant is the solution of the initial value problem

(25) |

The condition that and are strictly positive is important since the case of zero radial coordinate needs to be handled separately, which is discussed at length below. Note that we do not specify whether or is the greater or lesser of the two radii. The matricant allows us to express the state vector of partial modes in a cylinder as

(26) |

The pointwise elastodynamic energy balance is where is the energy density per unit volume and the energy flux vector. The pertinent form of for cylindrical elasticity is where is the time averaged radial component for azimuthal mode ,

(27) |

which together with the system equation (22) implies the symmetry (24) for (see also [12]).

The conditional impedance matrix relates traction and displacement at a particular value of , but specifically , according to eq. (6). The point requires a separate discussion, and indeed a newly defined impedance, introduced in the next section. For the moment we note that is contingent upon the definition of the (one-point) impedance at some radial coordinate, say . The traction at other values of is then unambiguously related to the local displacement by either the matricant or the two-point impedance matrices, using equation (8) or (17). By rewriting eq. (27) we see that the conditional impedance determines the pointwise flux,

(28) |

which is zero for all only if is Hermitian. This in turn is the case only if is Hermitian, i.e., if there is no flux across the surface . On the other hand, the 66 two-point impedance matrix of eq. (12) defines the global energy flow into or out of the finite region between the two radial coordinates . Let be the total energy in the shell cross-section per unit length of the cylinder for azimuthal mode . Its increment over one period of time harmonic motion is

(29) |

which is identically zero for real , and Hermitian parameters, i.e., when is Hermitian. If the material in the slab is lossy then should be positive definite in order that is not increasing with time.

The differential Riccati equation satisfied by the one-point impedance matrix follows from (10) as

(30) |

The initial value problem for is therefore

(31) |

where

(32) |

Equation (31) shows the explicit dependence upon , and the elastic moduli. The exclusion of the distinguished point at the cylinder centre is addressed next.

## 3 Wave impedance matrices for cylinders

In this section we describe typical uses of impedance matrices, and in the process introduce the solid-cylinder impedance and the radiation impedance . We consider the three distinct configurations depicted schematically in Figure 1.

### 3.1 Solid-cylinder impedance matrix

A solid cylinder, by definition, is one that includes the axis . A new impedance matrix is introduced to handle this situation. The solid-cylinder impedance is defined in the usual manner by its property of relating the traction and displacement vectors of (21) according to

(33) |

although this is not a conditional impedance matrix because of the absence of an auxiliary impedance condition at some other coordinate. Instead, the solidity of the cylinder at dictates the character of (and one could argue that it is ‘conditional’ in that sense). The limiting value of the solid impedance at plays a crucial role, and we accordingly define the central-impedance matrix:

(34) |

The properties of the central impedance are discussed in detail in §5 after we develop methods for finding the solid-cylinder impedance matrix in §4.

As an example application consider the task of finding the dispersion equation for guided waves of frequency and wavenumber . We suppose, quite generally, an interface condition on the level surface of the form

(35) |

where is considered as given. It could be zero (traction free condition), infinite (rigid boundary), or it could be defined by some surrounding material, whether finite or infinite in extent. For instance, if the solid cylinder is surrounded by a shell of cylindrically anisotropic material in lubricated contact at and free at , then where is the conditional impedance with the auxiliary condition . Assuming (35) describes the condition at the outer surface, the desired dispersion equation is

(36) |

It is instructive to compare (36) with the dispersion equation for a (possibly functionally graded) layer with the traction-free surface on a homogeneous substrate which may be written in the form [26]

where is a (constant) impedance of the substrate and is the conditional impedance of the layer satisfying the reference condition If the surrounding material beyond a rigid (say) interface is infinite, then there are Stoneley-like waves defined by the dispersion equation (35) with where is the radiation impedance discussed in §6.

The solid-cylinder impedance also provides a means to compute the modal displacement vector for all if the dispersion equation (36) is satisfied. By analogy with the case of an annulus [21], the unnormalized displacement follows from the system equation (22) and the definition of in (33) as the solution of the initial value problem

(37) |

where is the null vector of the surface impedance condition, . Note that the solution of (37) remains well behaved even as the matricant solution is numerically unstable, see §4.2.2.

### 3.2 Impedance matrices for cylinders of infinite radius

Consider a cylinder extending to infinity in the radial direction, with inner surface at , see Figure 1. A wave incident from results in a total field that can be expanded in terms of partial waves of the form (21). The amplitude of the azimuthal mode is

(38) |

where the scattered amplitude satisfies a radiation condition at . This in turn requires that the following condition prevails on the interface:

(39) |

where the radiation impedance matrix is defined by the radiation conditions, see §6. The scattered field is then uniquely determined by the condition at , which we assume is of the generalized form (35) with prescribed interface impedance . Then,

where is the impedance of the incident wave, which follows directly from the equations of motion. The scattered amplitude on the interface is therefore

(40) |

which provides the initial condition to determine the entire scattered field in . Further details on the radiation impedance matrix are provided in §6, including its asymptotic properties for large .

### 3.3 An annulus of finite thickness

The case of the annulus fits readily into the general theory. Again consider the task of finding the dispersion equation for guided waves, which may be found by simultaneous satisfaction of the conditions on the two radial surfaces. Suppose the conditions are both of the generalized form where are known quantities. The conditional impedance is determined (numerically) by integrating (25) from (say) with initial condition to give

(41) |

The interface condition at requires that

which implies the dispersion equation

(42) |

Variants on this equation may be obtained using the two-point impedance instead of the matricant. Thus, from eq. (16) we have the equivalent condition

The examples considered in this section illustrate the usefulness of wave impedance matrices for cylinders of finite and infinite radial extent. Solutions to problems of practical concern can be formulated concisely in terms of impedance matrices, such as the dispersion equation for guided waves, and the scattering of waves from a cylindrical region. Calculation of the impedance matrices is relatively straightforward using the matricant or two-point impedance matrices (see [12]), but only as long as the points or are not involved; otherwise the solid cylinder impedance and/or radiation impedance matrices are required. Determination of the solid-cylinder impedance matrix is discussed next.

## 4 The solid impedance matrix

In this section we develop methods to calculate the solid-cylinder impedance matrix for a radially inhomogeneous cylindrically anisotropic cylinder with material at . Two principle approaches are considered: a semi-explicit solution as a Frobenius series, and an implicit solution in terms of a differential Riccati equation.

Unlike the conditional impedance which can be determined directly from the matricant along with the prescribed reference value, the matricant is not of direct use here because of its divergence at . This introduces the need to identify ‘physical’ and ‘nonphysical’ constituents of the solution near , which is performed explicitly for the Frobenius solution. In the Riccati approach the displacement and traction fields are not considered explicitly and the divergence at is taken care of by the initial value of the impedance. The Frobenius solution is considered first.

### 4.1 Frobenius expansion

We take advantage of the fact that the fundamental solution can formally be written in terms of a Frobenius series, which is an explicit one-point solution valid at any (including ). As a result, the Frobenius-series approach provides a constructive definition of . The Frobenius series solution can be obtained via a recursive procedure with the number of numerically required terms increasing with . Before we present the formal solution for we review and develop some properties of the Frobenius series for cylindrically anisotropic materials, following the analysis in [20].

#### 4.1.1 Background material

The Frobenius solution is based on the integral matrix solution of Eq. (22), which can always be defined through the Frobenius series for any . The pivotal role in constructing this series belongs to the eigenspectrum of the matrix with the symmetry

(43) |

which follows from (24). Denote the eigenvalues and eigenvectors of by and (), and introduce the matrix . Barring extraordinary exceptions, if then (i) no two eigenvalues of differ by an integer, and (ii) all normally are distinct (and nonzero). Let us first consider this case , otherwise see §4.1.2. By virtue of (i), the integral matrix may be formulated as

(44) |

where is the Jordan form of the matrix which is diagonal when (ii) holds, and is defined recursively through [20, Eqs. (9)-(13)].

The arguments underlying eq. (2) imply that the matrix is a constant independent of , and according to eq. (27) this matrix defines the flux properties of the constituents , , , of . For the present purposes we wish to split them into a pair of triplets: a physical set and a nonphysical triplet , where the only non-zero flux interactions occur between and , thus ensuring the crucial property that has nonzero elements confined to the main diagonal of the off-diagonal blocks. The partitioning is accomplished through appropriate arrangement of the eigenspectrum of as

(45) |

[20, Eq. (44)]. Combining Eqs. (43) and (45) and adopting the normalization yields the orthogonality/completeness relation for the eigenvectors in the form

(46) |

It follows from Eqs. (43) through (46) that and hence the flux matrix at is ,

(47) |

Note that (46) yields .

In order to further clarify the structure of we represent the 66 matrices , and in terms of submatrices,

(48) |

where and is diagonal for . The integral matrix consequently has block structure

(49) |

Note in particular that the integral matrix consists of two distinct 63 matrices,

(50) |

the former with the columns tending to zero at and the latter with columns diverging at . The block structure of eqs. (46) and (47) is

The latter explicitly shows that the normal energy flux of the displacement-traction wave field comprising an arbitrary superposition of either the three modes or three modes with is zero at any . This specific arrangement of may be interpreted as the generalization of the isotropic case with solutions cast in terms of the cylinder functions and , corresponding to the physical and nonphysical triplets respectively, each of which yields zero flux individually. This partitioning will be crucial in developing an explicit solution for the solid impedance matrix.

#### 4.1.2 Overview of the cases and

Let us return to the two assumptions made above which are that (i) no two eigenvalues of differ by an integer and (ii) all are distinct, hence is semisimple (diagonalizable). Violating (i) invalidates the relatively simple form (44) of the Frobenius fundamental solution to the governing equation, see [23]. Violation of (ii), or more precisely, the occurrence of degenerate that makes non-semisimple, alters the orthogonality/completeness relations and the composition of given above for . The cases affected are (axisymmetric modes) and (lowest-order flexural modes): specifically, the property (i) does not hold for and the property (ii) does not hold for both and . From a physical point of view, the cases stand out because they are related to the rigid-body motions producing zero stresses [20, Eq. (19)]. Note also that admits a zero eigenvalue iff , [20, Eq. (30)] and that is always a double eigenvalue rendering non-semisimple.

Consider the axisymmetric case The six eigenvalues of are where for trigonal or tetragonal symmetry with [24, Eqs. (3.12), (3.13)].
It is seen that, whatever the symmetry, the
set of includes pairs different by
an integer. As a result, the integral matrix is now defined through in a rather intricate form elucidated in [20, Eqs. (A2), (A.4)].
This observation is essential for treating inhomogeneous
and low-symmetry homogeneous cylinders.
At the same time, if the cylinder is homogeneous and has orthorhombic
or higher symmetry with the exception of trigonal and tetragonal with ,
then decouples
into the solutions described by Bessel functions and/or by a simple
Frobenius form (44)^{3}^{3}3Orthorhombic or higher symmetry enables uncoupling of the pair of
torsional modes described by the Bessel solutions stemming from The four sagittal modes are associated with where for symmetry lower than the trigonal or tetragonal with . When so that the above quartet of involves pairs with an integer difference, the sagittal problem admits
explicit Bessel solutions for the isotropic or transverse isotropic symmetry
due to uncoupling of potentials. Note that double eigenvalues
at do not bring non-diagonal blocks into the Jordan form of
.

Consider the case The matrix has a doubly degenerate eigenvalue which makes non-semisimple [20, Eq. (36)]. This does not preclude taking in the form (44) but the matrix is now not diagonal. As a result, the triplet of physical modes (with one of the modes associated with ) retains its form (50), whereas the nonphysical triplet is no longer of the form (50) due to one of its modes involving both eigenvectors, the proper and the generalized ones and , associated with [20, Eqs. (51), (61)]. It is thus evident that the physical modes satisfy the same orthogonality/completeness relations as for ; moreover, subject to the optional condition , the nonphysical mod