Formulation and Numerical Solution of Plane Problems of the Theory of Elasticity in Strains

Abstract

This article is devoted to the formulation and numerical solution of boundary-value problems in the theory of elasticity with respect to deformations. Similar to the well-known Beltrami–Michell stress equations, the Saint-Venant compatibility conditions are written in the form of differential equations for strains. A new version of plane boundary-value problems in strains is formulated. It is shown that for the correctness of plane boundary value problems, in addition to the usual conditions, one more special boundary condition is required using the equilibrium equation. To discretize additional boundary conditions and differential equations, it is convenient to use the finite difference method. By resolving grid equations and additional boundary conditions with respect to the desired quantities at the diagonal nodal points, we obtained convergent iterative relations for the internal and boundary nodes. To solve grid equations, the elimination method was also used. By comparing with the Timoshenko–Goodyear solution on the tension of a rectangular plate with a parabolic load, the validity of the formulated boundary value problems in strains and the reliability of the numerical results are shown. The accuracy of the results has been increased by an average of 15%.

1. Introduction

The development of information technologies in the modern world, as well as their widespread application in various fields of scientific and technical applications, have set new goals and more complex innovative tasks for scientists and workers. They must adequately calculate the safety margins and reliability of structures and their elements, considering the influence of external factors. The mathematical and numerical modeling of linear and nonlinear processes of the deformation of engineering structures in mechanical engineering, nuclear power engineering, aircraft engineering, and astronautics; in the calculation of dams and the mining industry; as well as the study of the stress–strain state to determine the safety margins is an urgent problem in solid mechanics.

Typically, the boundary-value problems of the theory of elasticity are formulated with respect to displacements, and the necessary strains and stresses are calculated from the displacements.

The formulation of boundary value problems on stresses and strains is an urgent problem in solid mechanics. The formulation of boundary value problems is usually based on the conditions of compatibility of Saint-Venant deformations. From a mathematical point of view, the conditions for the compatibility of deformations are the result of the triviality of the components of the Riemann–Christophel tensor and provide an unambiguous determination of displacements along deformations for simply connected regions [1]. It is known that the compatibility conditions consist of six equations and the question of the dependence of the compatibility equations remains unresolved [2]. Note that to formulate plane boundary value problems in the theory of elasticity, one condition for the compatibility of deformations is sufficient. Plane problems usually reduce to solving a biharmonic equation with respect to the Airy stress function.

It is known that the conditions for compatibility of deformations, using Hooke’s law and the equilibrium equation, can be written with respect to the stress tensor in the form of the Beltrami–Michell equations [3]. The Beltrami–Michell equations, in combination with three equilibrium equations, represent a boundary value problem with nine equations and three boundary conditions [4]. The works of Borodachev [5] show that the first group of three Beltrami–Michell equations depends on the second group of equations. In the works of Pobedry [6], the compatibility conditions and equilibrium equations are reduced to a correct boundary value problem consisting of six equations [7]. In this case, the equilibrium equations on the boundary of a given region are considered as the three missing boundary conditions. In a particular case, the Beltrami–Michell equations follow from the Pobedry equations [8]. Issues of equivalence in the formulation of boundary value problems on displacements and stresses are considered in [9]. Questions of the existence and uniqueness of solutions to boundary value problems are considered in [10]. The Beltrami–Michell equations taking into account temperature are considered in the work of Nowatsky [8]. Coupled problems of thermoelasticity are considered in [11]. Dynamic boundary value problems in stresses are considered in the works of Konovalov [12].

The formulation of boundary value problems regarding deformations is a poorly studied area of solid mechanics. In this area, the works of Pobedra [3,4] and Borodachev [6,7] can be noted. In Pobedry’s works, the deformation compatibility equation, in combination with the equilibrium equation, is written in the form of six differential equations for the components of the deformation tensor. In [6,7], within the framework of the Beltrami–Michell equations, equations about deformations of an infinite half-plane are considered. Despite the existing effective methods for solving applied problems, such as the finite element method, FEM, and finite difference methods, there are few numerically solved boundary value problems regarding stresses. Let us note the works of Filonenko-Borodich [13]. The problem of equilibrium of a parallelepiped under stress was considered by the variational-difference method in [3,14,15].

In the studies of [16,17], the regular perturbation method is employed to solve the fundamental equations of fluid flow, encompassing continuity, momentum, mode, and energy. This approach models the distribution of velocity, mode, and temperature. The work’s novelty stems from treating the particle penetration speed as a perturbation parameter. Crucially, this rate of particle penetration at the base of the flow channel, which depends on the porosity of the gas diffusion layer (GDL) and the operational pressure within the channel, significantly influences the performance metrics, specifically the output voltage at a given current density, of proton exchange membrane fuel cells (PEMFCs).

This study is concerned with the formulation and numerical solution of boundary-value problems of the theory of elasticity in strains. Within the framework of the compatibility conditions, differential equations of deformations are expressed, which, in combination with the equilibrium equations and the corresponding boundary and additional boundary conditions, constitute the boundary problem of the theory of elasticity in strains. It is shown that the first two differential equations of deformation in the plane strain case are equivalent to the well-known condition ${\epsilon }_{11}{}_{,22}+{\epsilon }_{22}{}_{,11}=2{\epsilon }_{12}{}_{,12}$, and the third equation can also be considered as a new compatibility condition. Moreover, these two compatibility conditions combined with two equilibrium equations allow us to constitute two different plane strain boundary problems. In addition to the usual boundary conditions, additional equations are required for the correct formulation of boundary-value problems, which are obtained by considering the equilibrium equation on the domain’s boundary. It is shown that the differentiated equilibrium equations, in conjunction with the compatibility condition, can also be considered a boundary-value problem in strains.

Grid equations were compiled using the finite difference method for plane boundary-value problems in the strains, which were solved using the iterative and variable direction methods. By comparing the numerical results of the plane boundary-value problems with the well-known Timoshenk-Goodyear [18] solution for stretching a rectangular plane with a parabolic load applied to opposite sides, the validity of the formulated boundary-value problems of the theory of elasticity in strains and the reliability of the results obtained were substantiated [19,20].

2. Formulation of the Boundary-Value Problems of the Theory of Elasticity in Strains

Generally, Refs. [6,8] the boundary-value problem of elasticity theory consists of the following equilibrium equation:

$${\sigma }_{ij,j}+X_i=0,$$

