Научная статья на тему 'FINITE ELEMENT MODELS IN STRESSES FOR PLANE ELASTICITY PROBLEMS'

FINITE ELEMENT MODELS IN STRESSES FOR PLANE ELASTICITY PROBLEMS Текст научной статьи по специальности «Медицинские технологии»

CC BY
105
20
i Надоели баннеры? Вы всегда можете отключить рекламу.
Журнал
Magazine of Civil Engineering
Scopus
ВАК
RSCI
ESCI
Ключевые слова
FINITE ELEMENTS / STRESS APPROXIMATION / FUNCTIONAL OF ADDITIONAL ENERGY / LAGRANGE MULTIPLIERS

Аннотация научной статьи по медицинским технологиям, автор научной работы — Tyukalov Yu.Ya.

The solution of the plane problems of elasticity theory on the basis of stress approximation is considered. To construct the solution, the additional energy functional is used. With the help of the principle of possible displacements, algebraic equations of equilibrium of the nodes of the grid of finite elements are constructed. Equilibrium equations are included in the functional of additional energy by means of Lagrange multipliers. The necessary relations for rectangular and triangular finite elements are obtained. Variants with constant and piecewise-constant approximations of stresses in the region of the finite element are considered. The ribbon width of system of the solving linear equations is estimated. Calculations have been made for the bended beam and for stretched plate with the hole, for the different grids of finite element. It is made comparison of the solutions obtained in stresses with the solutions obtained by finite element method in displacements and with exact solutions. It is shown, that for plane problems in the theory of elasticity, solutions based on stress approximations make it possible to obtain convergence of displacements to exact values from above. For coarse grids, solutions based on piecewise constant stresses much more accurate results, but require large computational costs, since the width of the ribbon of non-zero elements of the resolving system of linear algebraic equations is approximately twice as large as in the other considered variants. Finite elements models in stresses allow constructing solutions, which are alternative to solutions obtained by finite element method in displacements.

i Надоели баннеры? Вы всегда можете отключить рекламу.
iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.
i Надоели баннеры? Вы всегда можете отключить рекламу.

КОНЕЧНО ЭЛЕМЕНТНЫЕ МОДЕЛИ В НАПРЯЖЕНИЯХ ДЛЯ ЗАДАЧ ПЛОСКОЙ ТЕОРИИ УПРУГОСТИ

Рассмотрено решение плоских задач теории упругости на основе аппроксимации напряжений. Для построения решения используется функционал дополнительной энергии. При помощи принципа возможных перемещений составляются алгебраические уравнения равновесия узлов сетки конечных элементов. Уравнения равновесия включаются в функционал дополнительной энергии при помощи множителей Лагранжа. Получены необходимые соотношения для прямоугольных и треугольных конечных элементов. Рассмотрены варианты с постоянными и кусочно-постоянными аппроксимациями напряжений по области конечного элемента. Даны оценки ширины ленты системы линейных разрешающих уравнений. Выполнены расчеты на изгиб консольной балки и на растяжение пластины с отверстием при различных сетках конечных элементов. Выполнено сравнение полученных решений с решениями по методу конечных элементов в перемещениях и точными решениями. Показано, что для плоских задач теории упругости решения на основе аппроксимации напряжений позволяют получить сходимость перемещений к точным значениям сверху. При грубых сетках решения на основе кусочно-постоянных напряжений дает существенно более точные результаты, но требуют больших вычислительных затрат, так как ширина ленты ненулевых элементов разрешающей системы линейных алгебраических уравнений примерно в два раза больше, по сравнению с другими рассмотренными вариантами. Модели конечных элементов в напряжениях позволяют строить решения, альтернативные решениям, полученным методом конечных элементов в перемещениях.

Текст научной работы на тему «FINITE ELEMENT MODELS IN STRESSES FOR PLANE ELASTICITY PROBLEMS»

doi: 10.18720/MCE.77.3

Finite element models in stresses for plane elasticity problems

Конечно элементные модели в напряжениях для задач плоской теории упругости

Yu.Ya. Tyukalov,

Vyatka State University, Kirov, Russia

Key words: finite elements; stress approximation; functional of additional energy; Lagrange multipliers

Д-р техн. наук, профессор Ю.Я. Тюкалов,

Вятский государственный университет, г. Киров, Россия

Ключевые слова: конечные элементы; аппроксимация напряжений; функционал дополнительной энергии; множители Лагранжа

Abstract. The solution of the plane problems of elasticity theory on the basis of stress approximation is considered. To construct the solution, the additional energy functional is used. With the help of the principle of possible displacements, algebraic equations of equilibrium of the nodes of the grid of finite elements are constructed. Equilibrium equations are included in the functional of additional energy by means of Lagrange multipliers. The necessary relations for rectangular and triangular finite elements are obtained. Variants with constant and piecewise-constant approximations of stresses in the region of the finite element are considered. The ribbon width of system of the solving linear equations is estimated. Calculations have been made for the bended beam and for stretched plate with the hole, for the different grids of finite element. It is made comparison of the solutions obtained in stresses with the solutions obtained by finite element method in displacements and with exact solutions. It is shown, that for plane problems in the theory of elasticity, solutions based on stress approximations make it possible to obtain convergence of displacements to exact values from above. For coarse grids, solutions based on piecewise constant stresses much more accurate results, but require large computational costs, since the width of the ribbon of non-zero elements of the resolving system of linear algebraic equations is approximately twice as large as in the other considered variants. Finite elements models in stresses allow constructing solutions, which are alternative to solutions obtained by finite element method in displacements.

Аннотация. Рассмотрено решение плоских задач теории упругости на основе аппроксимации напряжений. Для построения решения используется функционал дополнительной энергии. При помощи принципа возможных перемещений составляются алгебраические уравнения равновесия узлов сетки конечных элементов. Уравнения равновесия включаются в функционал дополнительной энергии при помощи множителей Лагранжа. Получены необходимые соотношения для прямоугольных и треугольных конечных элементов. Рассмотрены варианты с постоянными и кусочно-постоянными аппроксимациями напряжений по области конечного элемента. Даны оценки ширины ленты системы линейных разрешающих уравнений. Выполнены расчеты на изгиб консольной балки и на растяжение пластины с отверстием при различных сетках конечных элементов. Выполнено сравнение полученных решений с решениями по методу конечных элементов в перемещениях и точными решениями. Показано, что для плоских задач теории упругости решения на основе аппроксимации напряжений позволяют получить сходимость перемещений к точным значениям сверху. При грубых сетках решения на основе кусочно-постоянных напряжений дает существенно более точные результаты, но требуют больших вычислительных затрат, так как ширина ленты ненулевых элементов разрешающей системы линейных алгебраических уравнений примерно в два раза больше, по сравнению с другими рассмотренными вариантами. Модели конечных элементов в напряжениях позволяют строить решения, альтернативные решениям, полученным методом конечных элементов в перемещениях.

