R Chapko - On the numerical solution of a cauchy problem for an elastostatic equation - страница 1

1  2  3 



. . .

Ser. Appl. Math. Inform.

2009. Bun. 15. C. 135-148

2009. Is. 15. P. 135-148

УДК 519.6


R. Chapko, O. Sobeyko

Ivan Franko National University of L 'viv Universytetska st., 1, Lviv, 79000, e-mail: chapko@is.lviv.ua

In this paper a numerical solution of a Cauchy problem for an elastostatic equation in a double-connected domain in R2 is considered. It is shown how an iterative alternating procedure in conjunction with indirect boundary integral equations method could be applied in order to efficiently find a stable solution of the problem. Different techniques of dealing with singularities in kernels of the integral operators and numerical approximation are included.

Key words: elastostatics, Navier equation, integral equations method, alternating procedure, Cauchy problem, quadrature method, numerical solution.


In many fields of scientific research, be it economics, finances, physical simulation with non-destructive exploration, so-called inverse problems arise. Such problems are of great mathematical interest, but they are usually hard to solve, mostly because of their ill-posedness. In case of the Cauchy problem for an elastostatic equation the stability of the solution suffers. In order to obtain a stable numerical solution an appropriate method should be applied. Usually numerical methods for such kind of problems have regularizing properties such that they provide a stable approximation to the solution. An iterative alternating procedure [3, 5] proposed by Kozlov and Mazya is suitable for this purpose. This procedure has already been applied in case of the Laplace equation (see for example [3]) and has shown stable results. Each iteration of the alternating procedure requires solving of two mixed direct and well-posed boundary value problems for the elastostatic equation.

It is also very important not only to solve the problem but to do this efficiently as well. That is why an indirect integral equations method based on the elastic potentials theory [1] is used.


Let D <z R2 be a double-connected domain with smooth boundaries Г0 and Г. The boundaries have no common points and do not intersect themselves. An example of the domain is shown on Fig. 1.

Physically it may be assumed that the domain is filled with a solid isotropic elastic material. The contour Г could be treated as an accessible boundary of the body and the contour Г0 - as an inaccessible one.

© Chapko R., Sobeyko O., 2009

x: л

Fig. 1. Example of the domain

Consider the Cauchy boundary value problem for the elastostatic equation [1]. The statement is as follows: let a vector-function u є C2(D) n C1(D) satisfy the elastostatic equation (or the Navier equation)

jAu + + jj) graddivu = 0   in   D, (1) the Dirichlet boundary condition

u = gon Г (2) and the Neumann boundary condition

Tu = f on Г . (3)

Find values of u and Tu on the inaccessible boundary Г0. We may also consider a wider problem of reconstructing u and Tu inside of the domain and particularly on Г0.

The constants j and Я ( j > 0, Я>-/ ) are known as Lame coefficients which

characterize physical properties of the material. The operator T is a stress operator and it is defined by

Tu = Яdivuv + 2ji{v- grad )u + j div(Qu )Qv, (4) where v - is an outer ( with respect to the domain ) unit normal vector to the boundary and


Together the functions from the boundary conditions (2) and (3) are known as the Cauchy data.

From the physical point of view the vector-function u defines the displacements of the solid in any point, including the boundary points. Tu defines the stress which occurs on the boundary of the solid and which corresponds to the displacements.

This inverse problem is ill-posed since its solution is unstable. This means that a small perturbation of the input data could lead to a great perturbation of the solution.

Q is a matrix given by Q

The functions g and are known boundary

Therefore a suitable regularizing method should be applied to gain the approximate solution's stability. One of the most appropriate methods is a so-called alternating method proposed by Kozlov and Mazya [5].


In order to formulate an alternating iterative procedure for the Cauchy problem two mixed boundary value problems for the Navier equation have to be introduced.

The first problem is known as the Dirichlet-Neumann mixed boundary value

problem and it is formulated as follows: find a function u є C1 (D) n C (D) that solves the

elastostatic equation (1), satisfies the Dirichlet boundary condition

u = g і on T0 (5)

and the Neumann boundary condition

Tu = f1 on Г . (6) The second boundary value problem is known as the Neumann-Dirichlet mixed boundary  value  problem.   It  has  the  following  formulation:   find  a function

u є C2(D) n C1(D) that solves the elastostatic equation (1), satisfies the Neumann boundary condition

Tu = g2 on Го (7) and the Dirichlet boundary condition

u = f 2 on Г . (8)

Now based on the mixed boundary value problems an alternating procedure could be formulated as

the first approximation u0 to the solution u of (1)-(3) on Г0 could be chosen as an arbitrary initial guess;

having constructed the k-th approximation uk the mixed Dirichlet-Neumann problem with the boundary conditions g1 = u, L and f = f should be solved in order to find approximation to Tuk on Г0;

having found Tuk the mixed Neumann-Dirichlet problem with the boundary conditions g2 = Tuk and f2 = g should be solved in order to find the next iterative approximation uk+1.