(1)

Hooke’s law is expressed as follows:

$${\sigma }_i{}_j=\lambda \theta {\delta }_{ij}+2\mu {\epsilon }_{ij},$$

(2)

Cauchy ratio is expressed as follows:

$${\epsilon }_{ij}=\frac{1}{2}(u_{i,j}+u_{j,i}),$$

(3)

Boundary conditions are expressed as follows:

$$u_i{|}_{{\sum }_1}=u_i^0,$$

(4)

$${\sigma }_{ij}{n}_j{|}_{{\sum }_2}=S_i.$$

(5)

where ${\sigma }_{ij}$—denotes the stress tensor, ${\epsilon }_{ij}$—denotes the strain tensor, $u_i$—denotes the displacement, $\mathsf{\lambda },\mathsf{\mu }$ denotes the elastic Lame constants, ${\mathrm{e}}_{\mathrm{ij}},\theta$—denotes spherical part of the strain tensor, $S_i$—denotes the surface load, ${\mathrm{X}}_{\mathrm{i}}$ denotes the body forces, and ${\mathsf{\delta }}_{\mathrm{ij}}$ denotes the Kronecker symbol.

Substituting Equation (3) into Equation (2) from Equation (1), we can obtain the following differential equations for displacements in the form of the Lame equation:

$$\mu {\nabla }^2{u}_i+(\lambda +\mu )\theta {,}_i+X_i=0.$$

(6)

where the ${\nabla }^2$—denotes the Laplace operator and $\theta ={\epsilon }_{kk}$. If necessary, deformations and stresses can be calculated from the displacements.

The condition for the unique solvability of Equation (3) with respect to the displacements is the Saint-Venant compatibility condition, which is expressed as follows:

$${\epsilon }_{ij,kl}+{\epsilon }_{kl,ij}-{\epsilon }_{ik,jl}-{\epsilon }_{jl,ik}=0,$$

(7)

Multiplying the latter by ${\delta }_{kl}$ can be reduced to the following six equations:

$${\nabla }^2{\epsilon }_{ij}+{\theta }_{,ij}-{\epsilon }_{ik,jk}-{\epsilon }_{jk,ik}=0,$$

(8)

Using Hooke’s law (2) from the equilibrium Equation (1), one can receive

$${\epsilon }_{ij,j}=-\frac{1}{2\mu }(\lambda \theta {,}_i+X_i),$$

and substitute the last one into Equation (8) to find the differential equations in strain [8]:

$$\mu {\nabla }^2{\epsilon }_{ij}+(\lambda +\mu )\theta {,}_i{}_j+\frac{1}{2}(X_{i,j}+X_{j,i})=0,$$

(9)

Setting up a correct boundary value problem using Equation (9) is an unexplored complex mathematical problem. There are several reasons that prevent the correct formulation of the boundary value problem:

-

The boundary value problem consisting of Equation (9) with boundary conditions (4) and (5) does not describe the process of deformation of the solid bodies under study;

-

To formulate a correct boundary value problem, Equation (9) must be considered in combination with the equilibrium equation; then, the number of equations becomes equal to nine and the problem of choosing three independent equations from six (9) arises;

-

Boundary conditions (5) consist of three conditions, and for the correct formulation of the boundary value problem, three more boundary conditions will be required;

-

Equilibrium equations can be considered missing boundary conditions, but their numerical implementation is still unclear.

Equation (9), taking into account the equilibrium Equation (1) and boundary conditions (5), can be expressed in terms of deformations

$$\lambda \theta {,}_i{\delta }_{ij}+2\mu {\epsilon }_{ij,j}+X_i=0,$$

(10)

$${\left.(\lambda \theta {\delta }_{ij}+2\mu {\epsilon }_{ij})n_j\right|}_{{\Sigma }_2}=S_i,$$

(11)

and following additional boundary conditions

$$(\lambda \theta {,}_i{\delta }_{ij}+2\mu {\epsilon }_{ij,j}+X_i)\left|{}_{\Sigma }\right.=0$$

(12)

that represent the boundary-value problem of the theory of elasticity in strains [15].

Note that in the boundary-value problems of the theory of elasticity in strains, the boundary conditions do not depend on the derivatives of the desired quantities; that is, they are fulfilled exactly and, therefore, do not contain errors of numerical differentiation, in contrast to the boundary-value problems solved for displacements. The boundary-value problem (9)–(12) is discussed in the two-dimensional case in the next section.

3. Classical Plane Problems of Elasticity Theory in Stresses and Strains

Prior to the discussion of the boundary-value problem in strains (9)–(12), we first consider the typical plane problem of elasticity theory. In the absence of body forces, it consists of two equilibrium equations as follows:

$$\begin{array}{l}\frac{\partial {\sigma }_{11}}{\partial x} + \frac{\partial {\sigma }_{12}}{\partial y} =0,\\ \frac{\partial {\sigma }_{21}}{\partial x} + \frac{\partial {\sigma }_{22}}{\partial y} =0,\end{array}$$

(13)

It also consists of strain compatibility conditions [21]:

$$\frac{{\partial }^2{\epsilon }_{11}}{\partial y^2}+\frac{{\partial }^2{\epsilon }_{22}}{\partial x^2}=2\frac{{\partial }^2{\epsilon }_{12}}{\partial x\partial y}.$$

(14)

Hooke’s law for plane problems assumes the following form [8,18]:

$${\epsilon }_{11}=\frac{1}{E_1}{\sigma }_{11}-\frac{{\nu }_1}{E}{\sigma }_{22}, {\epsilon }_{22}=\frac{1}{E_1}{\sigma }_{22}-\frac{{\nu }_1}{E}{\sigma }_{11}, {\epsilon }_{12}=\frac{1}{2\mu }{\sigma }_{12}$$

(15)

where