1. Introduction

Many fundamental works have been devoted to the development of the theoretical foundations of the finite element method in displacements [1-4]. They present the basic variational principles based on which finite element solutions can be constructed. The principles of minimum potential energy and additional energy are considered. Hybrid and mixed approaches are also considered. In these studies, it is noted that solutions based on the principle of minimum potential energy give the lower limit of the solution, and based on the principle of additional energy can give an upper bound. Applying different approximations for displacements in the finite element region, we thereby reduce the number of degrees of freedom of the system under consideration, which leads to an increase in its rigidity. Therefore, the values of displacements determined by the finite element method in displacements. will always be less than the exact values. In addition, as a rule, the approximations used for displacements do not ensure the continuity of deformations, and therefore stresses, along the boundaries of finite elements. This leads to the appearance of gaps of stress fields along the boundaries of finite elements and, correspondingly, additional difficulties and inaccuracies in solving physically nonlinear problems. Special algorithms for calculating stresses are required [5]. In [3], methods based on the stress fields approximations are also considered.

Based on the finite element method, various algorithms for solving geometrically and physically nonlinear problems of rods, plates, shells, and bulk problems of the theory of elasticity are developed [611]. The finite element method in displacements is also successfully used for solving problems with allowance for geometric nonlinearity, problems with taking account of the shear deformations and for calculation the thin-walled constructions [12-14]. The finite element method also conveniently solves the problems of stability and the dynamics of structures [10, 15].

The mixed variants of the finite element method are considered in [17-22]. In mixed versions, approximations of both displacements and stresses (forces) are used. Therefore, the solution obtained by the mixed finite element method, when crushing the grid, can approach an exact solution both from below and from above and does not give either the lower or upper boundaries of the displacements. In [23], hybrid finite element models are considered, in which the equilibrium equations inside the elements are satisfied on average, and the resolving equations are reduced to the form, which allowing to apply the standard techniques used in the rigidity method. The paper [24] is devoted to the application of a mixed approach to the solution of plane problems in the theory of elasticity, in which piecewise constant approximations are used for displacements. Works [25, 26] are devoted to a semi-analytical method for solving the problems of building mechanics. The proposed method has a high solution accuracy, but has a limited application area.

In [27-32], solutions are constructed by the finite element method, which based on the approximation of stresses (forces). In [27], a combination of the principles of possible displacements and possible stressed states is used to obtain the solution. In [29-32], the solution is based on the use of the principle of minimum of the additional energy and the principle of possible displacements. This approach allows us to find solutions that are alternative solutions obtained by the finite element method in displacements, and can provide a lower bound of displacements.

The finite element method in displacements gives an approximate, one-sided solution of the problem under consideration. Therefore, despite the great successes in the use of the finite element method in displacements, the search for and development of additional, alternative solutions is topical.

The purpose of this paper is to develop the method for calculating planar rod systems based on various approximations of stresses over the region of the finite element.

2. Methods

The solutions of the plane problems of the theory of elasticity in stresses can be obtained using the functional of additional energy [3]:

nc = U* + V* =1j[a]T[E]-1{a}dn - J[T]T{A}dS ^ min,

(1)

U* - additional energy of the strains, V*- potential of boundary forces corresponding to the specified displacements [1]; {A} - vector given displacements of nodes; [T] - vector boundary forces; S -

boundary surface, on which the displacement nodes are given; 0 - subject area; (a) = \ay } -

stresses vector; [Я] - material stiffness matrix.

у

I ^xy)

The functions approximating the stress fields in (1) must satisfy the differential equations of equilibrium. Since, in the general case, it is practically impossible to find such functions, the following approach is proposed in [29-32]. The construction is dividing into rectangular or triangular finite elements. The stress fields in the region of the finite element can be approximated by linear, constant, or piecewise constant functions (Fig. 1). Linear approximating functions (Fig. 1a) ensure continuity of stresses fields over the entire subject area. Constant functions (Fig. 1b) are discontinuous along the boundaries of finite elements, but they satisfy differential equations of equilibrium in the absence of distributed loads. Piecewise-constant functions (Fig. 1c) also satisfy differential equations of equilibrium, are continuous along boundaries but have discontinuities inside finite elements.

Figure 1. Variants of the approximation of stresses over the region of the finite element: a) the stresses vary linearly; b) the stresses are constant; c) the stresses are piecewise constant.

For simplicity, we assume that there are no given node displacements. Then, using any variant of the approximating functions for the stresses (Fig. 1), the expression for the functional (1) can be written in the following matrix form:

nc = i(a)T[D](a) ^ min, (2)

(a) is the vector of unknown stresses for the whole system; [D] is the matrix of flexibility for the whole system. Then, using the principle of possible displacements, for all non-fixed directions of nodes along the coordinate axes, we compose algebraic equations of equilibrium of forces:

[Ci,x}T(oi) + Pi,x = 0, iESx, {Ci,y}T(oi) + Pi,y = 0, iEEy.

