Local Approximation Study of Johnson’s Method for 2D Elasticity Problems

V.G. Ganzha(*) and E.V. Vorozhtsov

Institute of Theoretical and Applied Mechanics, Russian Academy of Sciences, Novosibirsk 630090, Russia

(*) At the moment the University of Kassel, Kassel D 34127, Germany

Abstract

The Lagrangian finite-difference method, which was developed

by G.R. Johnson in 1976 for the numerical solution of two-dimensional high velocity impact problems within the framework of the elasticity and plasticity theories, has gained a widespread acceptance. We investigate the local approximation of this method with the aid of the computer algebra systems. The problem under consideration involves tremendous intermediate expressions. We present an algorithm of decomposition of a complex approximation problem into a finite set of simpler subproblems whose union provides the solution of the original approximation problem. The results obtained with the aid of this algorithm show that the Johnson’s method has no approximation in the case when the accelerations of the nodes are computed on the basis of four triangular elements surrounding a given node of the triangular mesh. We further show that the first-order approximation in the case of six irregular elements surrounding the given mesh node takes place only in the case of small deviations of the shapes of the elements from their regular shapes. The validity of the obtained results was checked by using two different computer algebra systems: REDUCE and Maple. It follows from the results obtained that the Johnson’s method as applied to high-velocity impact problems should be augmented by an algorithm for the local mesh rezoning in the subdomains of strong distortions of the Lagrangian mass elements to ensure the local approximation.

1 Introduction

The experimental and theoretical analyses of high velocity impact problems have been of interest for many years. The mathematical investigation of these problems had begun with the advent of sufficiently powerful computers. The high-velocity impact problems are characterized by the presence of large strains and displacements. In addition, it is necessary to take into account material strength and compressibility effects. This leads to the need in the development of big computer codes, which enable one to obtain the approximate numerical solutions of the class of problems under consideration.

The methods based on Lagrangian variables and explicit finite-difference schemes have gained a widespread acceptance for the numerical solution of problems on the elastoplastic deformation of solid bodies [1-5]. An extension of the Johnson’s method has been presented in [4] for the cases of one, two and three spatial variables. This extension consists of the fact that the Lagrangian elements can be, for example, in the two-dimensional case not only triangular, but also quadrangular or pentagonal, so that in the same problem the Lagrangian elements are used, whose vertices form an arbitrary set of points.

The complexity of the computational formulas of the Johnson’s method does not enable one to study the local approximation of this method with the aid of analytic calculations by hand. In this connection we propose in what follows to use the computer algebra means for this purpose.

It should be noted that the symbolic computations have previously been used for the determination of the local approximation order of difference operators on irregular grids. In particular, a symbolic algorithm has been presented in [6] for the determination of the approximation order of difference operators on a one-dimensional non-uniform mesh. The method of [6] was extended in [7] for the case of the study of difference operators on two-dimensional non-orthogonal spatial grids. The system of analytic computations REDUCE [8] was used in [7]. The method of the work [7] was extended for the case of any number of the spatial dimensions in [9]. According to [9], this method can be applied also for the computer-aided investigation of the local approximation of difference operators on a mesh of Dirichlet cells.

The method of [7] enables one to solve the following problem: let a difference operator P_{h}u be given, which approximates certain operator P u on a spatial grid; for example, P u = grad(u). Find the value k in the formula

P_{h}u = P u + O(h^{k}) + o(h^{k}),

where O(h^{k}) is the leading term of the approximation error, and h is a quantity, which characterizes the local dimensions of the grid cells in the neighborhood of the spatial point under consideration. Thus the authors of [7, 9] have not considered the problem on the computation of an explicit form of the leading term of approximation error. However, it is well known that information of such kind is important for understanding the structure of the approximation ”viscosity” of the numerical method [10,11], since this structure affects the computational stability, the presence of oscillations, dispersion.

We describe below a symbolic algorithm, which enables us to investigate the approximation of the Johnson’s method [2] for the elastic regime of deformations

[figure: m; j; k; l; i; 4; 1; 2; 3; x; y]

Figure 1: Arrangement of four triangles into a composite element.

and to obtain an explicit form of the leading term of the approximation error on the triangular spatial grid.

2 Description of Johnson’s Method

In the Johnson’s method the solid bodies are represented by finite sets of the triangular cells. The displacements of the vertices of these cells in the passage from the time level t_{n} to the time level t_{n+1} = t_{n} + τ_{n} are computed with the aid of the approximation of the elasticity equations on the triangular mesh. Consider a triangular element ijk (see Fig. 1). Let u and v be the given displacements along the x-axis and y-axis, respectively. Denote the values of the displacements at the vertices of the element ijk by (u_{i}, v_{i}), (u_{j}, v_{j}), etc. We further denote the coordinates (x, y) of the element ijk vertices at t = t_{n} by (x_{i}^{0}, y_{i}^{0}), (x_{j}^{0}, y_{j}^{0}), etc.