$$E_1=\left\{\begin{array}{l}\frac{E}{1-{\nu }^2} plane strain state\\ E \hspace{0.17em}plane stress state \end{array}\right. {\nu }_1=\left\{\begin{array}{l}\frac{\nu }{1-\nu } p.s.s\\ \nu \hspace{0.17em}p.s.s\end{array}\right.$$

It is known that the strain compatibility condition (14), with the assistance of the equilibrium Equation (13) and Hooke’s law (15), can be expressed as a harmonic equation as follows [22]:

$${\nabla }^2({\sigma }_{11}+{\sigma }_{22})=0$$

(16)

In this case, the boundary conditions have the following form:

$$\begin{array}{l}({\sigma }_{11}{n}_1+{\sigma }_{12}{n}_2){|}_{\Gamma }=S_1,\\ ({\sigma }_{21}{n}_1+{\sigma }_{22}{n}_2){|}_{\Gamma }=S_2,\end{array}$$

(17)

Equations (13), (16), and (17) represent the classical plane problem of the theory of elasticity in stress (Problem A). Problem A consists of three equations for the components of the stress tensor ${\sigma }_{11}, {\sigma }_{22}, {\sigma }_{12}$, with two boundary conditions. Problem A was typically reduced to solving a biharmonic equation with respect to the Airy stress function [18].

The classical plane problem can also be formulated with respect to the strain. To achieve this, using Hooke’s law (2), we express the equilibrium Equation (13) in terms of deformations, which, together with the compatibility condition (14), constitute a plane problem of the theory of elasticity in terms of strains (Problem B),

$$\begin{array}{l}(\lambda +2\mu )\frac{\partial {\epsilon }_{11}}{\partial x}+\lambda \frac{\partial {\epsilon }_{22}}{\partial x}+2\mu \frac{\partial {\epsilon }_{12}}{\partial y}=0,\\ (\lambda +2\mu )\frac{\partial {\epsilon }_{22}}{\partial y}+\lambda \frac{\partial {\epsilon }_{11}}{\partial y}+2\mu \frac{{\partial }^2{\epsilon }_{12}}{\partial x}=0,\\ \hspace{0.17em}\frac{{\partial }^2{\epsilon }_{11}}{\partial y^2}+\frac{{\partial }^2{\epsilon }_{22}}{\partial x^2}=2\frac{{\partial }^2{\epsilon }_{12}}{\partial x\partial y}\end{array}$$

(18)

with appropriate boundary conditions,

$$\begin{array}{l}({\sigma }_{11}{n}_1+{\sigma }_{12}{n}_2){|}_{\Gamma }=S_1,\\ ({\sigma }_{21}{n}_1+{\sigma }_{22}{n}_2){|}_{\Gamma }=S_2,\end{array}$$

(19)

where

$$\begin{array}{l}{\sigma }_{11 }=(\lambda +2\mu ){\epsilon }_{11}+\lambda {\epsilon }_{22},\\ {\sigma }_{22 }=\lambda {\epsilon }_{11}+(\lambda +2\mu ){\epsilon }_{22}, {\sigma }_{12}=2\mu {\epsilon }_{12}\end{array}$$

The boundary-value problem B also consists of three equations for the strain tensor components, ${\epsilon }_{11}, {\epsilon }_{22}, {\epsilon }_{12}$, and two boundary conditions.

In boundary-value problems A and B, an additional boundary condition is required for the correct formulation of the boundary-value problems. In the case of Problem A, the absence of a boundary condition is compensated for by introducing a stress function that identically satisfies the equilibrium equation, and the problem is reduced to solving a biharmonic equation with respect to the Airy stress function [10,21].

In the case of Problems B and A, the missing boundary condition, following [23,24,25], can be determined by considering the equilibrium equation on the boundary of the given domain.

4. New Plane Problems of the Theory of Elasticity in Strains

This section discusses the plane boundary-value problems of elasticity theory based on the boundary-value problem (9)–(12) (in the absence of body forces), that is,

$$(\lambda +2\mu )\frac{{\partial }^2{\epsilon }_{11}}{\partial x^2}+\mu \frac{{\partial }^2{\epsilon }_{11}}{\partial y^2}+(\lambda +\mu )\frac{{\partial }^2{\epsilon }_{22}}{\partial x^2}=0,$$

(20)

$$(\lambda +2\mu )\frac{{\partial }^2{\epsilon }_{22}}{\partial y^2}+\mu \frac{{\partial }^2{\epsilon }_{22}}{\partial x^2}+(\lambda +\mu )\frac{{\partial }^2{\epsilon }_{11}}{\partial y^2}=0,$$

(21)

$$\mu (\frac{{\partial }^2{\epsilon }_{12}}{\partial x^2}+\frac{{\partial }^2{\epsilon }_{12}}{\partial y^2})+(\lambda +\mu )(\frac{{\partial }^2{\epsilon }_{11}}{\partial x\partial y}+\frac{{\partial }^2{\epsilon }_{22}}{\partial x\partial y})=0,$$

(22)

$$(\lambda +2\mu )\frac{\partial {\epsilon }_{11}}{\partial x}+\lambda \frac{\partial {\epsilon }_{22}}{\partial x}+2\mu \frac{\partial {\epsilon }_{12}}{\partial y}=0,$$

(23)

$$(\lambda +2\mu )\frac{\partial {\epsilon }_{22}}{\partial y}+\lambda \frac{\partial {\epsilon }_{11}}{\partial y}+2\mu \frac{{\partial }^2{\epsilon }_{12}}{\partial x}=0.$$

(24)

The differential Equations (20)–(22) are a consequence of the deformation compatibility condition (7), and they can also, in principle, be considered the Saint-Venant compatibility conditions. However, the Saint-Venant compatibility condition in the plane case consists of one well-known condition (14). It appears that there must be some connection between these Equations (20)–(22) and (14).

To clarify this issue, we differentiate equilibrium Equations (23) and (24) with respect to x and y, respectively [26]:

$$(\lambda +2\mu )\frac{{\partial }^2{\epsilon }_{11}}{\partial x^2}+\lambda \frac{{\partial }^2{\epsilon }_{22}}{\partial x^2}=-2\mu \frac{{\partial }^2{\epsilon }_{12}}{\partial x\partial y},$$

(25)

$$(\lambda +2\mu )\frac{{\partial }^2{\epsilon }_{22}}{\partial y^2}+\lambda \frac{{\partial }^2{\epsilon }_{11}}{\partial y^2}=-2\mu \frac{{\partial }^2{\epsilon }_{12}}{\partial x\partial y},$$

(26)

Equations (20) and (21), considering relations (25) and (26), can be reduced, respectively, to the following form:

$$-2\mu \frac{{\partial }^2{\epsilon }_{12}}{\partial x\partial y}+\mu \frac{{\partial }^2{\epsilon }_{11}}{\partial y^2}+\mu \frac{{\partial }^2{\epsilon }_{22}}{\partial x^2}=0,$$

(27)

$$-2\mu \frac{{\partial }^2{\epsilon }_{12}}{\partial x\partial y}+\mu \frac{{\partial }^2{\epsilon }_{22}}{\partial x^2}+\mu \frac{{\partial }^2{\epsilon }_{11}}{\partial y^2}=0.$$

(28)

By adding these equations, we can obtain the well-known Saint-Venant compatibility condition (14) as follows:

$$\frac{{\partial }^2{\epsilon }_{11}}{\partial y^2}+\frac{{\partial }^2{\epsilon }_{22}}{\partial x^2}=2\frac{{\partial }^2{\epsilon }_{12}}{\partial x\partial y}.$$

It can be observed that the first two Equations (20) and (21) are equivalent to the compatibility condition (14). From this, we can conclude that Equation (22) can also be used as a compatibility condition instead of Equation (14).

Thus, the differential Equations (20)–(24) can be divided into two plane problems consisting of two equations in combination, in the first case with Equation (14), which is equivalent to Equations (20) and (21), and in the second case with Equation (24).

The first case, from the equilibrium Equations (23) and (24) in combination with (14), follows the plane boundary-value problem B, formulated in the previous paragraph by Equation (18):

$$\begin{array}{l}(\lambda +2\mu )\frac{\partial {\epsilon }_{11}}{\partial x}+\lambda \frac{\partial {\epsilon }_{22}}{\partial x}+2\mu \frac{\partial {\epsilon }_{12}}{\partial y}=0,\\ (\lambda +2\mu )\frac{\partial {\epsilon }_{22}}{\partial y}+\lambda \frac{\partial {\epsilon }_{11}}{\partial y}+2\mu \frac{{\partial }^2{\epsilon }_{12}}{\partial x}=0,\\ \frac{{\partial }^2{\epsilon }_{11}}{\partial y^2}+\frac{{\partial }^2{\epsilon }_{22}}{\partial x^2}=2\frac{{\partial }^2{\epsilon }_{12}}{\partial x\partial y}.\end{array}$$

The equilibrium Equations (23) and (24) are considered together with the third differential Equation (22):

$$\begin{array}{l}(\lambda +2\mu )\frac{\partial {\epsilon }_{11}}{\partial x}+\lambda \frac{\partial {\epsilon }_{22}}{\partial x}+2\mu \frac{\partial {\epsilon }_{12}}{\partial y}=0,\\ (\lambda +2\mu )\frac{\partial {\epsilon }_{22}}{\partial y}+\lambda \frac{\partial {\epsilon }_{11}}{\partial y}+2\mu \frac{{\partial }^2{\epsilon }_{12}}{\partial x}=0,\\ \mu (\frac{{\partial }^2{\epsilon }_{12}}{\partial x^2}+\frac{{\partial }^2{\epsilon }_{12}}{\partial y^2})+(\lambda +\mu )(\frac{{\partial }^2{\epsilon }_{11}}{\partial x\partial y}+\frac{{\partial }^2{\epsilon }_{22}}{\partial x\partial y})=0,\end{array}$$

(29)

which constitutes another plane strain problem (Problem C).

When deriving the differential Equation (9), the equilibrium equations are used in a differentiated form. Therefore, in boundary-value problem C, the use of differentiated equilibrium equations is beyond doubt. Thus, we obtain an additional boundary-value problem for the theory of elasticity in strain (Problem D):

$$\begin{array}{l}(\lambda +2\mu )\frac{{\partial }^2{\epsilon }_{11}}{\partial x^2}+\lambda \frac{{\partial }^2{\epsilon }_{22}}{\partial x^2}+2\mu \frac{{\partial }^2{\epsilon }_{12}}{\partial x\partial y}=0,\\ (\lambda +2\mu )\frac{{\partial }^2{\epsilon }_{22}}{\partial y^2}+\lambda \frac{{\partial }^2{\epsilon }_{11}}{\partial y^2}+2\mu \frac{{\partial }^2{\epsilon }_{21}}{\partial x\partial y}=0,\\ \mu (\frac{{\partial }^2{\epsilon }_{12}}{\partial x^2}+\frac{{\partial }^2{\epsilon }_{12}}{\partial y^2})+(\lambda +\mu )(\frac{{\partial }^2{\epsilon }_{11}}{\partial x\partial y}+\frac{{\partial }^2{\epsilon }_{22}}{\partial x\partial y})=0\end{array}$$

(30)

For boundary-value problems B, C, and D, the boundary conditions (19) have the following form:

$$\begin{array}{l}((\lambda +2\mu ){\epsilon }_{11}+\lambda {\epsilon }_{22})n_1+2\mu {\epsilon }_{12}{n}_2){|}_{\Gamma }=S_1,\\ (2\mu {\epsilon }_{12}{n}_1+(\lambda {\epsilon }_{11}+(\lambda +2\mu ){\epsilon }_{22})n_2){|}_{\Gamma }=S_2,\end{array}$$

(31)

with additional boundary conditions (12)

$$\begin{array}{l}\left[(\lambda +2\mu )\frac{\partial {\epsilon }_{11}}{\partial x}+\lambda \frac{\partial {\epsilon }_{22}}{\partial x}+2\mu \frac{\partial {\epsilon }_{12}}{\partial y}\right]{|}_{\Gamma }=0,\\ \left[(\lambda +2\mu )\frac{\partial {\epsilon }_{22}}{\partial y}+\lambda \frac{\partial {\epsilon }_{11}}{\partial y}+2\mu \frac{\partial {\epsilon }_{21}}{\partial x}\right]{|}_{\Gamma }=0,\end{array}$$

(32)

The boundary conditions (19) for a rectangular region (Figure 1) have the following form:

$$\begin{array}{l}for x=0,l_1: {\sigma }_{11}{|}_{x=0,l_1}=\phi , {\sigma }_{12}{|}_{x=0,l_1}=0,\\ for y=0,l_2: {\sigma }_{22}{|}_{y=0,l_2}=0, {\sigma }_{21}{|}_{y=0,l_2}=0.\end{array}$$

(33)

Considering Hooke’s law (15), the boundary conditions for deformation can be expressed in the following form:

$$\begin{array}{l}{\epsilon }_{22}{|}_{y=0}=\frac{1}{E_1}{\sigma }_{22}, {\epsilon }_{12}{|}_{y=0}=0, {\epsilon }_{22}{|}_{y=l_2}=-\frac{1}{E_1}{\sigma }_{22}, {\epsilon }_{12}{|}_{y=l_2}=0,\\ {\epsilon }_{11}{|}_{x=0}=0, {\epsilon }_{21}{|}_{x=0}=0, {\epsilon }_{11}{|}_{x=l_1}=0, {\epsilon }_{21}{|}_{x=l_1}=0.\end{array}$$

(34)

The additional boundary conditions (12) for a rectangular area can be obtained from (32) at y = 0, l₂ and x = 0, l₁ (Figure 1) for ${\epsilon }_{11}$ and ${\epsilon }_{22}$, respectively.

$$\begin{array}{l}\left[\frac{\partial {\epsilon }_{11}}{\partial y}\right]{|}_{y=0,l_2}=-\left[\frac{2\mu }{\lambda }\frac{\partial {\epsilon }_{21}}{\partial x}\right]{|}_{y=0,l_2},\\ \left[\frac{\partial {\epsilon }_{22}}{\partial x}\right]{|}_{x=0,l_1}=-\left[\frac{2\mu }{\lambda }\frac{\partial {\epsilon }_{12}}{\partial y}\right]{|}_{x=0,l_1} \end{array}$$

(35)

Boundary conditions (34) and (35) are valid for boundary-value problems $B, C$, and $D.$

5. Finite-Difference Equations of Plane Problems of the Theory of Elasticity in Strains and Methods for Their Solution

This section describes the construction of numerical models for the plane problems $B, C$, and $D$ considered in Section 3 and Section 4 and a comparison of their numerical results.

Let us consider the boundary-value problem B in the rectangular region of $\Omega =\{0\le x\le l_1 , 0\le y\le l_2\}$. To construct a finite-difference scheme, dividing the length of the sides $l_k$ of a rectangle by $N_k$, it can be observed that $h_k=l_k/N_k$, where $k=1,2$. Subsequently, the nodal points have the following form:

$$x_i=h_1\cdot i\hspace{0.33em},\hspace{1em}i=\overline{0,N_1}, y_j=h_2\cdot j\hspace{0.33em},\hspace{1em}j=\overline{0,N_2}$$

By replacing the derivatives with the corresponding finite difference relations, the difference equations for Problem B are obtained:

$$(\lambda +2\mu )\frac{{\epsilon }_{i+1,j}^{11}-{\epsilon }_{ij}^{11}}{2h_1}+\lambda \frac{{\epsilon }_{i+1,j}^{22}-{\epsilon }_{i-1,j}^{22}}{2h_1}+2\mu \frac{{\epsilon }_{i,j+1}^{12}-{\epsilon }_{i,j-1}^{12}}{2h_2}=0,$$

(36)

$$(\lambda +2\mu )\frac{{\epsilon }_{ij+1}^{22}-{\epsilon }_{i,j}^{22}}{2h_2}+\lambda \frac{{\epsilon }_{i,j+1}^{11}-{\epsilon }_{i,j-1}^{11}}{2h_2}+2\mu \frac{{\epsilon }_{i+1,j}^{12}-{\epsilon }_{i-1,j}^{12}}{2h_1}=0,$$

(37)

$$\frac{{\epsilon }_{i,j+1}^{11}-2{\epsilon }_{i,j}^{11}+{\epsilon }_{i,j-1}^{11}}{h_2^2}+\frac{{\epsilon }_{i+1,j}^{22}-2{\epsilon }_{i,j}^{22}+{\epsilon }_{i-1,j}^{22}}{h_1^2}=2\frac{1}{h_2}(\frac{{\epsilon }_{i,j+1}^{12}-{\epsilon }_{i-1,j+1}^{12}}{h_1}-\frac{{\epsilon }_{i,j}^{12}-{\epsilon }_{i-1,j}^{12}}{h_1})$$

(38)

Resolving these Equations (36)–(38) with respect to ${\epsilon }_{ij}^{11}, {\epsilon }_{ij}^{22} , {\epsilon }_{ij}^{12}$, we obtain the following:

$$\begin{array}{l}{\epsilon }_{ij}^{11}={\epsilon }_{i+1,j}^{11}+\frac{2h_1}{\lambda +2\mu }(\lambda \frac{{\epsilon }_{i+1,j}^{22}-{\epsilon }_{i-1,j}^{22}}{2h_1}+2\mu \frac{{\epsilon }_{i,j+1}^{12}-{\epsilon }_{i,j-1}^{12}}{2h_2}),\\ {\epsilon }_{ij}^{22}={\epsilon }_{i,j+1}^{22}+\frac{2h_2}{\lambda +2\mu }(\lambda \frac{{\epsilon }_{i,j+1}^{11}-{\epsilon }_{i,j-1}^{11}}{2h_2}+2\mu \frac{{\epsilon }_{i+1,j}^{12}-{\epsilon }_{i-1,j}^{12}}{2h_1}),\\ {\epsilon }_{ij}^{12}={\epsilon }_{i-1,j}^{12}+h_1(\frac{{\epsilon }_{i,j+1}^{12}-{\epsilon }_{i-1,j+1}^{12}}{h_1}-\frac{h_2}{2}(\frac{{\epsilon }_{i,j+1}^{11}-2{\epsilon }_{i,j}^{11}+{\epsilon }_{i,j-1}^{11}}{h_2^2}+\frac{{\epsilon }_{i+1,j}^{22}-2{\epsilon }_{i,j}^{22}+{\epsilon }_{i-1,j}^{22}}{h_1^2})).\end{array}$$

(39)

Using the following schemes for Equation (36),

$$\begin{array}{l}{\epsilon }_{ij}^{11}={\epsilon }_{i+1,j}^{11}+\frac{2h}{\lambda +2{\mu }_1}(\lambda \frac{{\epsilon }_{i+1,j}^{22}-{\epsilon }_{i-1,j}^{22}}{2h_1}+2\mu \frac{{\epsilon }_{i,j+1}^{12}-{\epsilon }_{i,j-1}^{12}}{2h_2}),\\ {\epsilon }_{ij}^{11}={\epsilon }_{i-1,j}^{11}-\frac{2h}{\lambda +2{\mu }_1}(\lambda \frac{{\epsilon }_{i+1,j}^{22}-{\epsilon }_{i-1,j}^{22}}{2h_1}+2\mu \frac{{\epsilon }_{i,j+1}^{12}-{\epsilon }_{i,j-1}^{12}}{2h_2}),\end{array}$$

(40)

after adding these two equations, we obtain the following relation:

$${\epsilon }_{ij}^{11}=\frac{{\epsilon }_{i+1,j}^{11}+{\epsilon }_{i-1,j}^{11}}{2},$$

(41)

Similarly, from (37) for ${\epsilon }_{ij}^{22}$, we can observe the following:

$${\epsilon }_{ij}^{22}=\frac{{\epsilon }_{i,j+1}^{22}+{\epsilon }_{i,j-1}^{22}}{2},$$

(42)

Instead of the first two Equations (39), Equations (41) and (42) can be used.

Equations (36)–(38), in combination with the boundary conditions (34) and (35), constitute a difference analog of the boundary-value problem B. Solving these equations with respect to ${\epsilon }_{11}, {\epsilon }_{22}, {\epsilon }_{12}$, we can obtain relation (39), which makes it possible to obtain the desired values at internal points using the iterative method. The additional boundary conditions in (35) at the nodal points have the following form:

for y = 0 and y = l₂

$$\begin{array}{l}{\epsilon }_{i,0}^{11} = {\epsilon }_{i,1}^{11} +\frac{\mu h_2}{\lambda }\frac{{\epsilon }_{i+1,0}^{12}-{\epsilon }_{i-1,0}^{12}}{h_1},\\ {\epsilon }_{i,N_2}^{11}={\epsilon }_{i,N_2-1}^{11}-\frac{\mu h_2}{\lambda }\frac{{\epsilon }_{i+1,N_2}^{12}-{\epsilon }_{i-1,N_2}^{12}}{h_1},\end{array}$$

(43)

for x = 0 and x = l₁

$$\begin{array}{l}{\epsilon }_{0,j}^{22}={\epsilon }_{1,j}^{22}+\frac{\mu h_1}{\lambda }\frac{{\epsilon }_{0,j+1}^{12}-{\epsilon }_{0,j-1}^{12}}{h_2},\\ {\epsilon }_{N_1,j}^{22}={\epsilon }_{N_1-1,j}^{22}-\frac{\mu h_1}{\lambda }\frac{{\epsilon }_{N_1,j+1}^{12}-{\epsilon }_{N_1,j-1}^{12}}{h_2}.\end{array}$$

(44)

Considering the boundary and additional conditions, the difference analog of Problem C (29) has the following form:

$$(\lambda +2\mu )\frac{{\epsilon }_{i+1,j}^{11}-{\epsilon }_{ij}^{11}}{2h_1}+\lambda \frac{{\epsilon }_{i+1,j}^{22}-{\epsilon }_{i-1,j}^{22}}{2h_1}+2\mu \frac{{\epsilon }_{i,j+1}^{12}-{\epsilon }_{i,j-1}^{12}}{2h_2}=0,$$

(45)

$$(\lambda +2\mu )\frac{{\epsilon }_{ij+1}^{22}-{\epsilon }_{i,j}^{22}}{2h_2}+\lambda \frac{{\epsilon }_{i,j+1}^{11}-{\epsilon }_{i,j-1}^{11}}{2h_2}+2\mu \frac{{\epsilon }_{i+1,j}^{12}-{\epsilon }_{i-1,j}^{12}}{2h_1}=0,$$

(46)

$$\begin{array}{c}\mu (\frac{{\epsilon }_{i+1,j}^{12}-2{\epsilon }_{i,j}^{12}+{\epsilon }_{i,j}^{12}}{h_1^2}+\frac{{\epsilon }_{i,j+1}^{12}-2{\epsilon }_{i,j}^{12}+{\epsilon }_{i,j-1}^{12}}{h_2^2})+(\lambda +\mu )(\frac{{\epsilon }_{i+1,j+1}^{11}-{\epsilon }_{i+1,j-1}^{11}-{\epsilon }_{i-1,j+1}^{11}+{\epsilon }_{i-1,j-1}^{11}}{4h_1{h}_2}\\ +\frac{{\epsilon }_{i+1,j+1}^{22}-{\epsilon }_{i+1,j-1}^{22}-{\epsilon }_{i-1,j+1}^{22}+{\epsilon }_{i-1,j-1}^{22}}{4h_1{h}_2})=0.\end{array}$$

(47)

Solving these equations for ${\epsilon }_{11}, {\epsilon }_{22}, {\epsilon }_{12}$, similar to Problem B, we can obtain the following expressions solved by the iteration method:

$${\epsilon }_{ij}^{11}=\frac{{\epsilon }_{i+1,j}^{11}+{\epsilon }_{i-1,j}^{11}}{2},$$

(48)

$${\epsilon }_{ij}^{22}=\frac{{\epsilon }_{i,j+1}^{22}+{\epsilon }_{i,j-1}^{22}}{2},$$

(49)

$$\begin{array}{l}{\epsilon }_{i,j}^{12}=((\lambda +\mu )(\frac{{\epsilon }_{i+1,j+1}^{11}-{\epsilon }_{i+1,j-1}^{11}-{\epsilon }_{i-1,j+1}^{11}+{\epsilon }_{i-1,j-1}^{11}}{4h_1{h}_2}+\frac{{\epsilon }_{i+1,j+1}^{22}-{\epsilon }_{i+1,j-1}^{22}-{\epsilon }_{i-1,j+1}^{22}+{\epsilon }_{i-1,j-1}^{22}}{4h_1{h}_2})\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} +\mu (\frac{{\epsilon }_{i+1,j}^{12}+{\epsilon }_{i-1,j}^{12}}{h_1^2}+\frac{{\epsilon }_{i,j+1}^{12}+{\epsilon }_{i,j-1}^{12}}{h_2^2}))/(\frac{2\mu }{h_1^2}+\frac{2\mu }{h_2^2}).\end{array}$$

