# A numerical procedure for regularized integral equations for elastostatic problems of flat cracks in opening mode

Anh Le Van, Jean Royer

► To cite this version:

Anh Le Van, Jean Royer. A numerical procedure for regularized integral equations for elastostatic problems of flat cracks in opening mode. Engineering Fracture Mechanics, 1991, 40 (2), pp.237 - 250. 10.1016/0013-7944(91)90260-8. hal-01004924

HAL Id: hal-01004924

https://hal.science/hal-01004924

Submitted on 22 Oct 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A NUMERICAL PROCEDURE FOR REGULARIZED INTEGRAL EQUATIONS FOR ELASTOSTATIC PROBLEMS OF FLAT CRACKS IN OPENING MODE

A. LE VAN and J. ROYER

Laboratoire de Mécanique des Structures, Ecole Nationale Supérieure de Mécanique,
1, rue de la Noë, 44072 Nantes Cedex, France

**Abstract**—The regularized integral equations for flat cracks in opening mode are solved numerically by means of a simple procedure using new interpolation functions. The convergence and the accuracy of the method are tested on problems of circular and elliptical cracks embedded in the infinite medium, which admit analytical solutions. Computational results ensure the validity of this new approach, then they are compared with those already obtained without regularization by the authors.

## INTRODUCTION

FEW TRIDIMENSIONAL crack problems admit an analytical solution. Treated in the infinite medium case, they mainly concern cracks of particular geometries and subjected to particular loadings. The most well-known among these problems is that of the penny-shaped crack under an axisymmetrical pressure[1]. More recently, the case of more general loadings is fully studied in refs [2] and [3]. Like the penny-shaped crack, the elliptical crack also attracts many authors' attention because of its particular geometry. The elliptical crack subjected to a uniform pressure is treated in ref. [4], whereas the case of a linear pressure is studied in ref. [5] and the case of a polynomial pressure in ref. [6]. For elliptical cracks subjected to shear loadings, the main contributions may be found in ref. [7] and more recently in refs [8] and [9]. The rectangular crack, though not very realistic, is another crack geometry amenable to analytical solutions[10, 11].

For more general crack configurations, one often has to turn to numerical computations. Integral equations derived from the potential theory have the common feature of being singular in the sense of the principal value. The numerical treatment of the principal value is difficult; a recent successful attempt to deal with this singularity in two dimensions can be found in ref. [12]. However, computations without particular treatment of the singular part are able to yield highly satisfactory results as shown in ref. [13].

The regularization task aims at lessening the principal value singularity. In theoretical respects, this yields weakly singular integral equations, which must be more tractable than principal value ones. In numerical aspects, the kernel of regularized equations is found to be lengthy and its discretization in the three-dimensional case remains for the time being inextricable. Recently, a numerical procedure based on the least squares global displacement gradient smoothing method was worked out in ref. [14] for flat cracks in opening mode.

The present paper describes a different approach to problems of flat cracks in opening mode, using new interpolation functions denoted by $M^{(a)}$ and $j$. The proposed numerical procedure has the advantage of being simple: similar to the non-regularization method, the solution of the discretized equation is carried out by the collocation method and the integration classically by Gaussian points.

## DISCRETIZATION OF THE REGULARIZED INTEGRAL EQUATION

Consider a planar crack S embedded in the infinite medium and subjected to an arbitrary internal pressure (Fig. 1). It is known[15] that under these conditions the opening mode is uncoupled from the other ones. The crack surface can be parametrized with Cartesian coordinates: $D = S \ni (y_1, y_2) \to y \in S$. The normal displacement discontinuity $\phi_3 \equiv u_3^+ - u_3^-$ is then related to
<figure><img src="image_2.png" /><figcaption>Fig. 1. Planar crack. $y_0$: singular point; $\partial S$: crack edge.</figcaption></figure>

the third component of the stress vector with respect to the normal $n$ of the crack by the following regularized integral equation [16]:

$$ 
\begin{aligned}
\forall y_0 \in S, t_3(y_0, n) &\equiv -p(y_0) \\
&= \frac{E}{8\pi(1-v^2)} \left\{ \cdot \int_S \frac{1}{r^2} \left( \frac{\partial\phi_3}{\partial y_\alpha}(y) - \frac{\partial\phi_3}{\partial y_\alpha}(y_0) \right) r_\alpha \, dy_1 \, dy_2 - \left( \int_{\partial S} \frac{v_\alpha}{r} \, d\gamma \right) \frac{\partial\phi_3}{\partial y_\alpha}(y_0) \right\}.
\end{aligned} 
\quad (1) $$

(Implicit sums on Greek indices are made from 1 to 2.) Denoting by $S_{y_0}$ a neighbourhood of $y_0$ one can split up the surface integral in eq. (1) as follows:

$$ 
\begin{aligned}
\forall y_0 \in S, -p(y_0) = & \frac{E}{8\pi(1-v^2)} \left\{ \cdot \int_{S_{y_0}} \frac{1}{r^2} \left( \frac{\partial\phi_3}{\partial y_\alpha}(y) - \frac{\partial\phi_3}{\partial y_\alpha}(y_0) \right) r_\alpha \, dy_1 \, dy_2 \right. \\
& \left. + \int_{S \setminus S_{y_0}} \frac{1}{r^2} \frac{\partial\phi_3}{\partial y_\alpha}(y) r_\alpha \, dy_1 \, dy_2 - \left( \int_{S \setminus S_{y_0}} \frac{r_{\alpha\alpha}}{r^2} \, dy_1 \, dy_2 + \left( \int_{\partial S} \frac{v_\alpha}{r} \, d\gamma \right) \frac{\partial\phi_3}{\partial y_\alpha}(y_0) \right) \right\}
\end{aligned} 
\quad (2) $$

$v = \tau \wedge n$, $\tau$ is the vector tangent to the crack edge $\partial S$.

Let us introduce the following notations for the discretization of eq. (2).

In a general way, the symbol $\langle \rangle$ relates to a row-vector, whereas $\{ \}$ relates to a column-vector, and $[]$ to a matrix.

The geometry transformation will be denoted: $e_0 \ni \xi \to y \in e$, where $e_0$ is the reference element of the physical element $e$. The coordinates of point $y$ are expressed in the usual way: $\forall \alpha \in \{1, 2\}$, $y_\alpha = \langle N(\xi) \rangle \{y_\alpha\}^e, \langle N(\xi) \rangle$ designates the shape function.

The Jacobian matrix of the geometry transformation is denoted by

$$ [J(\xi)] = \left[ J_{a\beta}(\xi) = \frac{\partial y_{\beta}}{\partial \xi_{\alpha}}(\xi) \right]_{\{\alpha, \beta\} \in \{1, 2\}}, $$