We must calculate the strains within each Lagrangian element. Since the strains are expressed in terms of the derivatives of the displacements u and v with respect to spatial variables, we must at first find the form of the displacements as the functions of x and y within the limits of each element. It is assumed in the Johnson’s method that the displacements vary linearly between the nodes (the element vertices), that is

u = α_{1} + α_{2}x + α_{3}y,  (1)
v = α_{4} + α_{5}x + α_{6}y,  (2)

where the constants α_{1}, …, α_{6} depend on the geometry and the displacements of the nodes. The constants α_{1}, α_{2}, α_{3} are found as a result of the solution of the algebraic system

\begin{matrix}α_{1} + α_{2}x_{i}^{0} + α_{3}y_{i}^{0} = u_{i}, \\ α_{1} + α_{2}x_{j}^{0} + α_{3}y_{j}^{0} = u_{j}, \\ α_{1} + α_{2}x_{k}^{0} + α_{3}y_{k}^{0} = u_{k}. \end{matrix}  (3)

The determinant of this system is

Δ = | \begin{matrix}1 & x_{i}^{0} & y_{i}^{0} \\ 1 & x_{j}^{0} & y_{j}^{0} \\ 1 & x_{k}^{0} & y_{k}^{0} \end{matrix} | ,  (4)

The constants α_{4}, α_{5} and α_{6} are determined in a similar way; for this purpose the x-displacements are replaced in (3) by y-displacements.

It follows from (1)-(2) that the derivatives ∂u/∂x, ∂u/∂y, ∂v/∂x and ∂v/∂y are constant within the Lagrangian element. For example, ∂u/∂x = α_{2}. This enables us to determine the strains from the following equations [12]:

ε_{x} = \frac{∂u}{∂x} + \frac{1}{2}[(\frac{∂u}{∂x})^{2} + (\frac{∂v}{∂x})^{2}]  (5)

ε_{y} = \frac{∂v}{∂y} + \frac{1}{2}[(\frac{∂u}{∂y})^{2} + (\frac{∂v}{∂y} )^{2}]  (6)
γ_{xy} = \frac{∂u}{∂y} + \frac{∂v}{∂x} + \frac{∂u}{∂x} \frac{∂u}{∂y} + \frac{∂v}{∂x} \frac{∂v}{∂y}  (7)
ε_{v} = \frac{A}{A^{0}} − 1,  (8)

where A^{0} is the area of the element,

A^{0} = \frac{1}{2}[x_{i}^{0}(y_{j}^{0} − y_{k}^{0}) + x_{j}^{0}(y_{k}^{0} − y_{i}^{0}) + x_{k}^{0}(y_{i}^{0} − y_{j}^{0})].  (9)

The area A is computed by a similar formula on the basis of the displaced nodes coordinates (x_{i}, y_{i}), (x_{j}, y_{j}), (x_{k}, y_{k}). The normal strains are represented by Eqs. (5)-(6), whereas the shear strain is represented by Eq. (7). The volumetric strain ε_{v} is determined by Eq. (8). The higher-order terms are included in (5)-(7) to represent the effect of element rotation.

It should be noted that Eqs. (5)-(7) involve only the ”differential” terms, whereas Eq. (8) involves the finite differences. In the local approximation study usually the truncated Taylor series are substituted into the equations of the numerical method instead of the dependent variables. As a result one finds the expression of the form

”original PDE + the leading term of the approximation error”.

All the terms in this expression involve certain partial derivatives of the dependent variables. The expression for ε_{v} is purely algebraic, therefore, it cannot be used directly for the check-up of the results of the approximation study of the Johnson’s method. On the other hand, it is well known that the volumetric strain ε_{v} may be expressed in terms of the normal strains as follows (see, for example, [13, p. 70]):

ε_{v} = \frac{A}{A^{0}} − 1 = ε_{x} + ε_{y}.  (10)

The expression on the right hand side already involves only the differential terms in accordance with Eqs. (5) and (6).

As regards the computer implementation of the Johnson’s method the formula (8) is more preferable than Eq. (10), because the difference approximations of the normal strains ε_{x} and ε_{y} would involve many difference approximations of the derivatives ∂u/∂x, ∂v/∂x, …, ∂v/∂y in accordance with Eqs. (5),(6). It is well known that the approximation of the derivatives by the divided differences always leads to errors. On the other hand, the formula (8) indeed involves the difference approximation of only one derivative, namely (τ_{n}/A)(∂A/∂t), where τ_{n} is the time step of the numerical method.

In the present work we consider only the regime of elastic deformation. The elastic stresses can be obtained directly from the strains by use of Hooke’s law [12].

σ_{x} = λε_{v} + 2Gε_{x} − Q  (11)
σ_{y} = λε_{v} + 2Gε_{y} − Q  (12)
τ_{xy} = Gγ_{xy}  (13)

where σ_{x} and σ_{y} are the normal stresses, and τ_{xy} the shear stress. λ and G are Lamé elastic constants. The artificial viscosity Q, which is introduced for the smoothing of spurious oscillations of the numerical solution in the problems with propagating waves, includes both a linear term and a quadratic term and is applied only when the volumetric strain rate is negative.