(50)

We now discuss the solution to Problem D (30). The finite-difference analog of Problem D (30) has the following form:

$$(\lambda +2\mu )\frac{{\epsilon }_{i+1,j}^{11}-2{\epsilon }_{i,j}^{11}+{\epsilon }_{i-1,j}^{11}}{h_1^2}+\lambda \frac{{\epsilon }_{i+1,j}^{22}-2{\epsilon }_{i,j}^{22}+{\epsilon }_{i-1,j}^{22}}{h_1^2}+2\mu \frac{{\epsilon }_{i+1,j+1}^{12}-{\epsilon }_{i-1,j+1}^{12}-{\epsilon }_{i+1,j-1}^{12}+{\epsilon }_{i-1,j-1}^{12}}{4h_1{h}_2}=0,$$

(51)

$$(\lambda +2\mu )\frac{{\epsilon }_{i,j+1}^{22}-2{\epsilon }_{i,j}^{22}+{\epsilon }_{i,j-1}^{22}}{h_2^2}+\lambda \frac{{\epsilon }_{i,j+1}^{11}-2{\epsilon }_{i,j}^{11}+{\epsilon }_{i,j-1}^{11}}{h_2^2}+2\mu \frac{{\epsilon }_{i+1,j+1}^{12}-{\epsilon }_{i-1,j+1}^{12}-{\epsilon }_{i+1,j-1}^{12}+{\epsilon }_{i-1,j-1}^{12}}{4h_1{h}_2}=0,$$