its determinant by either $\det J(\xi)$ or $D(y)/D(\xi)$, and its inverse by

$$ [j(\xi)] = \left[ j_{a\beta}(\xi) = \frac{\partial \xi_{\beta}}{\partial y_{\alpha}}(\xi) \right]_{\{\alpha, \beta\} \in \{1, 2\}} $$

The last three integrals are quite regular and thus require no particular investigation. For instance if use is made of isoparametric elements, one has for any element $e$ lying in $S \setminus S_{y_0}$:

$$ \int_{e \in S \setminus S_{y_0}} \frac{1}{r^2} \frac{\partial \phi_3}{\partial y_\alpha}(y) r_{\alpha} \, dy_1 \, dy_2 = \int_{S_0} \frac{1}{r^2} \left\langle \frac{\partial N}{\partial \xi_\beta}(\xi) \right\rangle j_{a\beta}(\xi) r_{\alpha} \frac{D(y)}{D(\xi)} d\xi_1 \, d\xi_2 \cdot \{\phi_3\}^e. $$

It has been shown in ref. [16] that the first integral in eq. (2) is ordinary improper. So its study is less critical than that of a principal value integral integral. However, for numerical purposes an adequate integration scheme remains to be devised in order to take into account the weakly singular behaviour of the type $r^{-2+\gamma}, 0 < \gamma \le 1$. In this study, a numerical procedure will be proposed which makes it possible to solve this weak singularity.
<figure><img src="image_3.png" /><figcaption>Fig. 2. Polar coordinates (λ, θ) in the reference element e₀. Case of a quadrilateral element.</figcaption></figure>

First we define from $\langle N(\xi) \rangle$, $[J(\xi)]$ and $[j(\xi)]$, the following five functions which will be referred to as “modified”. The first of them was introduced in ref. [17] and widely used in refs [14] and [18]. For each couple of points $(\xi, \xi_0)$ in $e_0$, we will consider the polar coordinates $(\lambda, \theta)$ of $\xi$, related to centre $\xi_0$, defined by (Fig. 2):

$$ \xi_1 = \xi_{01} + \lambda \cos \theta $$

$$ \xi_2 = \xi_{02} + \lambda \sin \theta. $$

## Definition

The modified shape function $\langle \hat{N}(\xi_0, \lambda, \theta) \rangle$ is such that:

$$ \langle N(\xi) \rangle - \langle N(\xi_0) \rangle = \lambda \langle \hat{N}(\xi_0, \lambda, \theta) \rangle \quad (3a) $$

and: $\lim_{\lambda \to \infty} \langle \hat{N}(\xi_0, \lambda, \theta) \rangle < \infty$.

(The notation $A < \infty$ signifies that $A$ is bounded if $A$ is a real, that all components (resp. elements) of $A$ are bounded if $A$ is a vector (resp. a matrix)).

The modified derivated shape functions $\langle M^{(a)}(\xi_0, \lambda, \theta) \rangle$, $\alpha \in \{1, 2\}$, are such that:

$$ \left\langle \frac{\partial N}{\partial \xi_{\alpha}}(\xi) \right\rangle - \left\langle \frac{\partial N}{\partial \xi_{\alpha}}(\xi_0) \right\rangle = \lambda \langle M^{(a)}(\xi_0, \lambda, \theta) \rangle \quad (3b) $$

and: $\lim_{\lambda \to \infty} \langle M^{(a)}(\xi_0, \lambda, \theta) \rangle < \infty$.

The modified Jacobian matrix $[\hat{J}(\xi_0, \lambda, \theta)]$ is such that:

$$ [J(\xi)] - [J(\xi_0)] = \lambda [\hat{J}(\xi_0, \lambda, \theta)] \quad (3c) $$

and: $\lim_{\lambda \to \infty} [\hat{J}(\xi_0, \lambda, \theta)] < \infty$.

The modified Jacobian determinant $\hat{d}(\xi_0, \lambda, \theta)$ is such that:

$$ \det J(\xi) - \det J(\xi_0) = \lambda \hat{d}(\xi_0, \lambda, \theta) \quad (3d) $$

and: $\lim_{\lambda \to \infty} \hat{d}(\xi_0, \lambda, \theta) < \infty$.

The modified inverse of the Jacobian matrix $[\hat{J}(\xi_0, \lambda, \theta)]$ is such that:

$$ [j(\xi)] - [j(\xi_0)] = \lambda [f(\xi_0, \lambda, \theta)] \quad (3e) $$

and: $\lim_{\lambda \to \infty} [f(\xi_0, \lambda, \theta)] < \infty$.
For simple elements, e.g. six-node triangles and eight-node quadrilaterals, the existence of functions $\langle \hat{N} \rangle$ and $\langle M^{(a)} \rangle$, $\alpha \in \{1, 2\}$, can easily be verified, and their expressions are also readily available. Concerning the other modified functions, the following proposition shows how to determine them:

**Proposition**

If there exist $\langle M^{(a)}(\xi_0, \lambda, \theta) \rangle$, $\alpha \in \{1, 2\}$, then there also exist the functions $[\hat{J}(\xi_0, \lambda, \theta)]$, $\hat{d}(\xi_0, \lambda, \theta)$, $[\hat{J}(\xi_0, \lambda, \theta)]$ and they are obtained by:

$$ \forall \alpha, \beta \in \{1, 2\}, \quad \hat{J}_{a\beta}(\xi_0, \lambda, \theta) = \langle M^{(a)}(\xi_0, \lambda, \theta) \rangle \{y_{\beta}\}^{\epsilon} \quad (4a) $$

$$ \hat{d}(\xi_0, \lambda, \theta) = \hat{J}_{11}J_{22}(\xi) + \hat{J}_{22}J_{11}(\xi_0) - \hat{J}_{12}J_{21}(\xi) - \hat{J}_{21}J_{12}(\xi_0) \quad (4b) $$

$$ \forall \alpha, \beta \in \{1, 2\}, \quad \hat{J}_{a\beta}(\xi_0, \lambda, \theta) = \frac{(-1)^{\alpha + \beta}}{\det J(\xi) \det J(\xi_0)} [\hat{J}_{\gamma\delta} \det J(\xi_0) - \hat{d} \cdot J_{\gamma\delta}(\xi_0)] \quad (4c) $$

where $(\gamma, \delta)$ is equal to $(2, 2)$ if $(\alpha, \beta) = (1, 1)$, to $(1, 1)$ if $(\alpha, \beta) = (2, 2)$, and to $(\alpha, \beta)$ otherwise.