(o"j) - vector of unknown stresses of finite elements, connected to node i; Ex, Ey - sets of nodes that have non-fixed displacements along the X, Y axes, respectively; Pix, Piy - generalized forces, corresponding to the potential of external loads, for possible unit displacements of node i along the axes X, Y; [Cix], [Ciy] - vectors containing coefficients for unknown nodal stresses in the equations of equilibrium of the node i along the axes X, Y.

Thus, we have obtained the problem of minimizing the quadratic function of several variables (2) in the presence of constraints in the form of the system of linear algebraic equations (3). The unknown parameters are either nodal stresses (Figs. 1a, 1c) or stresses in the finite element (Fig. 1b). The solution of this problem was considered in [31, 32] using the penalty function method (4):

nc = 1(a)T[D](a) + lj=Xiyliesj « ([Cijf (°i) + Pi,jf ^ min. (4)

a — penalty parameter. Equating the derivatives, along the unknown stresses, to zero, we obtain the system of linear algebraic equations. The matrix of coefficients of this system of equations will have a ribbon structure of non-zero elements for any variant of stresses approximations.

In this paper, to solve we use the method of Lagrange multipliers:

nc =1{a}T[D][a} + I,j=XiyIlieSjUij ({CiJ}T{ai} + PiJ) ^ min, (5)

Uij - displacement of node i in direction j. In this solution, it is appeared additional unknowns in the form of node displacements. But we must emphasize that the approximation of the displacement fields in the region of the finite element is not used in (5).

Expression (5) can be represented in a more convenient form for constructing the solution:

1

Uc = -{a}T[D]{a} + [u}T({F] - [L]{a}) ^ min, (6)

(u) - global vector of unknown nodal displacements; {F} - vector, whose elements are equal to the work of external forces on the corresponding unit displacements of the nodes; [L] - "equilibrium" matrix, whose rows are formed from the corresponding vectors {C^}. If we equate the derivatives nc along the vector {a} to zero, we obtain the equations of compatibility of the deformations in the stresses:

[D]{a} - [L]T{u) = 0. (7)

The derivatives nc along the vector {u} represent the system of equations of the equilibrium of

nodes

{F} - [L]{a} = 0. (8)

Combining (7) and (8), we obtain the following system of linear algebraic equations:

[D] -[L]T] —L] [0]

Expressing the vector {a} from the first matrix equation and substituting it into the second, we

obtain

[K] = [L][D]-1[L]T, (10)

[K]{u} = {F}, CM)

{a} = [D]-1[L]T{u}. (12)

Thus, solving the system of algebraic equations (11), we obtain the values of the node displacements {u}, and then the stress vector {a} from (12).

Next, we obtain the necessary expressions for the elements of the matrices [D], [L] and the vector {F}, entering in (6), when rectangular and triangular finite elements are used to discredit the subject area (Fig. 2).

If linear functions are used to approximate stresses (Fig. 1a), then the matrix will be filled and the solution of the system of linear algebraic equations (11) requires large computational costs. Therefore, below, two variants of approximating the stresses in the region of the finite element will be considered - constant and piecewise-constant functions. In these cases, the matrix will have a ribbon structure of non-zero elements.

ZM-M »

Figure 2. Triangular and rectangular finite elements

2.1. Variant 1. The stresses are constant in the region of the finite element

In the case of constant stresses in the region of the finite element, the vector of unknowns for the k-th finite element in the local coordinate system (~ok) will have three elements:

[ôt} =

'—к ■ Ox

—к

Oy

^Xyу

(13)

The symbols with the top dash denoted the parameters related to the local coordinate system, associated with the finite element. The superscript k is the serial number of the finite element. The vectors of stresses in the global { <rfe}and local {<7fc} coordinate systems are connected by the matrix of direction cosines:

к

■ 2

s i n2a

2

cos2a

2

с os2 a

2

s in2a

2s ina • с os a 2s ina • с os a

22

i—s ina • cos a s ina • cos a —s in2a — cos2a

, {0к} = №]{ок].

(14)

a - angle between the vertical axes of the global and local coordinate systems (Fig. 2).

The additional strain energy for a finite element of thickness tk from an isotropic material is determined by the area integral:

Пск = ик=^1к{о}т[Е]-1{о}с1А,

(15)

where:

O

{o} =

O

y

xy)

(16)

1 0 -n 1 0 0 0 2(1 +^)\

Substituting ( ak) in {14} instead of (a), we obtain an expression for the additional deformation energy in the global coordinate system

Uk = ±Ak tk(ak)T[E]-1(ak), (17)

where: Ak- area of the k-th finite element. The local matrix of a finite element flexibility has following form:

[Dк] =Ак гк[Е]-1.

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

(18)

Summing the additional deformation energies of all n finite elements, we obtain the following expression for the flexibility matrix [D] for the whole system:

[D] =

[D1]

0

[Dn]

(19)

The matrix [D] will have block-diagonal structure and is therefore easily reversible. This is important, since the inverse matrix [D]-1 is subsequently will used to obtain the matrix [tf] (10).

[ D]

-i _

[D1]-1

0

0

[Dn]

n -1

(20)

For rectangular finite element, we introduce the local coordinate system ^c], connected with its center (Fig. 2), and basic functions, which are expressed in normalized local coordinates in the following form:

Ni(x,y) =

(l + teXi + mv)

4 '

2x 2y

Ç = —, V=T~. i = 1,2,3,4.

ak bk

(21)

The index i denotes the local order number of node of the finite element; x,y - the coordinates along the axes X1 and Y1, respectively; - local normalized coordinates of node i, taking values 1 or -1. The nodes of finite element are numbered counterclockwise, beginning with the lower-left node.

As possible displacements of each node, consider the displacements SUi and Svi along the axes of the local coordinate system X1 and Y^ Possible displacements in the region of rectangular finite element are expressed by means of the linear approximation functions in the following form:

Su(x,y) = Ni(x,y)öUi, Sv(x,y) = Nt(x,y)Svi.

(22)

Figure 3. Possible displacements of node i in the global coordinate system

Since the possible displacements SUi and Svi can be any, we take them equal to unity and in subsequent transformations we omit them. Then, the deformations arising in the element k at the possible displacement SUi = 1, directed along the X1 axis will be as follows:

Ö£y =

d(SUU) ti(1 + vvi)

dx

2aL

ÖYxy =

d(SUU) _Vi(1 + Ki) dy 2bk .

(23)

The possible energy of deformations of the rectangular finite element with number k, arising at the possible displacement SUi = 1

—k

- f0k C tk(*xÔ£x + TxyÖYxy)dxdy

We write expression (24) in vector form:

tkbk c —k . tkak -k

■ >- —k , +

vi

xy-

(24)

' tkbk

—k