(52)

$$\begin{array}{c}\mu (\frac{{\epsilon }_{i+1,j}^{12}-2{\epsilon }_{i,j}^{12}+{\epsilon }_{i,j}^{12}}{h_1^2}+\frac{{\epsilon }_{i,j+1}^{12}-2{\epsilon }_{i,j}^{12}+{\epsilon }_{i,j-1}^{12}}{h_2^2})+(\lambda +\mu )(\frac{{\epsilon }_{i+1,j+1}^{11}-{\epsilon }_{i+1,j-1}^{11}-{\epsilon }_{i-1,j+1}^{11}+{\epsilon }_{i-1,j-1}^{11}}{4h_1{h}_2}\\ +\frac{{\epsilon }_{i+1,j+1}^{22}-{\epsilon }_{i+1,j-1}^{22}-{\epsilon }_{i-1,j+1}^{22}+{\epsilon }_{i-1,j-1}^{22}}{4h_1{h}_2})=0,\end{array}$$

(53)

To solve the difference Equations (51)–(53), taking into account the boundary conditions (34) and (35), it is convenient to use the elimination method. To solve this, we express Equation (51) in the following tridiagonal form [27]:

$$a_i{\epsilon }_{{}_{i+1,j}}^{11}+b_i{\epsilon }_{{}_{i,j}}^{11}+c_i{\epsilon }_{{}_{i-1,j}}^{11}=f_{ij}^x$$