\begin{matrix}Q = C_{L}\sqrt{(λ + 2G)ρA}|ε̇_{v}| + C_{0}^{2}ρA(ε̇_{v})^{2} for ε̇_{v} < 0 \\ Q = 0, ε̇_{v} ≥ 0, \end{matrix}  (14)

where the dot over ε_{v} denotes the differentiation with respect to t. C_{L} and C_{0}^{2} are the dimensionless constants the typical values of which are C_{L} = 0.5 and C_{0}^{2} = 4.0 [1,2].

After the element stresses are determined, the concentrated forces at the node i (cf. Fig. 1) are obtained. They are dependent on the coordinates x_{i}^{n+1}, y_{i}^{n+1}, etc., of the triangular elements vertices obtained after the displacement during the passage from the time t_{n} to the time t_{n+1}. The forces in the x and y-directions at node i of an element ”1” (see Fig. 1), are given by

F_{i,1}^{x} = (σ_{x})_{1}(y_{k} − y_{j}) + (τ_{xy})_{1}(x_{j} − x_{k});  (15)

F_{i,1}^{y} = (τ_{xy})_{1}(y_{k} − y_{j}) + (σ_{y})_{1}(x_{j} − x_{k}),  (16)

where (σ_{x})_{1} is the value of σ_{x} inside the element ijk, etc. The total force components

∑_{k=1}^{4} F_{i,k}^{x}, ∑_{k=1}^{4} F_{i,k}^{y}  (17)

in the node i of the element are obtained by summation of the contributions from each of the four elements having the vertex i (see Fig. 1).

The equations of motion are integrated numerically by assuming a constant velocity in the interval [t_{n}, t_{n} + τ_{n}]. The acceleration of the node i in the xdirection is

ü_{i} = (∑_{k=1}^{4} F_{i,k}^{x})/( ∑_{k=1}^{4} M_{i,k}).  (18)

The mass M_{i,k} in (18) is the contribution from the kth element to the lumped mass at the node i; this mass is computed in the Johnson’s method [2] by the formula

M_{i,k} = \frac{1}{3}A_{k}^{0}ρ_{k}^{0}, k = 1, 2, …, 4,  (19)

where ρ_{k}^{0} is the material density in the kth element at t = t_{n}, and only those four elements are taken, which have the common vertex i (see Fig. 1).

The new constant velocity in the x-direction is

u̇_{i}^{n+1} = u̇_{i}^{n} + ü_{i}τ_{n}.  (20)

Finally, the new displacement of the ith node at time t = t_{n+1} = t_{n} + τ_{n} is determined by the formula

u_{i}^{n+1} = u_{i}^{n} + u̇_{i}^{n+1}τ_{n}.  (21)

The equations of motion for the y-direction have a similar form.

We want to elucidate the following question: Do the computational formulas of the Johnson’s method approximate the fundamental equations of the elasticity theory, which have the following form in the notations given above:

ρü = \frac{∂σ_{x}}{∂x} + \frac{∂τ_{xy}}{∂y}  (22)

ρv̈ = \frac{∂τ_{xy}}{∂x} + \frac{∂σ_{y}}{∂y}  (23)

3 Algorithm of Decomposition for the Approximation Study

Our method of the local approximation study of the Johnson’s method is also based, as in [7,9], on the Taylor expansions of the dependent grid variables at a chosen spatial grid point. However, the approximation problem under consideration involves tremendous intermediate expressions in the case of an irregular mesh of triangular elements. We present below an algorithm of decomposition of a complex approximation problem into a finite set of simpler subproblems whose union provides the solution of the original approximation problem. Let us at first enumerate the basic steps of our decomposition algorithm.

Step 1. Representation of the differential operators on the right hand sides of Eqs. (22) and (23) as the sums of simpler differential operators.

Step 2. Identification of the difference operators approximating the differential operators identified at Step 1.

Step 3. Determination of the truncation error of each elementary difference operator identified at Step 2.

Step 4. Assembly of the results obtained at Step 3.

Let us now describe in detail each of the Steps 1–4.

Step 1. The difference approximations of the two equations (22), (23) for the nodes acceleration can be studied separately. The needed amount of computer memory thus reduces by a factor of two. In what follows we describe our decomposition strategy at the example of Eq. (22). Let u = (u, v). Then it is easy to see with regard for Eqs. (11)-(13) and (5)-(7) that

ρü = Lu,

where the differential operator L can be presented in accordance with (22), (5)-(7), (11)-(13) as the sum of the operators:

Lu = ∑_{μ=1}^{11} L_{μ}u,  (24)

where

L_{1}u = (λ + 2G)u_{xx}, L_{2}u = (λ + 2G)(\frac{1}{2}u_{x}^{2})_{x}, L_{3}u = (λ + 2G)(\frac{1}{2}v_{x}^{2})_{x},

L_{4}u = λv_{yx}, L_{5}u = λ(\frac{1}{2}u_{y}^{2})_{x}, L_{6}u = λ(\frac{1}{2}v_{y}^{2})_{x}, L_{7}u = Gu_{yy}, (25)
L_{8}u = Gv_{xy}, L_{9}u = G(u_{x}u_{y})_{y}, L_{10}u = G(v_{x}v_{y})_{y}, L_{11}u = −Q_{x}.