SUkx = [Ckxfin [CÜ H ~0

Si

(25)

k ak

Vi}

0

2

2

2

Similarly, for the possible displacement ÔVi = 1

—k

= {Ckyftfl {Cky} =

tkak

t^bk ^ 2

■Vi

'■St

(26)

The energy of deformations of all rectangular finite elements adjoined to the node under consideration is determined in the form of sums

SUiiX = ïk[ckx}T[ëk], SUiy = ïk[cky}T[ëk}.

(27)

The potential of external concentrated and uniformly distributed loads for possible displacements of node i along the global coordinate axes is determined by (28).

SVi,j = pi,j qfakbk = Rij, j = x,y.

(28)

Pij - forces, which concentrated in the node; qk - load, which evenly distributed over the element. The generalized forces Ry taken with the minus sign, are placed in the vector {F} (see (6)). Using the expressions (25) and (26), we introduce the following matrix notations:

r —k

SUi*

—k

SUl,y —k

SU2,X —k SU k,y

—k SU^x —k SU2,y —k SU4,x —k VS U4,yJ

Expressions for the energy of deformation for possible displacements of nodes of a finite element will be written in the matrix form:

И=

r—ki tk

, M=7

— bk 0 — ak

0 — ak — bk

bk 0 — ak

0 — ak bk

bk 0 ak

0 ak bk

— bk 0 ak

0 ak — bk

(29)

(5Uk) = [bk] {âk} = [bk] [Sk][ok].

(30)

To obtain the equilibrium equations for a node, it is necessary to consider the possible displacements ôu^ôVi along the global coordinate axes X and Y (Fig. 3). Displacements along global and local axes are connected by the matrix of direction cosines

[ ]=[ Li

s in a — с os a с os a s in a

],

(31)

Using (31), we obtain the expressions for the strain energy for possible displacements of the nodes along the axes of the global coordinate system:

[5Uk] = [Lk][ak], [Lk] = [Skv][Lk][Sk], [Skv] =

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

[ ] 0 0 0

0 [ ] 0 0

0 0 [ ] 0

0 0 0 [ ]

(32)

0

2

The matrix [Lk], conditionally, can be called the local matrix of equilibrium of the finite element. From matrices [Lk] for finite elements, in accordance with the numbering of the nodes and elements, the global equilibrium matrix [ L] is formed for the whole system.

To derive the relations for the triangular finite element, we use triangular coordinates [3], which allow us to obtain the matrix [Lk] directly in the global coordinate system. The triangular coordinate for the finite element Tk is determined by the formula

^^^(af + bfx + cty), i = 1,2,3.

where:

a

k _

xi+1yi+2 xi+2yi+1,

(33)

bk = yi+i

yi+2,

k

xi+2 xi + 1.

Ak-the area of the triangular finite element k; Xi, yi - coordinates of node i.

We can express possible displacements in the region of the triangular finite element by means of triangular coordinates in the following form:

Su(x,y) = Tjkôui, Sv(x,y) = TkSvi.

(34)

Sui,Svi - the possible displacements of the node i along the global coordinate axes. We take Sui = 1,Svi = 1. Then, the deformations arising in the element at the possible displacement of the node along the global X axis, will be as follows:

ôsx =

d(Su) dx

= xv = wu. =

2Ak' Yxy dy 2Ak'

Then

tkf,k

SU*X = ff tk(a£Sex + TkySyxy)dA =~ILaxk + Similarly, for the possible displacement along the Y axis, we get:

tkc?

2 'Txy'

Ô£y =

d(Sv) dy

k

2At

S Yxy =

d(Sv) bk

d x

2 A

k

SUjky = JJ tk(akSSy + TkySYxy)dA =

tkck tkbk

--ak +--— Tk

uy + 2 Lxy.

2

(35)

(36)

(37)

(38)

Using the expressions (36) and (38), we obtain the following matrix of equilibrium [Lk] of the triangular finite element in the global coordinate system:

\bk

m=v

0

0

bk o

k C1

bk

k

bk 0

k 2

k c3

2 b k

k c3

bk]

(39)

The potential of external concentrated and uniformly distributed loads for possible displacements of node i is determined by (40).

SVij = Pij + ilk q]jAk = Ri,p j = x,y.

(40)

The global equilibrium matrix [L] for the whole system will have ribbon structure of non-zero elements. Numbering of unknown stresses are assigned in accordance with the numbering of the finite elements. Therefore, the width of the ribbon of non-zero elements of the matrix [L] will be determined by the maximum difference of order numbers of the finite elements adjoined to node. The structure of the matrix [L] must be considered when storing its elements, as well as when calculating the matrix [tf] (10).

k

1

0

0

We note, that as a result, the width of the ribbon of non-zero elements of the matrix [tf] equal with the width of the ribbon of the system of equations of the finite element method in displacements.

2.2. Variant 2. The stresses are piecewise-constant in the region of the finite

element

In this case, nodal stresses are taken as unknowns. The stresses fields in the region of the finite

element are discontinuous and piecewise-constant (Fig. 1c). We Introduce the notation: {<7i} = { &i,y

\T i,xy, ( °i,x ]

- stresses in node i in the local coordinate system; { ai} = { °i,y } - stresses in node i in the global

\T i,xy)

coordinate system. The expression for the additional strains energy in the global coordinate system can be written in the form of a simple sum:

tiAi{ai]T[E]-l{ai],

¿^-4=1

л . — yUr ZAS + упт 1 лs Ai=bs=i4A + bs=itA .

пт

^=13

(41)

(42)

m - the total number of nodes; ^ - the thickness of the plate at node i, nR - the number of rectangular elements adjoined to node i; nT - number of triangular elements adjoined to node i; As - the area of s -the finite element.

We introduce the notation for the flexibility matrix of "neighborhoods" of the node i:

[Di] = tiAi[E]-1.

(43)

The global flexibility matrix for the entire system [D] consists of matrices of flexibility [Dt] for all nodes of the system and has the following block-diagonal form:

[D] =

[Di] 0

0

[Dm]\

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

(44)

Consider a possible displacement of the node Sut = 1 along the local axis Xt of the rectangular finite element. The possible energy of deformations of the rectangular finite element with number k