(54)

$$\left\{\begin{array}{l}{\alpha }_{01}{\epsilon }_{{}_{0,j}}^{11}+{\beta }_{01}{\epsilon }_{{}_{1,j}}^{11}={\gamma }_{01}\\ {\alpha }_{02}{\epsilon }_{{}_{n-1,j}}^{11}+{\beta }_{02}{\epsilon }_{{}_{n,j}}^{11}={\gamma }_{02}\end{array}\right.$$

(55)

where

$$a_i=\frac{\lambda +2\mu }{h_1^2}, b_i=\frac{-2(\lambda +2\mu )}{h_1^2}, c_i=\frac{\lambda +2\mu }{h_1^2},$$

$$f_{ij}^x=-\lambda \frac{{\epsilon }_{i+1,j}^{22}-2{\epsilon }_{i,j}^{22}+{\epsilon }_{i-1,j}^{22}}{h_1^2}-2\mu \frac{{\epsilon }_{i+1,j+1}^{12}-{\epsilon }_{i-1,j+1}^{12}-{\epsilon }_{i+1,j-1}^{12}+{\epsilon }_{i-1,j-1}^{12}}{4h_1{h}_2}$$

Considering (55) and the boundary conditions in (34), we observe the following [28]:

$$\begin{array}{l}{\alpha }_{01}=1, {\beta }_{01}=0, {\gamma }_{01}=0,\\ {\alpha }_{02}=0, {\beta }_{02}=1, {\gamma }_{02}=0\end{array}$$

(56)

Equations (52) and (53), similarly to (54), can be reduced to a tridiagonal form with different coefficients as follows:

$$\begin{array}{l}{a}_i^{\prime }{\epsilon }_{i,j+1}^{22}+b_i^{\prime }{\epsilon }_{i,j}^{22}+c_i^{\prime }{\epsilon }_{i,j-1}^{22}=f_{ij}^y,\\ {\dot{a}}_i{\epsilon }_{{}_{i+1,j}}^{12}+{\dot{b}}_i{\epsilon }_{{}_{i,j}}^{12}+{\dot{c}}_i{\epsilon }_{{}_{i-1,j}}^{12}=f_{ij}^{xx},\\ {\tilde{a}}_i{\epsilon }_{{}_{i,j+1}}^{12}+{\tilde{b}}_i{\epsilon }_{{}_{i,j}}^{12}+{\tilde{c}}_i{\epsilon }_{{}_{i,j-1}}^{12}=f_{ij}^{yy}\end{array}$$

(57)

From Equations (54) and (57), the solution of the difference Equations (51) and (53) follows a successive fourfold application of the elimination method. The first two equations were solved using the elimination method over the indices $i$ and $j$, respectively, and the third equation was solved using $i, j$. According to [29], this solution method is called the variable-direction method.

6. Numerical Examples

This section describes the numerical solution of the plane boundary-value problems B, C, and D in strains and a comparison of the results with each other as well as with the well-known Timoshenko solution [18].

Let a rectangular plate with dimensions (2a, 2b) be under the action of a uniaxial parabolic load applied from opposite sides perpendicular to the OX axis [18]. The remaining sides were free from loads as follows Figure 2:

$$for\hspace{1em}x=\pm a:\hspace{1em}{\sigma }_{11}=S_0(1-\frac{y^2}{a^2}), {\sigma }_{12}=0,$$

(58)

$$for\hspace{1em}y=\pm b:\hspace{1em}{\sigma }_{22}=0, {\sigma }_{21}=0.$$

(59)

For the problem under consideration in the study by Timoshenko–Goodier [18], based on the condition for minimizing the strain energy using the Airy stress function, the following expressions were obtained for the components of the stress tensor:

$$\begin{array}{l}{\sigma }_{11}=S(1-\frac{y^2}{a^2})-0.1702S(1-\frac{3y^2}{a^2}){(1-\frac{x^2}{a^2})}^2,\\ {\sigma }_{22}=-0.1702S(1-\frac{3x^2}{a^2}){(1-\frac{y^2}{a^2})}^2,\\ {\sigma }_{12}=-0.6805S\frac{xy}{a^2}(1-\frac{x^2}{a^2})(1-\frac{y^2}{a^2}).\end{array}$$

(60)

The initial data have the following dimensionless values:

$$\lambda =0.8, \mu =0.5, l_1=2a, l_2=2b, a=b=1, N_1=N_2=10.$$

Table 1. Stress values ${\sigma }_{11}$ according to Problem A: (60) [18].

	x = −1	x = −0.8	x = −0.6	x = −0.4	x = −0.2	x = 0
y = −1	0.0000	0.0441	0.1394	0.2402	0.3137	0.3404
y = −0.8	0.3600	0.3803	0.4241	0.4705	0.5043	0.5166
y = −0.6	0.6400	0.6418	0.6456	0.6496	0.6525	0.6536
y = −0.4	0.8400	0.8285	0.8037	0.7776	0.7584	0.7515
y = −0.2	0.9600	0.9406	0.8987	0.8543	0.8220	0.8102
y = 0	1.0000	0.9779	0.9303	0.8799	0.8431	0.8298

lists the stress values ${\sigma }_{11}$ in one-quarter of a rectangular plate based on the results of Timoshenko–Goodier [18].

Table 2. Strain values ${\epsilon }_{11}$ according to Problem A: (60) [18].

	x = −1	x = −0.8	x = −0.6	x = −0.4	x = −0.2	x = 0
y = −1	0.0000	0.0300	0.0948	0.1633	0.2132	0.2314
y = −0.8	0.2297	0.2516	0.2877	0.3237	0.3494	0.3586
y = −0.6	0.3876	0.4144	0.4369	0.4644	0.4680	0.4644
y = −0.4	0.4893	0.5256	0.5431	0.5498	0.5515	0.5516
y = −0.2	0.5459	0.5903	0.6066	0.6084	0.6056	0.6041
y = 0	0.5640	0.6115	0.6277	0.6282	0.6240	0.6219

lists the strain distribution ${\epsilon }_{11}$ of the slabs. Deformations were calculated from the stress (60) based on Hooke’s law.

Table 3. Stress distribution ${\sigma }_{11}$ in the plate at y = 0, according to Problems A, B, and C.

y = 0	x = −1	x = −0.8	x = −0.6	x = −0.4	x = −0.2	x = 0
Problem A	1.0000	0.9779	0.9303	0.8799	0.8431	0.8298
Problem B (k = 68)	1.0000	0.9714	0.9434	0.8751	0.8522	0.8378
Problem C (k = 62)	1.0000	0.9691	0.9424	0.8769	0.8542	0.8404

lists the stress values ${\sigma }_{11}$ in the section y = 0 of a rectangular plate, obtained as a result of the numerical solution of boundary-value problems A, B, and C. Boundary-value problems B and C were solved by the iterative method and required 68 and 62 iterations, respectively.

The values of the stresses ${\sigma }_{11}$ and strains ${\epsilon }_{11}$ in section y = 0, provided in the first two rows of

Table 4. Stress distribution ${\sigma }_{11}$ in the plate at y = 0, according to Problems B, C, and D.

y = 0	x = −1	x = −0.8	x = −0.6	x = −0.4	x = −0.2	x = 0
Problem B (k = 80)	1.0000	0.9928	0.9869	0.9789	0.9766	0.9789
Problem C (k = 84)	1.0000	0.9908	0.9813	0.9722	0.979	0.9751
Problem D	1.0000	0.9818	0.9818	0.9818	0.9818	0.9818

and

Table 5. Strain distribution ${\epsilon }_{11}$ in the plate at y = 0, according to Problems B, C, and D.

y = 0	x = −1	x = −0.8	x = −0.6	x = −0.4	x = −0.2	x = 0
Problem B (80)	0.6797	0.6709	0.6621	0.6618	0.6605	0.6589
Problem C (84)	0.6797	0.6737	0.6689	0.6640	0.6622	0.6603
Problem D	0.6797	0.6797	0.6797	0.6797	0.6797	0.6797

, were obtained by solving Problems B and C using the iterative method, where k denotes the number of iterations. The stress values from ${\sigma }_{11}$ are provided in the third line of the table. Four is determined by solving Problem D using the variable direction method [30]. According to the variable direction method [21], the solution of finite difference Equations (51)–(53) was reduced to the sequential application of the elimination method to solve Equations (54)–(57). As can be observed from Table 4, the stress values ${\sigma }_{11}$ at y = 0 tend to be the maximum value of the specified load $S_0=1.$ Figure 3 and Figure 4 show the distribution of stresses in the plate based on the results of Timoshenko–Goodier [18] (Problem A) and the solution of the boundary-value problem D. Figure 4 shows that the stress distribution based on the results of task D is more accurate and closer to the maximum value of the given load.

7. Conclusions

To formulate a spatial boundary value problem during deformation, it is sufficient to consider the first or second group of differential deformation equations in combination with three equilibrium equations with the corresponding three surfaces and three additional boundary conditions;

• Equilibrium equations expressed with respect to deformations can be considered as additional boundary conditions at the boundary of a given region;
• The correct formulation of plane boundary value problems, consists of two equilibrium equations and compatibility conditions ${\epsilon }_{11}{}_{,22}+{\epsilon }_{22}{}_{,11}=2{\epsilon }_{12}{}_{,12}$ or a third equation of the two-dimensional strain differential equations with two-boundary and one additional boundary conditions;
• The plane boundary value problems formulated using two equilibrium equations and third strain differential equation with a corresponding boundary condition are more suitable for numerical solutions direct the strain tensor components;
• In the formulation of the plane boundary value problems in strains the equilibrium equations expressed with a strain may be used in a differentiated form, which allows to increase in the order of approximation of finite-difference equations;
• The finite difference method is convenient for satisfying additional boundary conditions;
• Grid equations for plane problems (B, C, and D) in strains were compiled using the finite difference method and solved using the iteration and variable direction methods;
• The validity of the formulated plane boundary-value problems in strains and the reliability of the results obtained were substantiated by comparing the numerical results of Problems B, C, and D and with the well-known Timoshenko–Goodier solution for the tension of a rectangular plate with a parabolic load applied on opposite sides.
• The considered methodology can be used in formulating and solving coupled thermoelasticity and thermoplasticity problems, as well as considering strain rates and specifying boundary conditions regarding strains.