Step 2. The structure of the differencing formulas of the Johnson’s method enables us to easily identify the difference operators Λ_{μ}u approximating the differential operators L_{μ}u, μ = 1, …, 11.

Step 3. We now present the sequence of steps needed for the local approximation study of each elementary difference operator Λ_{μ}u, μ = 1, …, 11.

Step 3.1. The derivatives (∂u/∂x)_{ν}, (∂u/∂y)_{ν}, (∂v/∂x)_{ν}, ∂v/∂y)_{ν}, ν = 1, 2, 3, 4, which are constant within the νth element, are given the temporary names: (∂u/∂x)_{1} = bu21, (∂v/∂x)_{1} =bv21, (∂u/∂y)_{1} =bu31, (∂v/∂y)_{1} =bv31, etc.

Step 3.2. These temporary representations of the derivatives u_{x}, u_{y}, v_{x}, v_{y} for all four elements are substituted into the expression of the rational form of the type (18) for the difference operator Λ_{μ}u. As a result we obtain a preliminary expression for this operator.

Step 3.3. The algebraic system (3) is solved in a symbolic form for each of the four elements having their vertices in the node i, see Fig. 1. One can make the following table establishing a correspondence between the number ν of the element, 1 ≤ ν ≤ 4, and the indices of the vertices of the νth element:

\begin{matrix}ν & | & i, & i_{1}, & i_{2} \\ – & – & – & – & – \\ 1 & | & i, & j, & k \\ 2 & | & i, & k, & l \\ 3 & | & i, & l, & m \\ 4 & | & i, & m, & j \end{matrix}  (26)

As a result of the symbolic form solution of system (3) we obtain the rational expressions for the α_{1}, α_{2}, α_{3}, which involve the grid values u_{i}, u_{j}, u_{k} of the displacements.

Step 3.4. The temporary names buξη, bvξη, ξ = 2, 3, η = 1, 2, 3, 4, in the expression for Λ_{μ}u are replaced by the corresponding expressions for α_{2}, α_{3}, α_{5}, α_{6} obtained at Step 3.3. This leads to a considerable increase of the length of the expression for Λ_{μ}u.

Step 3.5. The displacements u_{ν} and v_{ν} entering the obtained expression for Λ_{μ}u are expanded in accordance with the Taylor formula:

u_{ν} = u^{∗} + \frac{∂u(x^{∗}, y^{∗})}{∂x}(x_{ν} − x^{∗}) + \frac{∂u(x^{∗}, y^{∗})}{∂y}(y_{ν} − y^{∗}) +

+\frac{1}{2} \frac{∂^{2}u(x^{∗}, y^{∗})}{∂x^{2}}(x_{ν} − x^{∗})^{2} + \frac{∂^{2}u(x^{∗}, y^{∗})}{∂x∂y}(x_{ν} − x^{∗})(y_{ν} − y^{∗}) +

\frac{1}{2} \frac{∂^{2}u(x^{∗}, y^{∗})}{∂y^{2}}(y_{ν} − y^{∗})^{2},  (27)

where ν = i, j, k, l, m and (x^{∗}, y^{∗}) is a point with respect to which the expansions are performed.

Step 3.6. The specification of the (x^{∗}, y^{∗}) point. The computational results presented below have been obtained at x^{∗} = x_{i}^{0}, y^{∗} = y_{i}^{0}.

Step 3.7. Representation of the initial coordinates of the four nodes (x_{j}^{0}, y_{j}^{0}), …, (x_{m}^{0}, y_{m}^{0}) as the deviations from the point (x^{∗} = x_{i}^{0}, , y^{∗} = y_{i}^{0}):

x_{j}^{0} = x^{∗} + C_{j}^{x} · h, x_{k}^{0} = x^{∗} + C_{k}^{x} · h, x_{l}^{0} = x^{∗} + C_{l}^{x} · h,
x_{m}^{0} = x^{∗} + C_{m}^{x} · h, y_{j}^{0} = y^{∗} + C_{j}^{y} · h, y_{k}^{0} = y^{∗} + C_{k}^{y} · h,  (28)
y_{l}^{0} = y^{∗} + C_{l}^{y} · h, y_{m}^{0} = y^{∗} + C_{m}^{y} · h.

Here h is a quantitative measure, which locally characterizes the density of a grid of triangular elements in the neighborhood of the (x^{∗}, y^{∗}) point. Following [9] let us denote by S(i) the stencil of the ith node, that is the vertices of those five elements, which are involved in the computation of the expression Λ_{μ}u. In the notations of Fig. 1 these are the vertices i, j, k, l, m. If we take x^{∗} = x_{i}^{0}, y^{∗} = y_{i}^{0}, then h may be computed by the formula

h = (1/4){ ∑_{ν∈S(i)}[(x_{ν} − x^{∗})^{2} + (y_{ν} − y^{∗})^{2}]}^{0.5},