SUi>x = ¡Цк С tk(ax5sx + TxySYxy)dxdy = (l + oXiJ + tf^m (l +

+

(45)

2 ) ХУ'У

Similarly, for the possible displacement Svi = l

-—k _ rbk гак k

SUi,y = С С tk(bS£y + TxySYxy)dxdy = ZU^rn (l + + Zj^fi (l +

+-г)т

(46)

xy,]-

In expressions (45) and (46), j is the local order number of node in the finite element. Using (45) and (46), we write the expression for the matrix of equilibrium for the rectangular finite element [L ].

t-k

-3bk 0 -3ak -3bk 0 -ak -bk 0 -ak -bk 0 -3ak

0 -3ak -3bk 0 -ak -3bk 0 -ak -bk 0 -3ak -bk

3bk 0 -ak 3bk 0 -3ak bk 0 -3ak bk 0 -ak

0 — ak 3bk 0 -3ak 3bk 0 -3ak bk 0 -ak bk

bk 0 ak bk 0 3ak 3bk 0 3ak 3bk 0 ak

0 ak bk 0 3ak bk 0 3ak 3bk 0 ak 3bk

-bk 0 3ak -bk 0 ak -3bk 0 ak -3bk 0 3ak

0 3ak -bk 0 ak -bk 0 ak -3bk 0 3ak -3bk

(47)

We write expressions for the strain energy, at possible unit displacements of the nodes of the finite element k in the local coordinate system in a matrix form analogous to (30):

[8Uk}=[Lk]{ak} = [Lk] [Cka][ok].

(48)

[Ck] - the matrix of transformation of the unknown stress vector for the finite element from the local to the global coordinate system. This matrix consists of the matrices [S^] (14) and has the following form:

k

№ 0 0 0

0 №] 0 0

0 0 [sk] 0

0 0 0 [sk]

(49)

Going to possible displacements in the global coordinate system, we obtain the matrix of equilibrium for the finite element

[Lk] = [Skv] [Lk] [C*].

(50)

From the matrices [Lk], for all finite elements, a global matrix of equilibrium of the whole system [L] is formed.

For a triangular finite element, expressions like expressions (45) and (46), but in the global coordinate system, will have the following form:

suk suk

tkb*

i,x = 6 yj=l°x,j +

tkck

l-yj=1T.

xyj'

tkck

i,y = 6 ¿Jj = luy,j +

tkbk

^yj T ■

Ai j = l Lxy,j i

(51)

(52)

In a triangular element, the stresses are constant in each third of the area adjoined to the node, and equal to the stresses in the given node. Since the possible deformations, in accordance with expressions (35) and (37), are constant values throughout the entire element, the strain energy is equal to the sum of multiplies the stresses, deformations and the area of one third of the element. Using (51) and (52), we obtain the following matrix of equilibrium [Lk] of the finite element in the global coordinate system:

k l

k k

[Lk]=T

bk 0 k i bk 0 k i bk 0 i

0 k i bk 0 k i bk 0 k i bi

bk 0 rk 2 bk 0 rk 2 bk 0 2

0 k 2 bk 0 k 2 bk 0 k 2 b

bk 0 k bk 0 k bk 0

0 k j j bk 0 k j j bk 0 k j j b

(53)

2

2

k j

k

Expressions of the potential of external concentrated and uniformly distributed loads for possible displacements of the node (28) and (40) do not depend on the type of stress approximation.

The global equilibrium matrix [ L] will also have a ribbon structure of non-zero elements. When determining the width of the ribbon of non-zero elements, the minimum and maximum node numbers of all finite elements adjacent to the node under consideration are determined for each node. The maximum difference of these numbers, multiplied by three, will determine the width of the tape of nonzero elements of the matrix [ L]. Note, that the width of the ribbon of non-zero elements of the matrix [K] is approximately twice the width of the ribbon of the system of equations of the finite element method in displacements.

6

6

3. Results and Discussion

By program developed in Mathcad 14.0, calculations of the cantilever beam for the action of a uniformly distributed vertical load were performed (Fig. 4b). The thickness of the beam - 1m., the modulus of elasticity - E = 10000 kN/m2, the coefficient of transverse deformations - 0.25. The results of calculations of the cantilever beam are presented in Figures 4-7 and in Tables 1-2. In Figures,

solutions based on piecewise constant approximations of stresses are denoted by digit 1; solutions based on constant stresses - denoted by digit 2; solutions based on FEM in displacements - denoted by digit 3.

The results of calculations for five variants of finite element grids show that solutions based on stress approximation allow obtaining convergence of displacements to the exact values from above and are more flexibility as compared to solutions obtained by FEM in the LIRA-SAPR displacements (Fig. 4a).

The stresses in the clamp of the beam with the use of piecewise constant stresses approximations are always greater than the stresses obtained using both constant stresses and with the use of FEM in displacements (Table 2). This is since for piecewise constant stresses approximations we use, as unknowns, the stresses at the nodes of the finite element grid, which allows us to obtain stress values directly in external fibers and better model the edge effect. The difference in the stress values (for the smallest grid), in comparison with the FEM in the displacements, is 15.6 %, 35.6 % and 26.4 %, respectively, for the stresses ax

oy and тху.

When using constant stresses in the region of the finite element, the stress ax is also always greater than the stresses obtained by the FEM in displacements. For the coarsest grid, the values differ by 9.3 %, and for the smallest by 1.5 %. The stresses ay and Txy are smaller in the case of using constant stresses. The stresses ay differ by 39 % for the coarsest grid and by 29 % for the shallowest grid. The stresses Txy differ by 5 % and 15 %, respectively.

Note that the stress ax in this case is much greater than the stresses ay and Txy. In addition, due to the edge effect, when we crash the grid, all the stresses increase, and it is not possible to compare the values obtained with precise values. When the solution of the problem is constructed in a physically nonlinear formulation, then solutions based on the approximation of stresses, giving large stresses for the same grids, should provide a greater reserve of strength in comparison with the FEM decisions in displacements.

Figure 4. Cantilever beam: a) vertical movements; b) finite element grid 12x4

Figure 5. The stress ax in the upper fiber in the clamp

Figure 6. Stresses ay in the upper fiber in the clamp

Figure 7. Stresses rxy in the upper fiber in the clamp

Table 1. Vertical displacements of the console

N Grid Piecewise constant stresses Constant stresses FEM in displacements LIRA-SAPR

Ribbon width w, mm Ribbon width w, mm Ribbon width w, mm

1 12x4 26 24.5166 14 24.5578 14 22.4760

2 24x8 42 23.6121 22 23.5639 22 23.0436

3 48x16 74 23.3552 38 23.3325 38 23.1990

4 96x32 138 23.2845 70 23.2762 70 23.2410

5 192x64 266 23.2653 134 23.2624 134 23.2531

Table 2. Stresses in the upper fiber of the section in the clamp

Piecewise constant stresses Constant stresses FEM in displacements LIRA-SAPR

N Grid ay, T xy, Oy, T xy, ay, T xy,

kN/m2 kN/m2 kN/m2 kN/m2 kN/m2 kN/m2 kN/m2 kN/m2 kN/m2