In the first or last two terms in eq. (4b), $\xi$ and $\xi_0$ are interchangeable. Clearly, $\hat{d} \neq \det \hat{J}$. The foregoing proposition, which can easily be shown, enables us to obtain systematically the expressions for $[\hat{J}]$, $\hat{d}$ and $[\hat{J}]$, provided the existence of $\langle M^{(a)} \rangle$, $\alpha \in \{1, 2\}$. Expressions for the modified function $\langle \hat{N} \rangle$ can be found in ref. [18]; for the sake of completeness, we give them in the Appendix, together with the new functions $\langle \hat{M}^{(a)} \rangle$ and $[\hat{J}]$. Now we divide the neighbourhood $S_{y_0}$ of the singular point into elements $e \ni y_0$, and transform the improper integral in eq. (2) over such an element as follows:

$$ \int_{e \in S_{y_0}} \frac{1}{r^2} \left\{ \left( \left( \frac{\partial \phi_3}{\partial \xi_\beta}(\xi) - \frac{\partial \phi_3}{\partial \xi_\beta}(\xi_0) \right) j_{a\beta}(\xi) + (j_{a\beta}(\xi) - j_{a\beta}(\xi_0)) \frac{\partial \phi_3}{\partial \xi_\beta}(\xi_0) \right) \{e, e_\alpha\} \frac{D(y)}{D(\xi)} d\xi_1 d\xi_2 \right. \quad (5) $$

$$ \text{From eq. (3a): } y_{\alpha} - y_{0\alpha} = (\langle N(\xi) \rangle - \langle N(\xi_0) \rangle) \{y_{\alpha}\}^{\epsilon} = \lambda \langle \hat{N}(\xi_0, \lambda, \theta) \rangle \{y_{\alpha}\}^{\epsilon} $$

$$ \text{hence: } r^2 = \sum_{\alpha=1}^{2} (y_{\alpha} - y_{0\alpha})^2 = \lambda^2 \cdot \sum_{\alpha=1}^{2} (\langle \hat{N}(\xi_0, \lambda, \theta) \rangle \{y_{\alpha}\})^2. \quad (6a) $$

Furthermore, the difference in the kernel in eq. (5) can be transformed, using eqs (3b) and (3c):

$$ \frac{\partial \phi_3}{\partial \xi_\beta}(\xi) - \frac{\partial \phi_3}{\partial \xi_\beta}(\xi_0) = (\langle N_\beta(\xi) \rangle - \langle N_\beta(\xi_0) \rangle) \{\phi_3\}^\epsilon = \lambda \langle M^{(\beta)}(\xi_0, \lambda, \theta) \rangle \{\phi_3\}^\epsilon \quad (6b) $$

$$ j_{a\beta}(\xi) - j_{a\beta}(\xi_0) = \lambda \hat{J}_{a\beta}(\xi_0, \lambda, \theta). \quad (6c) $$

Lastly, by noting that: $D(\xi_1, \xi_2)/D(\lambda, \theta) = \lambda$, the integral (5) can be rewritten hereafter where the square of polar radius $\lambda^2$ is now *cancelled* out:

$$ \int_{e \in S_{y_0}} = \left( \int_{(\lambda, \theta)} \frac{1}{\sum_{\alpha=1}^{2} (\langle \hat{N}(\xi_0, \lambda, \theta) \rangle \{y_{\alpha}\})^2} \times (\langle M^{(\beta)}(\xi_0, \lambda, \theta) \rangle j_{a\beta}(\xi) + \langle N_\beta(\xi_0) \rangle \hat{J}_{a\beta}(\xi_0, \lambda, \theta)) (e, e_\alpha) \frac{D(y)}{D(\xi)} d\lambda \, d\theta \right) \{\phi_3\}^\epsilon. \quad (7) $$

This time the integration is to be carried over $(\lambda, \theta)$ and eq. (7) involves only regular integrals which should yield no difficulties for the numerical purposes. It should be noticed that the common feature between the approaches in refs [14] and [18] and the method herein proposed is that regular integrals proceed merely from a numerical procedure. The kernel in eq. (5) actually remains (weakly) singular of the type $r^{-2+\gamma}$, $0 < \gamma \le 1$. Here the use of the interpolation type functions has made it possible to force the kernel to behave as a regular function in eq. (7).
<figure><img src="image_4.png" /><figcaption>Fig. 3. Change of variables: square $[-1, +1]^2 \ni v = (v_1, v_2) \rightarrow (\lambda, \theta) \in \text{subtriangle in } e_0$.</figcaption></figure>

In the following study, the discretization of eq. (5) is carried out considering only two usual elements: the six-node triangle and the eight-node quadrilateral. Nevertheless, the principle remains the same for any more general element.

### Use of eight-node quadrilaterals

The integration over $(\lambda, \theta)$ in the reference element $e_0$ is effected by dividing this element into subtriangles and by making another change of variables in each of these subtriangles. This operation will allow us to bring again the integration domain onto the square $[-1, +1]^2$ (Fig. 3).

According to whether the antecedent $\xi_0$ of point $y_0$ is located at a corner or on a side of element $e_0$, the subtriangle number is 2 or 3.

Case when $\xi_0$ is located at a corner of element $e_0$:

Let point $\xi_0$ be located at node $i \in \{1, 3, 5, 7\}$ of element $e_0$ (Fig. 4). There are two subtriangles and it can be checked for the expressions of the polar coordinates $(\lambda, \theta)$ and the Jacobian of the transformation $(v_1, v_2) \rightarrow (\lambda, \theta)$:

**Subtriangles of type 1:**

$$ \lambda = \frac{1 + v_1}{\cos\left(\frac{\pi}{8}(1 + v_2)\right)} $$

$$ \theta = \frac{\pi}{8}(v_2 + 2i - 1) $$

$$ \frac{D(\lambda, \theta)}{D(v)} = \frac{\pi}{8 \cos\left(\frac{\pi}{8}(1 + v_2)\right)} $$

**Subtriangles of type 2:**

$$ \lambda = \frac{1 + v_1}{\cos\left(\frac{\pi}{8}(1 - v_2)\right)} $$

$$ \theta = \frac{\pi}{8}(v_2 + 2i + 1) $$

$$ \frac{D(\lambda, \theta)}{D(v)} = \frac{\pi}{8 \cos\left(\frac{\pi}{8}(1 - v_2)\right)} $$

<figure><img src="image_5.png" /><figcaption>Fig. 4. $\xi_0$ located at a corner of a quadrilateral element.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Fig. 5. $\xi_0$ located at a mid-side of a quadrilateral element.</figcaption></figure>

Case when $\xi_0$ is located at a mid-side of $e_0$:

The four possible subcases are shown in Fig. 5. The node number of $\xi_0$ is $i \in \{2, 4, 6, 8\}$. Here there are three subtriangles. By putting $\theta_0 = \tan^{-1}(2)$, the following relations can easily be established:

Subtriangles of type 1:

$$ \lambda = \frac{1 + v_1}{2 \cos\left(\frac{\theta_0}{2}(1 + v_2)\right)} $$

$$ \theta = \frac{\theta_0}{2}(v_2 + 1) + \left(\frac{i}{2} - 1\right) \frac{\pi}{2} $$

$$ \frac{D(\lambda, \theta)}{D(v)} = \frac{\theta_0}{4 \cos\left(\frac{\theta_0}{2}(1 + v_2)\right)} $$

Subtriangles of type 2:

$$ \lambda = \frac{1 + v_1}{\cos\left(\left(\frac{\pi}{2} - \theta_0\right)v_2\right)} $$

$$ \theta = \frac{1}{2} \left[ (\pi - 2\theta_0)v_2 + \frac{i}{2}\pi \right] $$

$$ \frac{D(\lambda, \theta)}{D(v)} = \frac{\pi - 2\theta_0}{2 \cos\left(\left(\frac{\pi}{2} - \theta_0\right)v_2\right)} $$

Subtriangles of type 3:

$$ \lambda = \frac{1 + v_1}{2 \cos\left(\frac{\theta_0}{2}(1 - v_2)\right)} $$

$$ \theta = \frac{\theta_0}{2}(v_2 - 1) + \left(\frac{i}{2} + 1\right) \frac{\pi}{2} $$

$$ \frac{D(\lambda, \theta)}{D(v)} = \frac{\theta_0}{4 \cos\left(\frac{\theta_0}{2}(1 - v_2)\right)} $$

Use of six-node triangles

Here use is made of the same technique of dividing element $e_0$ into subtriangles and making the change of variables $[-1, +1]^2 \ni v = (v_1, v_2) \to (\lambda, \theta)$. If $\xi_0$ is at a corner (resp. at a mid-side) of $e_0$, the number of subtriangles is 1 (resp. 2). For triangular elements, there is no general rule and the six possible cases are to be considered separately. We recall the notation $\theta_0 = \tan^{-1}(2)$.

Case of $\xi_0$ at node 1 (Fig. 6a)

$$ \lambda = \frac{\sqrt{2(1 + v_1)}}{4 \cos \frac{\pi v_2}{4}} $$

$$ \theta = \frac{\pi}{4}(v_2 + 1) $$

$$ \frac{D(\lambda, \theta)}{D(v)} = \frac{\pi \sqrt{2}}{16 \cos \frac{\pi v_2}{4}} $$
<figure><img src="image_7.png" /><figcaption>Fig. 6. Positions of $ξ_0$ in a six-node triangle. Definition of $λ$ and $θ$.</figcaption></figure>

**Case of $ξ_0$ at node 2 (Fig. 6b) in subtriangle 1**

in subtriangle 2

$$ \lambda = \frac{1 + v_1}{4\sqrt{2 \cos\left(\theta - \frac{\pi}{4}\right)}} $$

$$ \theta = \frac{1}{2}(\pi - \theta_0)(1 + v_2) \qquad \theta = \frac{\theta_0}{2}(v_2 - 1) + \pi $$

$$ \frac{D(\lambda, \theta)}{D(v)} = \frac{\pi - \theta_0}{8\sqrt{2 \cos\left(\theta - \frac{\pi}{4}\right)}} $$

$$ \frac{D(\lambda, \theta)}{D(v)} = \frac{\theta_0}{8 \cos(\pi - \theta)} $$

**Case of $ξ_0$ at node 3 (Fig. 6c)**

$$ \lambda = \frac{1 + v_1}{2 \cos\left(\frac{\pi}{8}(1 - v_2)\right)} $$

$$ \theta = \frac{\pi}{8}(v_2 + 7) $$

$$ \frac{D(\lambda, \theta)}{D(v)} = \frac{\pi}{16 \cos\left(\frac{\pi}{8}(1 - v_2)\right)} $$
Case of $\zeta_0$ at node 4 (Fig. 6d)
in subtriangle 1

in subtriangle 2

$$
\lambda = \frac{1 + v_1}{4 \cos \frac{\pi v_2}{4}}
$$

$$
\lambda = \frac{1 + v_1}{4 \cos \frac{\pi v_2}{4}}
$$

$$
\theta = \frac{\pi}{4}(v_2 + 4)
$$

$$
\theta = \frac{\pi}{8}(v_2 + 6)
$$

$$
\frac{D(\lambda, \theta)}{D(v)} = \frac{\pi}{16 \cos \frac{\pi v_2}{4}}
$$

$$
\frac{D(\lambda, \theta)}{D(v)} = \frac{\pi}{16 \cos \frac{\pi v_2}{4}}
$$

Case of $\xi_0$ at node 5 (Fig. 6e)

$$
\lambda = \frac{1 + v_1}{2 \cos\left(\frac{\pi}{8}(1 + v_2)\right)}
$$

$$
\theta = \frac{\pi}{8}(v_2 + 13)
$$

$$
\frac{D(\lambda, \theta)}{D(v)} = \frac{\pi}{16 \cos\left(\frac{\pi}{8}(1 + v_2)\right)}
$$

Case of $\xi_0$ at node 6 (Fig. 6f)
in subtriangle 1

in subtriangle 2

$$
\lambda = \frac{1 + v_1}{4\sqrt{2 \cos\left(\theta - \frac{\pi}{4}\right)}}
$$

$$
\lambda = \frac{1 + v_1}{4 \cos \left( \theta - \frac{3\pi}{2} \right)}
$$

$$
\theta = \frac{1}{2}[(\pi - \theta_0)v_2 + \theta_0]
$$

$$
\theta = \frac{\theta_0}{2}(v_2 + 1) + \frac{3\pi}{2}
$$

$$
\frac{D(\lambda, \theta)}{D(v)} = \frac{\pi - \theta_0}{8\sqrt{2} \cos\left(\theta - \frac{\pi}{4}\right)}
$$

$$
\frac{D(\lambda, \theta)}{D(v)} = \frac{\theta_0}{8 \cos\left(\frac{\theta_0}{2}(1 + v_2)\right)}
$$

Discretized regularized equation

Eventually, the integral in eq. (5) assumes the following form:

$$
\begin{equation}
\begin{split}
\int_{e \in S_{y_0}} \frac{1}{r^2} \left( \frac{\partial \phi_3}{\partial y_\alpha} (y) - \frac{\partial \phi_3}{\partial y_\alpha} (y_0) \right) r_\alpha \, dy_1 \, dy_2 = \sum_{\text{v-squares}} \int_{(v_1, v_2)} \frac{1}{\sum_{\alpha=1}^2 \left( \left\langle \hat{N}(\xi_0, \lambda, \theta) \right\rangle \{y_\alpha\} \right)^2} \\
& \quad \times \left( \langle M^{(\beta)}(\xi_0, \lambda, \theta) \rangle j_{\alpha\beta}(\xi) + \langle N_{,\beta}(\xi_0) \rangle j_{\alpha\beta}(\xi_0, \lambda, \theta) \right) (e, e_\alpha) \frac{D(y)}{D(\xi)} \frac{D(\lambda, \theta)}{D(v)} \, dv_1 \, dv_2 \cdot \{\phi_3\}^e.
\end{split}
\tag{8}
\end{equation}
$$

The integration over the variables ($v_1, v_2$) is made in the $v$-squares $[-1, +1]^2$. Since the kernels involved are regular, the Gaussian scheme for numerical integrations is quite acceptable. Figure 7 presents the successive steps to pass from the former square $[-1, +1]^2 \ni (v_1, v_2)$ to the final physical element in which the point $y = (y_1, y_2)$ lies.
<figure><img src="image_8.png" /><figcaption>Fig. 7. Change from $(v_1, v_2)$ to $(y_1, y_2)$. Case when the element is a quadrilateral and $y_0$ lies at a corner of element $e$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 8. Coarsest and finest circular crack meshes and the computational results on $\phi_3$.</figcaption></figure>
Equation (8) together with eq. (2) constitutes the discretized regularized equation for a flat crack in opening mode. As remarked above, the weak singularity in eq. (2) has been removed with the aid of functions $\tilde{N}$, $M^{(a)}$ and $j$; of the numerical scheme there remains only regular integrals.

# NUMERICAL PARTS

## Convergence test: the circular crack problem

Consider the problem of a circular crack centred at O, with radius $a$, imbedded in the infinite medium and subjected to a uniform internal pressure $p$.

Like any flat crack in the infinite medium, the opening mode is uncoupled from the other modes. Four crack mesh configurations are chosen, containing respectively 85, 161, 209 and 321 interior nodes; the first and the last of them are presented in Fig. 8. Use is made of both six-node triangular and eight-node quadrilateral elements; those adjacent to the crack edge are quarter-point elements[19, 20]. All elements have straight edges; apart from the corner and quarter-point nodes the others are located at the mid-side of the edges.

For elements not containing the collocation point (i.e. lying in $S \setminus S_{0}$, see notations in eq. (2)), the integration is performed with 7 points over triangles and 9 points over quadrilaterals. Integration over elements containing the collocation point must be carried out using the change of variables which leads to a square integration region in the plane $(v_{1}, v_{2})$ (Fig. 7). The Gaussian quadrature scheme requires $3 \times 3 = 9$ points. Here the number of Gaussian points is twice as high as that used in the same study without regularization[13]. This is accounted for by the fact that the original physical element is more “distorted” through the transformation (see Fig. 7) leading to the square element in the $v$-plane.

As expected, the numerical results show a revolution symmetry for the values of the displacement discontinuity $\phi_{3}$ (Fig. 8). It should be noticed that the Young modulus $E$ datum is unnecessary for computing the stress intensity factor $K_{1}$. On the contrary, the Poisson ratio $\nu$ must be specified; here it is set to be equal to 0.3.

Table 1 shows minimal and maximal errors on the displacement discontinuity $\phi_{3}$ as well as errors on the stress intensity factor $K_{1}$, with respect to the well-known analytical solution[1]:

$$ \frac{E\phi_3(\rho)}{8\pi(1-\nu^2)} = \frac{pa}{\pi^2} \sqrt{[1-(\rho/a)^2]}, \quad \rho = |Oy|; \quad \text{which yields: } K_1 = 2p\sqrt{a/\pi}. $$

As remarked in the study of notched strips in ref.[19], although the stiffness of quarter-point elements is reasonably accurate, the local values of stress and displacement in these elements are poor, thus our error estimations on $\phi_{3}$ do not take into account quarter-point node values. Furthermore, the $\phi_{3}$-error is not maximal at the crack centre but at some radial distance varying according to the mesh configuration; in present cases it remains about 5% in magnitude. As regards the $K_{1}$-error, it markedly decreases as the mesh becomes denser. On the strength of these results, another investigation concerning an elliptical crack is treated in the next part.

## Elliptical cracks

Consider an elliptical crack of semi-axes $a$, $b$, imbedded in the infinite medium and subjected to a uniform internal pressure $p$. Figure 9 shows the meshes for two cracks with semi-axes ratio $a/b$ respectively equal to 1.5 and 2; they include 161 (resp. 113) interior nodes and 64 (resp. 48) elements. Here use is also made of quarter-point elements along the crack edge.

<table><caption>Table 1. Circular crack. Errors on $\phi_{3}$ and $K_{1}$, with respect to the analytical solution</caption><thead><tr><th>Number of d.o.f.</th><th>Minimal errors on φ<sub>3</sub></th><th>Maximal errors on φ<sub>3</sub></th><th>Errors on K<sub>1</sub></th></tr></thead><tbody><tr><td>85</td><td>0.28%</td><td>5.12%</td><td>6.06%</td></tr><tr><td>161</td><td>-0.05%</td><td>3.05%</td><td>2.29%</td></tr><tr><td>209</td><td>0.42%</td><td>-4.47%</td><td>0.78%</td></tr><tr><td>321</td><td>-0.02%</td><td>-4.72%</td><td>0.11%</td></tr></tbody></table>
<figure><img src="image_10.png" /><figcaption>Fig. 9. Elliptical crack meshes and the computational results on $\phi_3$.</figcaption></figure>

The mesh generation technique is rudimentary: the crack is divided up with elliptical curves shaped in the crack contour, and radial straight lines issued from the crack centre. This inevitably gives rise to some distorted (elongated) elements which behave badly in numerical purposes. Moreover, this does not allow an accurate representation of the crack edge, whenever the semi-axes ratio becomes large. Figure 9 clearly shows that the mesh is particularly ill-suited to shape the vertex of the major axis. These reasons explain why computations are limited to small semi-axes ratio geometries. They also justify the use of 12 (resp. 16) Gaussian points on triangles (resp. quadrilaterals). Numerical results will be compared with the theoretical solutions [4]:

$$ \frac{E\phi_3(y)}{8\pi(1-v^2)} = \frac{pb}{2\pi E(k)} \sqrt{[1-(y_1/a)^2-(y_2/b)^2]} $$

where $E(k) = \int \sqrt{[1 - k^2 \sin^2 \varphi]} d\varphi$, $k^2 = 1 - b^2/a^2$, denotes the complete elliptic integral of the second kind.

$$ K_1(A) = K_{\text{Imax}} = p\sqrt{(\pi b)/E(k)} $$

$$ K_1(B) = K_{\text{Imin}} = p\sqrt{(\pi b^2/a)/E(k)}. $$

Table 2 gives minimal and maximal errors on the displacement discontinuity $\phi_3$ together with errors on the stress intensity factor $K_1$ at the vertices A, B of the semi-axes (Fig. 9).
<table><caption>Table 2. Elliptical crack. Errors on $\phi_3$ and $K_1$</caption><thead><tr><th>a/b</th><th>Number of d.o.f.</th><th>Minimal errors on $\phi_3$</th><th>Maximal errors on $\phi_3$</th><th>Errors on $K_1(A)$</th><th>Errors on $K_1(B)$</th></tr></thead><tbody><tr><td>1.5</td><td>161</td><td>-0.10%</td><td>-5.36%</td><td>4.21%</td><td>-2.25%</td></tr><tr><td>2.0</td><td>113</td><td>0.16%</td><td>8.22%</td><td>8.71%</td><td>3.13%</td></tr></tbody></table>

As predicted, the values of the stress intensity factor at vertex B are always worse than those at vertex A, since the mesh is not smooth enough at B, especially in the second mesh. The maximal $\phi_3$-error in the last mesh occurs at a point of the major axis, close to vertex B; it undoubtedly contributes to the large error on $K_1(A)$. In order to improve the results, more sophisticated mesh configurations must be devised, which result in an adequate refinement at the crack front and at the same time less elements at the crack centre where there is a small displacement gradient.

## CLOSURE

The results show that the regularization method allows a high convergence and is quite operational for plane crack problems. Nevertheless, when computing rectangular cracks, for unknown reasons the results obtained by the authors are too poor to be published here. For circular cracks, the errors on $K_1$ are quickly reduced to less than 1%. For elliptical cracks, despite a coarse mesh and a bad crack contour representation, $K_1$-values remain quite acceptable where they are larger. Obviously, the circular mesh type (with radial rays) leads to too many elements at the crack centre without enough elements along the crack front. Thus more sophisticated mesh generation techniques are required to obtain better results.

Moreover, compared to the classical transformation $(\xi_1, \xi_2) \to (y_1, y_2)$, the transformation $(v_1, v_2) \to (y_1, y_2)$ entails a greater distortion of the physical element and implies a denser integration scheme as a compensation. The Gaussian point number used here is twice or four times as high as that for the method without regularization (9 or 16 points instead of 4). The above distortion may also account for the fact that the elements seem to be very sensitive to their aspect ratio (ratio of their longer to shorter side): the longer elements are, the worse the results. In connection with the non-regularization method, there appear two questions which call for answers.

Can one avoid using quarter-point elements? Indeed, the use of quarter-point elements condemned us to the range of six-node triangles and eight-node quadrilaterals. The non-regularization method, as shown in ref.[13], gave good results without using special elements.

The regularized eq. (8) was obtained and used in Cartesian coordinates. Is it possible to obtain a regularized equation in other coordinate systems? With reference to the discretized equation for the non-regularization method, it enabled us to cope with more general coordinates and in particular to compute the stress intensity factor on a local basis.

Under present conditions, the results obtained with regularization are encouraging and should encourage investigations with an improved mesh generation technique, and in a further stage to numerical developments with the view to studying three-dimensional cracks.

## REFERENCES

[1] I. N. Sneddon, The distribution of stress in the neighbourhood of a crack in an elastic solid. Proc. R. Soc. London, Series A 187, 229-260 (1946).

[2] J. T. Guiders and R. W. Lardner, Penny-shaped cracks. J. Elasticity 5, 59-73 (1975).

[3] J. C. Bell, Stresses from arbitrary loads on a circular crack. Int. J. Fracture 15, 85-104 (1979).

[4] A. E. Green and I. N. Sneddon, The distribution of stress in the neighbourhood of a flat crack in an elastic solid. Proc. Camb. Phil. Soc. Math. Phys. Sci. 46, 159-163 (1950).

[5] M. K. Kassir and G. C. Sih, Geometric discontinuities in elastostatics. J. Math. Mech. 16, 927-948 (1967).

[6] R. C. Shah and A. S. Kobayashi, Stress intensity factor for an elliptical crack under arbitrary normal loading. Engng Fracture Mech. 3, 71-96 (1971).

[7] M. K. Kassir and G. C. Sih, Three-dimensional stress distribution around an elliptical crack under arbitrary loading. J. appl. Mech. 60-611 (1966).

[8] F. W. Smith and D. R. Sorensen, The elliptical crack subjected to non uniform shear loading. J. appl. Mech., Series E 41, 502-506 (1974).

[9] A. N. Borodachev, Elliptical crack under the action of shearing stresses of polynomial form. Prikl. Mekh. 18, 84-87 (1983).
[10] M. K. Kasir, Stress intensity factor for three-dimensional rectangular cracks. *J. appl. Mech.* **48**, 309-312 (1981).

[11] M. K. Kasir, Three-dimensional rectangular cracks subjected to shear loading. *Int. J. Solids Struct.* **18**, 1075-1082 (1982).

[12] V. Sladek and J. Sladek, The calculation of singular integrals in the boundary integral formulation of 2D elastostatics. *Engng Anal.* **3**, 25-35 (1986).

[13] A. Le van and B. Peseux, Boundary element analysis of an integral equation for three-dimensional crack problems. *Int. J. numer. Meth. Engng* **26**, 2383-2402 (1988).

[14] E. Z. Polch, T. A. Cruse and C. J. Huang, Traction B.I.E. solutions for flat cracks. *Comput. Mech.* **2**, 253-267 (1987).

[15] H. D. Bui, Applications des potentiels élastiques à l'étude des fissures planes de forme arbitraire en milieu tridimensionel. *Comptes Rend. Acad. Sci.*, Paris, t.280, Series A, 1157-1160 (1975).

[16] A. Levan and J. Royer, Theoretical basis of regularized integral equations for elastostatic crack problems. *Int. J. Fracture* **44**, 155-166 (1990).

[17] F. J. Rizzo and D. J. Shippy, An advanced BIEM for three-dimensional thermoelasticity. *Int. J. numer. Meth. Engng* **11**, 1753-1768 (1977).

[18] M. Bonnet, Méthode des équations intégrales régularisées en élastodynamique. *Thesis*. Ecole Nationale des Ponts et Chaussées (November 1986).

[19] R. D. Henshell and K. G. Shaw, Crack tip elements are unnecessary. *Int. J. numer. Meth. Engng* **9**, 495-509 (1975).

[20] R. S. Barsoum, On the use of isoparametric finite elements in linear fracture mechanics. *Int. J. numer. Meth. Engng* **10**, 25-38 (1976).

# APPENDIX

For the six-node triangle and the eight-node quadrilateral, we give functions $\langle N \rangle$ together with the modified functions $\langle \hat{N} \rangle$, $\langle M^{(a)} \rangle$ and $\langle J \rangle$. Function $\langle J \rangle$ can be obtained from relations (4).

## Six-node triangular elements

<table><caption>Table 1A. Function $N$</caption><thead><tr><th></th><th>N<sub>1</sub>(ξ)</th></tr></thead><tbody><tr><td>1</td><td>-ξ'(1 - 2ξ')</td></tr><tr><td>2</td><td>4ξ'ξ<sub>1</sub></td></tr><tr><td>3</td><td>-ξ<sub>1</sub>(1 - 2ξ<sub>1</sub>)</td></tr><tr><td>4</td><td>4ξ<sub>1</sub>ξ<sub>2</sub></td></tr><tr><td>5</td><td>-ξ<sub>2</sub>(1 - 2ξ<sub>2</sub>)</td></tr><tr><td>6</td><td>4ξ'ξ<sub>2</sub></td></tr></tbody></table>