where S(i) = {j, k, l, m}. The constants C_{j}^{x}, …, C_{m}^{y} in (28) do not depend on h. The representations (28) characterize the arbitrary distortions, which the Lagrangian mass elements can undergo in the process of high-velocity interaction of solid bodies.

Step 3.8. The displaced coordinates (x_{ν}, y_{ν}), ν = i, j, k, l, m, are also involved in the computation of the difference operator Λ_{μ}u. In the Johnson’s method these coordinates are related to the coordinates (x_{ν}^{0}, y_{ν}^{0}) by the formulas

x_{ν} = x_{ν}^{0} + u_{ν}, y_{ν} = y_{ν}^{0} + v_{ν}, ν = i, j, k, l, m.  (29)

If the displacements u(x, y) and v(x, y) are smooth functions in the vicinity of the node (x_{i}^{0}, y_{i}^{0}), then we can write that

u(x_{ν}^{0}, y_{ν}^{0}) = u(x_{i}^{0}, y_{i}^{0}) + \frac{∂u(x̄, ȳ)}{∂x}(x_{ν}^{0} − x_{i}^{0}) + \frac{∂u(x̄, ȳ)}{∂y}(y_{ν}^{0} − y_{i}^{0}), (30)

where

x̄ = x_{i}^{0} + θ_{1}(x_{ν}^{0} − x_{i}^{0}), ȳ = y_{i}^{0} + θ_{2}(y_{ν}^{0} − y_{i}^{0}),

0 < θ_{1}, θ_{2} < 1. Substituting the representations (28) into (30) we can obtain the formula

u_{ν} = u_{i} + A_{ν}^{x}h, ν = i, j, k, l, m,  (31)

where

A_{ν}^{x} = \frac{∂u(x̄, ȳ)}{∂x}C_{ν}^{x} + \frac{∂u(x̄, ȳ)}{∂y}C_{ν}^{y}.

Aluminum plate

Figure 2: Uniform triangular mesh. The case of two diagonals in each rectangular cell.

In a similar way we can obtain the representations

v_{ν} = v_{i} + A_{ν}^{y}h, ν = i, j, k, l, m.  (32)

Substituting the formulas (31) and (32) into (29) we obtain that

x_{ν} = x_{ν}^{0} + u_{i} + A_{ν}^{x}h, y_{ν} = y_{ν}^{0} + v_{i} + A_{ν}^{y}h, ν = i, j, k, l, m.  (33)

The representations (33) are substituted into the expression for the operator Λ_{μ}u.

At Step 4, the assembly of the expansions of the difference operators Λ_{μ}u, μ = 1, …, 11, is performed in accordance with (24) simply by the concatenation of the files containing the expansion results for the individual operators Λ_{μ}u.

4 Results and Discussion

Let us consider similarly to [2] a problem on the high-velocity impact of a nickel rectangular projectile onto a thick aluminum plate. Let us introduce an uniform rectangular grid in the both of the bodies. Let h_{1} and h_{2} be the sizes of the cells in the directions of the x- and y-axes, respectively. Let us denote by κ the cell aspect ratio, i.e. κ = h_{2}/h_{1}. Now let us draw two diagonals in each of the rectangular cells. Then we obtain an uniform triangular mesh shown in Fig. 2.

Let us take an individual rectangular cell (see Fig. 3). We can see that four elements have a common vertex at node i. Let us consider equation (22) for the

[figure: m; j; k; l; i; 4; 1; 2; 3]

Figure 3: Individual rectangular cell and its triangulation.

x-acceleration. In what follows we present the results of the local approximation study of the sum

Λ_{l}u = Λ_{1}u + Λ_{4}u + Λ_{7}u + Λ_{8}u  (34)

of the linear operators present in the complete differential operator (24) corresponding to Eq. (22).

We have implemented the algorithm of decomposition presented above in Section 3 in a REDUCE program. To verify the results obtained by the RE-DUCE program we have written a similar program in Maple [14]. The symbolic results obtained by the two programs for the case of the investigation of the operator (34) are the same and have the following form:

Λ_{l}u = \frac{3}{2}[(λ + 2G)u_{xx} + 2λv_{xy} + 2κ^{2}Gu_{yy} +

κ^{2}λu_{yy} + Gu_{yy} + 2Gv_{xy} + \frac{1}{κ^{2}} Gu_{xx}] + O(h).  (35)

The exact result should however have in accordance with Eqs. (25) the form

L_{l}u = L_{1}u + L_{4}u + L_{7}u + L_{8}u =
(λ + 2G)u_{xx} + λv_{yx} + Gu_{yy} + Gv_{xy}.  (36)

Comparing formulas (35) and (36) we can see that there is no approximation of the operator L_{l}u by the Johnson’s method. To understand the reason for this we have considered only the operator Λ_{1}u. It turns out that the Johnson’s method approximates instead of the operator (λ + 2G)u_{xx} the operator

\frac{3}{2}(λ + 2G)(u_{xx} + κ^{2}u_{yy}).  (37)