1 12x4 17.556 3.118 1.745 13.957 0.954 1.113 12.652 1.571 1.172

2 24x8 19.548 3.521 2.402 16.527 1.433 1.447 16.005 2.120 1.631

3 48x16 22.063 4.031 3.215 18.911 1.783 1.931 18.584 2.565 2.239

4 96x32 25.160 4.645 4.055 21.447 2.102 2.486 21.123 2.994 2.912

5 192x64 28.909 5.375 4.927 24.399 2.439 3.081 24.020 3.460 3.623

Also, calculations were made for the stretched square plate with the hole. In Figure 8 the finite-element grids are shown for the one-fourth of the plate. The size of the side of plate is 10m., the diameter of hole is 1m., the thickness is 1m., the modulus of elasticity is E = 10000 kN/m2, the coefficient of transverse deformations is 0.25. The stretching load is q = 10 kN/m2. For this problem, the exact values of the stresses at points 1 and 2 are known (Fig. 8). At the point 1 ax = 30 kN/m2, at the point 2 oy = -10 kH/m2 [3]. Table 2 shows the stress values obtained at points 1 and 2. In the variant in Fig. 8a only triangular finite elements are used and the number of grid nodes is 322, in the variant in Fig. 8b both rectangular and triangular elements are used, and the number of nodes is 2579.

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

Figure 8. One quarter of the square plate with the hole. Variants of finite element grid and stress distribution ax, obtained by FEM in displacements

The results obtained show that for the coarser grid (scheme a), a solution based on piecewise constant stresses give substantially more accurate results. The solution differs from the exact one by 7.7 % at point 1 and by 20 % at point 2. For variants with constant stresses and for FEM in displacements, the differs from the exact solution for point 1 is 13.8 % and 47.6 % for point 2. For the triangular finite element, the FEM in displacements also simulates constant stresses, so the solutions for the scheme in Figure 8a, which obtained by the FEM in displacements, and obtained on the basis of constant stresses, coincided.

For the fine grid (Fig. 8b), the solutions obtained by all methods are very close. In this case, solutions based on the approximation of stresses, are also closer to an exact solution.

4. Conclusions

1. For plane problems of the theory of elasticity, solutions based on the approximation of stresses make it possible to obtain convergence of displacements to the exact values from above.

2. For coarse grids, solutions based on piecewise constant stresses much more accurate results, but require large computational costs, since the width of the ribbon of non-zero elements of the resolving system of linear algebraic equations is approximately twice as large as in the other considered variants.

3. Solutions based on the approximation of stresses make it possible, with the same grids, to obtain more accurate stress values in comparison with the FEM solutions in displacements. By difference of solutions, by the method based on the approximation of stresses and FEM in displacements for identical grids, one can estimate the accuracy of the approximate solutions obtained.

4. Finite elements models in stresses allow constructing solutions, which are alternative to solutions obtained by finite element method in displacements.

Table 2. Stresses in the square plate with the hole

Scheme Point 1 - ах, kN/m2 Point 2 - ay, kN/m2

Piecewise constant stresses Constant stresses FEM in displacements, LIRA-SAPR Piecewise constant stresses Constant stresses FEM in displacements, LIRA-SAPR

a) 27.697 25.853 25.853 -8.003 -5.244 -5.244

b) 29.871 29.865 29.845 -9.290 -9.577 -9.559

References

1. Zenkevich O. Metod konechnykh elementov v tekhnike [The finite element method in techique]. Moscow: Mir, 1975. 541 p.(rus)

2. Zenkevich O., Morgan K. Konechnyye elementy i approksimatsiya [Finite Elements and Approximation]. Moscow: Mir, 1986. 318 p.(rus)

3. Gallager R. Metod konechnykh elementov. Osnovy [Finite element method. Basics]. Moscow: Mir, 1984. 428 p.(rus)

4. Rozin L.A. Metod konechnykh elementov v primenenii k uprugim sistemam [Finite element method in application to elastic systems]. Moscow: Stroyizdat, 1977. 129 p.(rus)

5. Stein E., Ahmad R. An equilibrium method for stress calculation using finite element displacement models. Comput. Meth. Appl. Mech. And Eng. 1977. Vol. 10. No. 2. Pp. 175-198.

6. Nguen-Xuan H. A polygonal finite element method for plate analysis. Computers & Structures. 2017. Vol. 188. Pp. 45-62.

7. Di Paola M., Fileccia Scimemi G. Finite element method on fractional visco-elastic frames. Computers & Structures. 2016. Vol. 164. Pp. 15-22.

8. Lyashko A.D., Tayupov Sh.I., Timerbayev M.R. Skhemy metoda konechnykh elementov vysokogo poryadka tochnosti dlya sistemy ellipticheskikh uravneniy s vyrozhdayushchimisya koeffitsiyentami na interval [Schemes of the method of finite elements of high order of accuracy for a system of elliptic equations with degenerate

Литература

1. Зенкевич О. Метод конечных элементов в технике. М.: Мир, 1975. 541 с.

2. Зенкевич О., Морган К. Конечные элементы и аппроксимация. М.: Мир, 1986. 318 с.

3. Галлагер Р. Метод конечных элементов. Основы. М.: Мир, 1984. 428 с.

4. Розин Л. А. Метод конечных элементов в применении к упругим системам. М.: Стройиздат, 1977. 129 с.

5. Stein E., Ahmad R. An equilibrium method for stress calculation using finite element displacement models // Comput. Meth. Appl. Mech. And Eng. 1977. Vol. 10. № 2. Pp. 175-198.

6. Nguen-Xuan H. A polygonal finite element method for plate analysis // Computers & Structures. 2017. Vol. 188. Pp. 45-62.

7. Di Paola M., Fileccia Scimemi G. Finite element method on fractional visco-elastic frames // Computers & Structures. 2016. Vol. 164. Pp. 15-22.

8. Ляшко А.Д., Таюпов Ш.И., Тимербаев М.Р. Схемы метода конечных элементов высокого порядка точности для системы эллиптических уравнений с вырождающимися коэффициентами на интервале // Известия высших учебных заведений. Математика. 2009. № 7. С. 22-34.

9. Лалин В.В., Рыбаков В.А., Морозов С.А.. Исследование конечных элементов для расчета тонкостенных стержневых систем // Инженерно-строительный

coefficients on the interval]. Izvestiya vysshikh uchebnykh zavedeniy. Matematika. 2009. No. 7. Pp. 22-34.(rus)