where $ξ' = 1 - ξ_1 - ξ_2$.

<table><caption>Table 2A. Function $\hat{N}$</caption><thead><tr><th></th><th>\hat{N}_i(ξ<sub>0</sub>, λ, θ)</th></tr></thead><tbody><tr><td>1</td><td>(1 - 4ξ'<sub>0</sub>)(cos θ + sin θ) + 2λ(cos θ + sin θ)<sup>2</sup></td></tr><tr><td>2</td><td>4(ξ'<sub>0</sub> cos θ - ξ<sub>01</sub>(cos θ + sin θ)) - 4λ cos θ (cos θ + sin θ)</td></tr><tr><td>3</td><td>(4ξ'<sub>01</sub> - 1)cos θ + 2λ cos<sup>2</sup> θ</td></tr><tr><td>4</td><td>4(ξ'<sub>02</sub> cos θ + ξ<sub>01</sub> sin θ) + 4λ cos θ sin θ</td></tr><tr><td>5</td><td>(4ξ'<sub>02</sub> - 1)sin θ + 2λ sin<sup>2</sup> θ</td></tr><tr><td>6</td><td>4(ξ'<sub>0</sub> sin θ - ξ<sub>02</sub>(cos θ + sin θ)) - 4λ sin θ (cos θ + sin θ)</td></tr></tbody></table>