This fact can in turn be explained as follows. In the particular case of a rectangular cell shown in Fig. 3 we must specify the following values of the constants

m  jk l  i 4 2  1 3

Figure 4: The partition of a rectangular cell into four cells with the stepsizes h_{1}/2 and κh_{1}/2.

C_{ν}^{x} and C_{ν}^{y} in Eqs. (28):

\begin{matrix}C_{j}^{x} = \frac{1}{2}, & C_{j}^{y} = −\frac{κ}{2}, & C_{k}^{x} = \frac{1}{2}, & C_{k}^{y} = \frac{κ}{2}, \\ C_{l}^{x} = −\frac{1}{2}, & C_{l}^{y} = \frac{κ}{2}, & C_{m}^{x} = −\frac{1}{2}, & C_{m}^{y} = −\frac{κ}{2}. \end{matrix}  (38)

Then we can obtain that the difference operator Λ_{1}u has the form

Λ_{1}u = (λ + 2G)\frac{(u_{j} + u_{k})κ − 4u_{i} + u_{m} + u_{l}}{(1/3)h_{1}^{2}} .  (39)

Let us now subdivide the elementary rectangular cell of Fig. 3 into four rectangular cells shown in Fig. 4. Then we can see from Eq. (38) that on the mesh with the steps h_{1}/2 and h_{2}/2 the difference operator Λ_{1}u can be represented schematically as the discrete convolution

Λ_{1}u = H_{1} ∗ F_{3},  (40)

where

H_{1} = \frac{3(λ + 2G)}{h_{1}^{2}} [ \begin{matrix}1 & 0 & κ \\ 0 & −4 & 0 \\ 1 & 0 & κ \end{matrix} ] , F_{3} = [ \begin{matrix}u_{l} & a & u_{k} \\ b & u_{i} & c \\ u_{m} & d & u_{j} \end{matrix} ] .  (41)

We have denoted by a, b, c, d in (41) the missing values of u at points 1,2, 3,4; these values are indeed not needed, because they are multiplied by zero, see the expression for the mask H_{1} in (41). On the other hand, it is well known that in the particular case of a square grid with the stepsize h the Laplacian Δu can be approximated with the aid of equation [15]

Δu = H_{2} ∗ F_{3},

Aluminum plate

Figure 5: Uniform triangular mesh. The case of one diagonal in each rectangular cell.

where

H_{2} = \frac{1}{2h̃^{2}} [ \begin{matrix}1 & 0 & 1 \\ 0 & −4 & 0 \\ 1 & 0 & 1 \end{matrix} ] .

In our case h̃ = h_{1}/2, so that

H_{2} = \frac{2}{h^{2}} [ \begin{matrix}1 & 0 & 1 \\ 0 & −4 & 0 \\ 1 & 0 & 1 \end{matrix} ] .  (42)

Comparing the expressions for the masks H_{1} in (41) and H_{2} (42) we can see that (i) the elements of the masks coincide in the particular case of κ = 1; (ii) the factor ”2” is absent in the numerator of the coefficient of H_{1}. This explains the appearance of the ”2” in the denominator of the formula (37).

Let us now consider the case when only one diagonal is drawn in each rectangular cell of an uniform grid. To be specific let us consider the case of a diagonal drawn at the angle α = arctan(h_{2}/h_{1}) to the positive direction of the x-axis. Then we obtain a triangular mesh shown in Fig. 5. It is easy to see from Fig. 5 that in this case six triangular elements have a common vertex at point i. It is not difficult to write down the difference formulas similar to Eqs. (15), (16) and (18) for this case. We will not present these formulas for the sake of brevity. We present only an analog of the formulas (28). Let us denote the vertices of the six elements surrounding the node i as in Fig. 6. Then we can write that

x_{ν}^{0} = x_{i}^{0} + C_{ν}^{x}h, y_{ν}^{0} = y_{i}^{0} + C_{ν}^{y}h, ν = j, k, l, m, n.  (43)

[figure: n; m; l; k; i; j; p; x; y]

Figure 6: The environment of the node i.

In the particular case of an uniform rectangular grid we should substitute the values

\begin{matrix}C_{j}^{x} = 1, & C_{j}^{y} = 0, & C_{k}^{x} = 1, & C_{k}^{y} = κ, \\ C_{l}^{x} = 0 & C_{l}^{y} = κ, & C_{m}^{x} = −1, & C_{m}^{y} = 0, \\ C_{n}^{x} = −1, & C_{n}^{y} = −κ \end{matrix}  (44)

into the formulas (43). Then the approximation of the linear operator (36) by the Johnson’s method takes place only when the multiplier 2/3 is taken in the formula (18) instead of the multiplier 1/3 presented in [2]. In the case of the multiplier 2/3 we have obtained the second order of approximation in space, that is the order O(h^{2}), where h = h_{1}.