As it could be noticed the alternating method requires solving of two mixed direct boundary value problems on each iteration, which are well-posed. This is the cause why we get a stable approximation of the searched solution on every iteration. However the more iterations are being performed the less stable the solution becomes. This means that increasing the number of iterations leads to loss of the regularizing properties of the algorithm. As a stopping criterion of the iterative process the Morozov's discrepancy principle may be applied. However in some cases this criterion causes the algorithm to stop too early.


As it has already been mentioned, each iteration of the alternating method assumes solving two direct mixed boundary value problems. Technically the problems could be solved by using any of the well known numerical methods. But not all of them are suitable from the computational point of view. It should be emphasized that we are searching the solution of the Cauchy problem particularly on the interior boundary Г0. That is why

applying such method as FEM is possible but not preferable as it might lead to a significant increase of required computational resources. Therefore, a boundary integral equations method is applied to solve the Dirichlet-Neumann and the Neumann-Dirichlet direct problems. This approach is based on the potential theory in case of elastostatics. The fundamental solution for the Navier equation is given by

Ф (x, y) = A. W (x, y)I + A J (x - y), (9) 2n 2n



Я + 3ji Л + jU

+ 2 ц)        ц(Я + 2ц)

W (x, y) = ln,-- x Ф y.

x - y

The matrix I is an identity 2 by 2 matrix and J is defined as


J(h) :=т,со* 0.

In order to apply the integral equations method an integral representation of the searched solution is needed. Using the potential theory approach we may seek the solution in the form of a single-layer elastic potential

u(x) = JФ(x, yp (y)ds(y) + ]Ф(x, y)<p(y)ds(y),   x є D. (10)

Here p0 and p are some vector densities. The integral representation (10) of the

solution is perfectly suitable for both the Dirichlet-Neumann and the Neumann-Dirichlet cases. One of the main advantages of such approach is that it allows us to decrease the number of dimensions of the problem by one as all of the integrals in the solutions representation are contour integrals. The boundary curves Г0 and Г are assumed to have parametric representation:

Г0:= M(') = WW))|   0 <' < 2п} (11)


Г := W) = WIM'))|   0 <' < 2n}. (12)



In case of the Dirichlet-Neumann boundary value problem we introduce the following integral operators:

(Sy)(x) := JФ(x,y)y/(y)ds(y),   xёГ0,

У)( x) := ]Ф( x, y)ys( y)ds( y),   x є Г о


(N¥)( x):= J (Tx Ф)( x, y)¥(y)ds(y)


(H¥)( x):= J(Tx Ф)( x, y)¥( y)ds( y)

x є Г, x є Г.


(14) (15) (16)

Taking into account the integral representation (10) and the well known boundary properties of the elastic potentials (see [1]) the problem (1)-(5)-(6) could be reduced to the system of integral operator equations of the form:

(0%)(x) + (C<p)(x) = gl(x),   x єГо,

[(Np)(x) + ((/ + H)p)(x) = f1(x),   x є Г. () Here the functions <ро and p are unknown and have to be computed.

It should be noticed that the integral operators S and H contain singularities in their kernels. Thus S should be understood as an improper integral and H - in sense of a Cauchy principal value.

Well-posedness of the obtained integral equations system in corresponding Holder spaces is established by the following theoretical result [2].

Theorem 1. For m єгЧ, m < p , where p describes smoothness of the boundary curves Го and Г, for g1 є Cma (Го) and for f є Ст-1а(Г) with а є (0,1) the integral equations system (17) has exactly one solution po є Ст-1,а(Го) and рє Ст-1,а(Г) which continuously depends on the input data.

Proof. The system of integral equations (17) may be rewritten in the operator form

's    о W о с vo i+h){n о

In terms of the theorem and according to [6] the integral operators s and i + h may be shown bounded with bounded inverse operators s 1 and (I + h)_1 in corresponding spaces. Also the integral operators c and n are compact since their kernels are continuous. Therefore, the system may be rewritten as




\ ґ +


(I + H)-1N 0

S 1 g1 (I + H)-1 f


Since the direct Dirichlet-Neumann boundary value problem is uniquely solvable and the corresponding boundary value problem with homogeneous boundary conditions has only a

trivial solution we may now apply the main result of Riesz theory to the operator equation of the second kind with a compact operator. This finally proves the theorem.Having obtained the system of boundary integral equations and taking into consideration the parametric representation of the boundary curves the system is parameterized and the integral operators have the following parametric representation:


(SpXYW) = J Ф (\,(0,\,(т))Р,(т) \y\TT)\ dT, (18)



(CpXYXO) = J Ф(\0«,\(т)ЖТ)|\'(t)| dT, (19)



(NftXYO) = J (Tm ФШ)\0ТТ)Шт)\\\(t)\ dT, (20)



(HpW)) = J(Tr(t), \(т))р(т) |\'(t)| dT. (21)

1  2  3 

Похожие статьи

R Chapko - On the numerical solution of a cauchy problem for an elastostatic equation