where $ξ' = 1 - ξ_01 - ξ_02$.

By noting that: $ξ' - ξ_0' = -λ(cos θ + sin θ)$, we obtain Table 3A.

<table><caption>Table 3A. Function $M^{(a)}$</caption><thead><tr><th></th><th>M<sup>(1)</sup>(ξ<sub>0</sub>, λ, θ)</th><th>M<sup>(2)</sup>(ξ<sub>0</sub>, λ, θ)</th></tr></thead><tbody><tr><td>1</td><td>4(cos θ + sin θ)</td><td>4(cos θ + sin θ)</td></tr><tr><td>2</td><td>-4(2 cos θ + sin θ)</td><td>-4 cos θ</td></tr><tr><td>3</td><td>4 cos θ</td><td>0</td></tr><tr><td>4</td><td>4 sin θ</td><td>4 cos θ</td></tr><tr><td>5</td><td>0</td><td>4 sin θ</td></tr><tr><td>6</td><td>-4 sin θ</td><td>-4(cos θ + 2 sin θ)</td></tr></tbody></table>

Here the $M^{(a)}$ are independent of $λ$.

## The modified Jacobian inverse $\hat{J}$

Application of relations (4) gives, after performing all calculations:

$$[J(\xi_0, \lambda, \theta)] = \frac{1}{\det J(\xi) \det J(\xi_0)} \begin{bmatrix} A_{apq} & A_{-d, -e, -f} \\ A_{-b, -g, -h} & A_{abc} \end{bmatrix}$$
where for any triplet $(r, s, t) \in \{a, b, c, d, e, f, g, h, p, q\}^3$:

$$
\begin{aligned}
A_{rrt} ={}& (ae - bd)[-(r\xi_1 - \xi_0)\cos\theta + s(\sin\theta - \xi_1^2 - 2\xi_0\xi_2\cos\theta - \lambda\xi_0\cos\theta) - t(2\xi_0 + \lambda\cos\theta)\cos\theta] \\
& + (bp - ge)[r(\cos\theta - \xi_2^2 - 2\xi_0\xi_2\sin\theta - \lambda\xi_0\sin\theta) - s\xi_2\xi_2\sin\theta - t(2\xi_2 + \lambda\sin\theta)\sin\theta] \\
& + (ap - gd)[-(r\xi_1 - \xi_0)\sin\theta - s\xi_0\xi_2\cos\theta - t(\xi_0\sin\theta + \xi_2\cos\theta + \lambda\sin\theta)\cos\theta] \\
& + (aq + ce - bf - hd)[s(-\xi_2\cos\theta + \xi_0\sin\theta) - t\cos\theta] \\
& + (bq + cp - gf - he)[r(\xi_2\cos\theta - \xi_0\sin\theta) - t\sin\theta] \\
& + (cq - hf)[r\cos\theta + s\sin\theta]
\end{aligned}
$$

with:

$$a = 4(y_1 - 2y_1^2 + y_1^3)$$