9. Lalin V.V., Rybakov V.A., Morozov S.A.. Issledovaniye konechnykh elementov dlya rascheta tonkostennykh sterzhnevykh sistem [Finite element analysis for the calculation of thin-walled rod systems]. Magazine of Civil Engineering. 2012. Vol. 27. No. 1. Pp. 53-73.(rus)

10. Lalin V.V., Rozin L.A., Kushova D.A. Variatsionnaya postanovka ploskoy zadachi geometricheski nelineynogo deformirovaniya i ustoychivosti uprugikh sterzhney [Variational formulation of the plane problem of geometrically nonlinear deformation and stability of elastic rods]. Magazine of Civil Engineering. 2013. No. 1(36). Pp. 87-96.(rus)

11. Lalin V.V., Zdanchuk Ye.V., Kushova D.A., Rozin L.A. Variatsionnyye postanovki nelineynykh zadach s nezavisimymi vrashchatelnymi stepenyami svobody[Variational formulations of nonlinear problems with independent rotational degrees of freedom]. Magazine of Civil Engineering. 2015. No. 4(56). Pp. 54-65.(rus)

12. Lalin V.V., Belyayev M.O. Izgib geometricheski nelineynogo konsolnogo sterzhnya. Resheniye po teoriyam Kirkhgofa i Kossera - Timoshenko[Bending of a geometrically nonlinear cantilever rod. Decision by the theories of Kirchhoff and Kosser - Tymoshenko]. Magazine of Civil Engineering. 2015. No. 1(53). Pp. 39-55.(rus)

13. Lalin V., Rybakov V., Alexander S. The finite elements for design of frame of thin-walled beams. Applied Mechanics and Materials. 2014. Vol. 578-579. Pp. 858-863.

14. Dyakov S.F., Lalin V.V. Postroyeniye i analiz konechnogo elementa tonkostennogo sterzhnya s uchetom deformatsiy sdviga dlya resheniya zadach dinamiki[Construction and analysis of a finite element of a thin-walled rod with allowance for shear deformations for solving problems of dynamics]. Internet-zhurnal Naukovedeniye. 2013. No. 5(18). Pp. 93. (rus)

15. Perelmuter A.V., Slivker V.I. Ustoychivost ravnovesiya konstruktsiy i rodstvennyye problem [Stability of structural equilibrium and related problems]. Vol. 1. Moscow: SKAD SOFT, 2010. 704 p.

16. Gordeyeva A.O., Vatin N.I. Raschetnaya konechno-elementnaya model kholodnognutogo perforirovannogo tonkostennogo sterzhnya v programmno-vychislitelnom komplekse SCAD Office [Estimated finite element model of a cold-bended perforated thin-walled rod in the software complex SCAD Office]. Magazine of Civil Engineering. 2011. No. 3(21). Pp. 36-46.(rus)

17. Mu L., Wang J., Ye X. A hybridized formulation for the weak Galerkin mixed finite element method. Journal of Computational and Applied Mathematics. 2016. Vol. 307. Pp. 335-345.

18. Ignatyev A.V., Yefremova N.S. Raschet kosougolnoy plastiny po metodu konechnykh elementov v forme klassicheskogo smeshannogo metoda [Calculation of an skew plate by the finite element method in the form of a classical mixed method]. Stroitelnaya mekhanika i konstruktsii. 2016. Vol. 1. No. 12. Pp. 39-44.(rus)

19. Ignatyev A.V., Ignatyev V.A., Bochkov M.I. Primeneniye metoda konechnykh elementov v forme klassicheskogo smeshannogo metoda k raschetu sistem s odnostoronnimi svyazyami [The application of the finite element method in the form of a classical mixed method to the calculation of systems with unilateral constraints]. Stroitelnaya mekhanika i raschet sooruzheniy. 2017. No. 2(271). Pp. 52-61.(rus)

20. Verbitskiy V.V., Verbitskiy A.V. Smeshannyy metod konechnykh elementov v zadache o pologoy obolochke s funktsiyey napryazheniy [A mixed finite element method in the problem of a shallow shell with a stress function]. Izvestiya vysshikh uchebnykh zavedeniy. Matematika. 2009. No. 5. Pp. 13-23. (rus)

21. Kulikov G.M., Plotnikova S.V. A hybrid-mixed four-node

журнал. 2012. № 1(27). С. 53-73.

10. Лалин В.В., Розин Л.А., Кушова Д.А. Вариационная постановка плоской задачи геометрически нелинейного деформирования и устойчивости упругих стержней // Инженерно-строительный журнал. 2013. № 1(36). С. 87-96.

11. Лалин В.В., Зданчук Е.В., Кушова Д.А., Розин Л.А. Вариационные постановки нелинейных задач с независимыми вращательными степенями свободы // Инженерно-строительный журнал. 2015. № 4(56). С. 54-65.

12. Лалин В.В., Беляев М.О. Изгиб геометрически нелинейного консольного стержня. Решение по теориям Кирхгофа и Коссера - Тимошенко // Инженерно-строительный журнал. 2015. № 1(53). С. 39-55.

13. Lalin V., Rybakov V., Alexander S. The finite elements for design of frame of thin-walled beams // Applied Mechanics and Materials. 2014. Vol. 578-579. Pp. 858-863.

14. Дьяков С.Ф., Лалин В.В. Построение и анализ конечного элемента тонкостенного стержня с учетом деформаций сдвига для решения задач динамики // Интернет-журнал Науковедение. 2013. № 5(18). С. 93.

15. Перельмутер А.В., Сливкер В.И. Устойчивость равновесия конструкций и родственные проблемы. Т. 1. М.: СКАД СОФТ, 2010. 704 с.

16. Гордеева А.О., Ватин Н.И. Расчетная конечно-элементная модель холодногнутого перфорированного тонкостенного стержня в программно-вычислительном комплексе SCAD Office // Инженерно-строительный журнал. 2011. № 3(21). С. 36-46.

17. Mu L., Wang J., Ye X. A hybridized formulation for the weak Galerkin mixed finite element method // Journal of Computational and Applied Mathematics. 2016. Vol. 307. Pp. 335-345.

18. Игнатьев А.В., Ефремова Н.С. Расчёт косоугольной пластины по методу конечных элементов в форме классического смешанного метода // Строительная механика и конструкции. 2016. Т. 1. № 12. С. 39-44.