Let us finally consider the approximation of the linear operator L_{1}u = (λ + 2G)u_{xx} by the corresponding operator Λ_{1}u in the case of the six distorted triangular elements surrounding the node i. Such severe distortions of the triangular mass elements occur in the region of the high-velocity interaction of a projectile with an obstacle. Fig. 7 gives some idea of these distortions. In this case we should use the general formulas (28) in the approximation study. The expression for the Taylor expansion of the operator Λ_{1}u is sufficiently large, it takes about 800 Kb. Therefore, we present only the general structure of this expression together with a few monomials involved:

Λ_{1}u = 3h(λ + 2G)(A_{1}u_{xx} + 2A_{2}u_{xy} + A_{3}u_{yy})/(2hB),  (45)

where

A_{1} = (C_{j}^{x})^{3}C_{k}^{x}C_{l}^{x}C_{m}^{x}C_{n}^{x}(C_{k}^{y})^{2}C_{l}^{y}C_{m}^{y}C_{n}^{y}(C_{p}^{y})^{2} + …
−C_{k}^{x}C_{l}^{x}C_{m}^{x}C_{n}^{x}(C_{p}^{x})^{3}(C_{j}^{y})^{2}C_{k}^{y}C_{l}^{y}C_{m}^{y}(C_{n}^{y})^{2},

i

Figure 7: Distorted six elements surrounding the node i.

A_{2} = (C_{j}^{x})^{2}C_{k}^{x}C_{l}^{x}C_{m}^{x}C_{n}^{x}C_{j}^{y}(C_{k}^{y})^{2}C_{l}^{y}C_{m}^{y}C_{n}^{y}(C_{p}^{y})^{2} + …

−C_{k}^{x}C_{l}^{x}C_{m}^{x}C_{n}^{x}(C_{p}^{x})^{2}(C_{j}^{y})^{2}C_{k}^{y}C_{l}^{y}C_{m}^{y}(C_{n}^{y})^{2}C_{p}^{y},

A_{3} = −(C_{j}^{x})^{2}C_{k}^{x}C_{l}^{x}C_{m}^{x}C_{k}^{y}C_{l}^{y}C_{m}^{y}(C_{n}^{y})^{4}(C_{p}^{y})^{2} + …

+C_{l}^{x}C_{m}^{x}C_{n}^{x}(C_{p}^{x})^{2}(C_{j}^{y})^{2}(C_{k}^{y})^{4}C_{l}^{y}C_{m}^{y}C_{n}^{y},
B = (C_{j}^{x})^{3}C_{k}^{x}C_{l}^{x}C_{m}^{x}C_{n}^{x}(C_{k}^{y})^{2}C_{l}^{y}C_{m}^{y}C_{n}^{y}(C_{p}^{y})^{2} +
… − C_{k}^{x}C_{l}^{x}C_{m}^{x}C_{n}^{x}(C_{p}^{x})^{3}(C_{j}^{y})^{2}C_{k}^{y}C_{l}^{y}C_{m}^{y}(C_{n}^{y})^{2}.  (46)

In each of the expressions (45) we have presented only the first and the last monomial. Note that the terms of similar structures were collected in these expressions by the REDUCE system. The expressions (45), (46) are presented under the assumption that the correct multiplier (2/3) is used instead of (1/3) in Eq. (18).

It follows from Eq. (45) that there is no approximation of the operator (λ + 2G)u_{xx} by the Johnson’s method in the case of large distortions of the triangular cells.

Let us now consider another case, when the shapes of the six triangular elements surrounding the node i differ little from the regular shape shown in Fig. 6, in the sense that the relationships (cf. (44))

\begin{matrix}C_{j}^{x} = 1 + D_{j}^{x}h, & C_{j}^{y} = D_{j}^{y}h, & C_{k}^{x} = 1 + D_{k}^{x}h, \\ C_{k}^{y} = κ + D_{k}^{y}h, & C_{l}^{x} = D_{l}^{x}h, & C_{l}^{y} = κ + D_{l}^{y}h, \\ C_{m}^{x} = −1 + D_{m}^{x}h, & C_{m}^{y} = D_{m}^{y}h, & C_{n}^{x} = −1 + D_{n}^{x}h, \\ C_{n}^{y} = −κ + D_{n}^{y}h, & C_{p}^{x} = D_{p}^{x}h, & C_{p}^{y} = −κ + D_{p}^{y}h \end{matrix}  (47)

are satisfied, where the constants D_{ν}^{x}, D_{ν}^{y} have the order of smallness O(1) and do not depend on h. In this case the operator Λ_{1}u approximates the operator L_{1}u with the order O(h):

Λ_{1}u = L_{1}u + O(h).

Substituting the formulas (46) into Eqs. (28), we can see that, for example,

x_{j}^{0} = x_{i}^{0} + h + D_{j}^{x}h^{2}.

This means that the positions of the vertices of the triangular elements should differ by an amount O(h^{2}) from their regular positions to ensure a local approximation.

From this investigation we can draw a conclusion that it is reasonable to perform a local mesh rezoning in the result of which the regions of large distortions of the triangular cells are replaced by the regions of ”nearly regular cells”.

5 Conclusions

The obtained results on the local approximation of the Johnson’s method can be summarized as follows.