$$d = 4(y_1^2 - 2y_1^2 + y_1^3)$$

$$b = 4(y_1 - y_1^2 + y_1^3 - y_1^4)$$

$$e = 4(y_1^2 - y_1^2 + y_1^3 - y_1^4)$$

$$c = -3y_1^2 + 4y_1^2 - 3y_1^3$$

$$f = -3y_1^2 + 4y_1^2 - 3y_1^3$$

$$g = 4(y_1 + y_1^2 - 2y_1^3)$$

$$p = 4(y_1^2 + y_1^3 - 2y_1^4)$$

$$h = -3y_1^2 - y_1^3 + 4y_1^4$$

$$q = -3y_1^2 - y_1^3 + 4y_1^4$$

Eight-node quadrilaterals

<table><caption>Table 4A. Function N</caption><thead><tr><th></th><th>N<sub>1</sub>(ξ)</th></tr></thead><tbody><tr><td>1</td><td>-(1 - ξ<sub>1</sub>)(1 - ξ<sub>2</sub>)(1 + ξ<sub>1</sub> + ξ<sub>2</sub>)</td></tr><tr><td>2</td><td>(1 - ξ<sub>1</sub>)(1 - ξ<sub>2</sub>)</td></tr><tr><td>3</td><td>-(1 + ξ<sub>1</sub>)(1 - ξ<sub>2</sub>)(1 - ξ<sub>1</sub> + ξ<sub>2</sub>)</td></tr><tr><td>4</td><td>(1 + ξ<sub>1</sub>)(1 - ξ<sub>2</sub>)</td></tr><tr><td>5</td><td>-(1 + ξ<sub>1</sub>)(1 + ξ<sub>2</sub>)(1 - ξ<sub>1</sub> - ξ<sub>2</sub>)</td></tr><tr><td>6</td><td>(1 - ξ<sub>1</sub>)(1 + ξ<sub>2</sub>)</td></tr><tr><td>7</td><td>-(1 - ξ<sub>1</sub>)(1 + ξ<sub>2</sub>)(1 + ξ<sub>1</sub> - ξ<sub>2</sub>)</td></tr><tr><td>8</td><td>(1 - ξ<sub>1</sub>)(1 - ξ<sub>2</sub>)</td></tr></tbody></table>