19. Игнатьев А.В., Игнатьев В.А., Бочков М.И. Применение метода конечных элементов в форме классического смешанного метода к расчету систем с односторонними связями // Строительная механика и расчет сооружений. 2017. № 2(271). С. 52-61.

20. Вербицкий В.В., Вербицкий А.В. Смешанный метод конечных элементов в задаче о пологой оболочке с функцией напряжений // Известия высших учебных заведений. Математика. 2009. № 5. С. 13-23.

21. Kulikov G.M., Plotnikova S.V. A hybrid-mixed four-node quadrilateral plate element based on sampling surfaces method for 3D stress analysis // International journal for numerical methods in engineering. 2016. Vol. 108. No. 1. Pp. 26-54.

22. Stupishin L.U., Nikitin K.E. Mixed finite element for geometrically nonlinear orthotropic shallow shells of revolution // Advanced Materials Research. 2014. Vol. 919. Pp. 1299-1302.

23. Wolf J.P. Alternate hybrid stress finite element models // International Journal for Numerical Methods in Engineering. 1975. Vol. 9. № 3. Pp. 601-615.

24. Qiu W., Demkowicz L. Mixed hp-finite element method for linear elasticity with weakly imposed symmetry: stability analysis // SIAM J. Numer. Anal. 2011. Vol. 49. № 2. Pp. 619-641.

25. Elakkad A., Bennani M.A., Mekkaoui J., Elkhalfi A. A mixed finite element method for elasticity problem // International Journal of Advanced Computer Science and Applications. 2013. Vol. 4. No. 2. Pp. 161-166.

26. Akimov P.A., Negrozov O.A. Semianalytical structural analysis based on combined application of finite element method and discrete-continual finite element method part 1:

quadrilateral plate element based on sampling surfaces method for 3D stress analysis. International journal for numerical methods in engineering. 2016. Vol. 108. No. 1. Pp. 26-54.

22. Stupishin L.U., Nikitin K.E. Mixed finite element for geometrically nonlinear orthotropic shallow shells of revolution. Advanced Materials Research. 2014. Vol. 919. Pp. 1299-1302.

23. Wolf J.P. Alternate hybrid stress finite element models. International Journal for Numerical Methods in Engineering. 1975. Vol. 9. No. 3. Pp. 601-615.

24. Qiu W., Demkowicz L. Mixed hp-finite element method for linear elasticity with weakly imposed symmetry: stability analysis. SIAM J. Numer. Anal. 2011. Vol. 49. No. 2. Pp. 619-641.

25. Elakkad A., Bennani M.A., Mekkaoui J., Elkhalfi A. A mixed finite element method for elasticity problem. International Journal of Advanced Computer Science and Applications. 2013. Vol. 4. No. 2. Pp. 161-166.

26. Akimov P.A., Negrozov O.A.. Semianalytical structural analysis based on combined application of finite element method and discrete-continual finite element method part 1: two-dimensional theory of elasticity. Procedia Engineering. 2016. Vol. 153. Pp. 8-15.

27. Akimov P.A., Negrozov O.A.. Semianalytical structural analysis based on combined application of finite element method and discrete-continual finite element method part 2: three-dimensional theory of elasticity. Procedia Engineering. 2016. Vol. 153. Pp. 16-23.

28. Gienke E. A simple "mixed" method for plate and shell problems. Nuclear Engineering and Design. 1974. Vol. 29. No. 1. Pp. 141-155.

29. Kuo Y.-L. Stress-based finite element analysis of sliding beams. Appl. Math. Inf. Sci. 2015. Vol. 9. No. 2. Pp. 609-616.

30. Tyukalov Yu.Ya. Stress finite element models for determining the frequencies of free oscillations. Magazine of Civil Engineering. 2016. No. 7(67). Pp. 39-54.

31. Tyukalov Yu.Ya. The functional of additional energy for stability analysis of spatial rod systems. Magazine of Civil Engineering. 2017. No. 2(70). Pp. 18-32.

32. Tyukalov Yu.Ya. Raschet obolochek proizvolnoy formy metodom konechnykh elementov v napryazheniyakh [Calculation of shells of arbitrary shape by the finite element method in stresses]. Stroitelnaya mekhanika i raschet sooruzheniy. 2006. No. 1. Pp. 65-74. (rus)

33. Tyukalov Yu.Ya. Resheniye ploskoy zadachi teorii uprugosti metodom konechnykh elementov v napryazheniyakh [Solution of the plane elasticity problems by the finite element method in stresses]. Stroitelnaya mekhanika i raschet sooruzheniy. 2006. No. 2. Pp. 34-38.(rus)

two-dimensional theory of elasticity // Procedia Engineering. 2016. Vol. 153. Pp. 8-15.

27. Akimov P.A., Negrozov O.A.. Semianalytical structural analysis based on combined application of finite element method and discrete-continual finite element method part 2: three-dimensional theory of elasticity // Procedia Engineering. 2016. Vol. 153. Pp. 16-23.

28. Gienke E. A simple "mixed" method for plate and shell problems. Nuclear Engineering and Design. 1974. Vol. 29. № 1. Pp. 141-155.

29. Kuo Y.-L. Stress-based finite element analysis of sliding beams // Appl. Math. Inf. Sci. 2015. Vol. 9. № 2. Pp. 609-616.

30. Тюкалов Ю.Я. Определение частот свободных колебаний методом конечных элементов в напряжениях // Инженерно-строительный журнал. 2016. № 7(67). С. 39-54

31. Тюкалов Ю.Я. Функционал дополнительной энергии для анализа устойчивости пространственных стержневых систем // Инженерно-строительный журнал. 2017. № 2(70). С. 18-32.

32. Тюкалов Ю.Я. Расчет оболочек произвольной формы методом конечных элементов в напряжениях // Строительная механика и расчет сооружений. 2006. № 1. С. 65-74.

33. Тюкалов Ю.Я. Решение плоской задачи теории упругости методом конечных элементов в напряжениях // Строительная механика и расчет сооружений. 2006. № 2. С. 34-38.

Yury Tyukalov,

+7(912)821-89-77; [email protected]

Юрий Яковлевич Тюкалов, +7(912)821-89-77; эл. почта: [email protected]

© Tyukalov Yu.Ya., 2018

i Надоели баннеры? Вы всегда можете отключить рекламу.