1. There is no approximation in the case when the accelerations of the nodes are computed on the basis of four triangular elements surrounding a given node of the triangular mesh.

2. In the case of six elements and of the uniform triangular mesh, the approximation takes place only if the lumped mass in the given node is equal to two thirds of the sum of masses of the six elements having a common vertex at the given node. In this particular case, the Johnson’s method has the second order of approximation in space.

3. In the case of irregular, distorted Lagrangian mass elements the Johnson’s method has no approximation of the original elasticity equations also in the case of six triangular elements surrounding the given node of the mesh.

4. The first-order approximation in the case of six irregular elements takes place only in the case of small deviations of the shapes of the elements from their regular shapes.

5. It is reasonable to augment the Johnson’s method by an algorithm for the local mesh rezoning in the subdomains of severe distortions of the Lagrangian mass elements to ensure the local approximation.

In concluding we wish to formulate the following problem: Is it possible to modify the Johnson’s method in such a way that it possesses the local approximation also in the case of severe distortions of the triangular elements?

A possible way for the solution of this problem is to refuse the linear approximations (1), (2) inside the elements and to use the finite volume method and the Green’s theorem for the approximation of the second derivatives u_{xx}, u_{xy},…, v_{yy} as this was described, for example, in [16].

Acknowledgments

We express our gratitude to Prof. V.M. Fomin, the director of the Institute of Theoretical and Applied Mechanics of the Siberian Division of the Russian Academy of Sciences, Novosibirsk, who has drawn our attention to the problem of the approximation study of the Johnson’s method.

References

1. M.L. Wilkins. Calculation of elastic-plastic flow. In: Methods in Computational Physics, Vol. 3 / Eds. B. Alder, S. Fernbach and M. Rotenberg. Academic Press, New York, 1964, p. 211-263.

2. G.R. Johnson. Analysis of elastic-plastic impact involving severe distortions. J. Appl. Mech., vol. 43, No. 3, Trans. ASME, Vol. 98, Series E, Sept. 1976, p. 439-444.

3. G.R. Johnson. Liquid-solid impact calculations with triangular elements. Trans. ASME, ser. I, vol. 99, No. 3, p. 598-600 (1977).

4. A.I. Gulidov. Organization of the computational process and the data structures at the numerical solution of problems of the mechanics of deformable media. Modeling in Mechanics (in Russian), 1991, vol. 5, No. 3, p. 127-141.

5. V.V. Bashurov and S.K. Buruchenko. Modeling of high-velocity penetration in Lagrangian coordinates. Mathematical Modeling (in Russian), 1992, vol. 4, No. 9, p. 37-42.

6. V.G. Ganzha, M.Yu. Shashkov and V.P. Shapeev. Local approximation study of difference operators in symbolic form. Preprint of the Keldysh Institute of Applied Mathematics of the USSR Academy of Sciences (in Russian), No. 79, Moscow, 1987, p. 1-12.

7. V.G. Ganzha and M.Yu. Shashkov. Local approximation study of difference operators by means of REDUCE system. In: ISSAC’90 Proceedings /Eds. S. Watanabe and M. Nagata. Addison-Wesley Publishing Company, Reading, 1990, p. 185-192.

8. A.C. Hearn. REDUCE 3.5, User’s Manual. Rand Corporation, Santa Monica, Cal., 1993.

9. V.G. Ganzha, A.V. Solovyov, and M.Yu. Shashkov. A new algorithm for the symbolic computation of the analyses of local approximation of difference operators. In: Advances in Computer Methods for Partial Differential Equations – VII. Proceedings of the Seventh IMACS Internat. Conf. on Computer Methods for Partial Differential Equations, New Brunswick, New Jersey, USA, June 22-24, 1992. / Eds. R.Vichnevetsky, D. Knight and G.Richter. IMACS Publ., New Brunswick, NJ 08903, 1992, p. 280-286.

10. Y.I. Shokin. The Method of Differential Approximation. (Translated from the Russian by K.G. Roesner). Springer-Verlag, Berlin, Heidelberg, New York, 1983.

11. Yu.I. Shokin and N.N. Yanenko. The Method of Differential Approximation. Application to Gas Dynamics. (in Russian). Nauka, Siberian Division, Novosibirsk, 1985.

12. A.P. Boresi. Elasticity in Engineering Mechanics. Prentice-Hall, Inc., Englewood Cliffs, N.J., 1965.

13. B.L. Glushak, V.F. Kuropatenko, and S.A. Novikov. Investigation of Materials Strength at Dynamic Loadings. (in Russian). Nauka, Novosibirsk, 1992.

14. B.W. Char, K.O. Geddes, G.H. Gonnet, M.B. Monagan, and S.M. Watt. Maple Reference Manual, 5th Edition. Waterloo Maple Publishing, Waterloo, Ontario, Canada, 1990.

15. P.J. Roache. Computational Fluid Dynamics. Hermosa, Albuquerque, NM, 1976.

16. C.A.J. Fletcher. Computational Techniques for Fluid Dynamics 1. Fundamental and General Techniques. Springer-Verlag, Berlin, Heidelberg, New York, 1988.