<table><caption>Table 5A. Function N̂</caption><thead><tr><th></th><th>N̂<sub>1</sub>(ξ<sub>0</sub>, λ, θ)</th></tr></thead><tbody><tr><td>1</td><td>{-cos θ(a<sub>1</sub> + b<sub>1</sub>)a<sub>2</sub> - sin θ(a<sub>2</sub> + b<sub>1</sub>)a<sub>1</sub> + λ[a<sub>2</sub> cos<sup>2</sup> θ + (a<sub>1</sub> + a<sub>2</sub> + b<sub>1</sub>)cos θ sin θ + a<sub>1</sub> sin<sup>2</sup> θ] - λ<sup>2</sup> cos θ sin θ(cos θ + sin θ)}</td></tr><tr><td>2</td><td>{a<sub>2</sub>(a<sub>1</sub> - a<sub>3</sub>)cos θ - a<sub>1</sub>a<sub>3</sub> sin θ + λ[(a<sub>3</sub> - a<sub>1</sub>)cos θ sin θ - a<sub>2</sub> cos<sup>2</sup> θ] + λ<sup>2</sup> cos<sup>2</sup> θ sin<sup>2</sup> θ}</td></tr><tr><td>3</td><td>obtained from N̂<sub>1</sub>, N̂<sub>2</sub> by changing: (a<sub>1</sub>, a<sub>2</sub>, a<sub>3</sub>) into (a<sub>2</sub>, a<sub>3</sub>, a<sub>4</sub>)</td></tr><tr><td>4</td><td>b<sub>1</sub> into b<sub>2</sub> (cos θ, sin θ) into (sin θ, -cos θ)</td></tr><tr><td>5</td><td>obtained from N̂<sub>1</sub>, N̂<sub>2</sub> by changing: (a<sub>1</sub>, a<sub>2</sub>, a<sub>3</sub>) into (a<sub>3</sub>, a<sub>4</sub>, a<sub>1</sub>)</td></tr><tr><td>6</td><td>b<sub>1</sub> into b<sub>3</sub> (cos θ, sin θ) into (-cos θ, -sin θ)</td></tr><tr><td>7</td><td>obtained from N̂<sub>1</sub>, N̂<sub>2</sub> by changing: (a<sub>1</sub>, a<sub>2</sub>, a<sub>3</sub>) into (a<sub>4</sub>, a<sub>1</sub>, a<sub>2</sub>)</td></tr><tr><td>8</td><td>b<sub>1</sub> into b<sub>4</sub> (cos θ, sin θ) into (-sin θ, cos θ)</td></tr></tbody></table>

where:

$$a_1 = 1 - \xi_{01}$$

$$b_1 = -1 - \xi_{01} - \xi_{02}$$

$$a_2 = 1 - \xi_{02}$$

$$b_2 = -1 + \xi_{01} - \xi_{02}$$

$$a_3 = 1 + \xi_{01}$$

$$b_3 = -1 + \xi_{01} + \xi_{02}$$

$$a_4 = 1 + \xi_{02}$$

$$b_4 = -1 - \xi_{01} + \xi_{02}$$

<table><caption>Table 6A. Function M<sup>(a)</sup></caption><thead><tr><th></th><th>M<sup>(1)</sup>(ξ<sub>0</sub>, λ, θ)</th><th>M<sup>(2)</sup>(ξ<sub>0</sub>, λ, θ)</th></tr></thead><tbody><tr><td>1</td><td>[(1 - ξ<sub>02</sub>)(2 cos θ + sin θ) - sin θ(2ξ<sub>1</sub> + ξ<sub>2</sub>)]</td><td>[(1 - ξ<sub>01</sub>)(2 sin θ + cos θ) - cos θ(2ξ<sub>2</sub> + ξ<sub>1</sub>)]</td></tr><tr><td>2</td><td>ξ<sub>01</sub> sin θ - cos θ(1 - ξ<sub>2</sub>)</td><td>cos θ(ξ<sub>01</sub> + λ cos θ/2)</td></tr><tr><td>3</td><td>[(1 - ξ<sub>02</sub>)(2 cos θ - sin θ) - sin θ(2ξ<sub>1</sub> - ξ<sub>2</sub>)]</td><td>[(1 + ξ<sub>01</sub>)(2 sin θ - cos θ) + cos θ(2ξ<sub>2</sub> - ξ<sub>1</sub>)]</td></tr><tr><td>4</td><td>-sin θ(ξ<sub>02</sub> + λ sin θ/2)</td><td>-(ξ<sub>02</sub> cos θ + sin θ(1 + ξ<sub>1</sub>))</td></tr><tr><td>5</td><td>[(1 + ξ<sub>02</sub>)(2 cos θ + sin θ) + sin θ(2ξ<sub>1</sub> + ξ<sub>2</sub>)]</td><td>[(1 + ξ<sub>01</sub>)(cos θ + 2 sin θ) + cos θ(2ξ<sub>2</sub> + ξ<sub>1</sub>)]</td></tr><tr><td>6</td><td>-(ξ<sub>01</sub> sin θ + cos θ(1 + ξ<sub>2</sub>))</td><td>-M<sup>(2)</sup></td></tr><tr><td>7</td><td>[(1 + ξ<sub>02</sub>)(2 cos θ - sin θ) + sin θ(2ξ<sub>1</sub> - ξ<sub>2</sub>)]</td><td>[(1 - ξ<sub>01</sub>)(2 sin θ - cos θ) - cos θ(2ξ<sub>2</sub> - ξ<sub>1</sub>)]</td></tr><tr><td>8</td><td>-M<sup>(1)</sup></td><td>ξ<sub>02</sub> cos θ - sin θ(1 - ξ<sub>1</sub>)</td></tr></tbody></table>

(Received 22 May 1989)