ARBITRARY QUADRANGULAR FINITE ELEMENT FOR PLATES WITH SHEAR DEFORMATIONS Текст научной статьи по специальности «Физика»

Magazine of Civil Engineering. 2021. 107(7). Article No. 10707

Magazine of Civil Engineering

journal homepage: http://engstroy.spbstu.ru/

ISSN 2712-8172

DOI: 10.34910/MCE.107.7

Arbitrary quadrangular finite element for plates with shear deformations

Yu.Ya. Tyukalov

Vyatka State University, Kirov, Russian Federation E-mail: [email protected]

Keywords: finite element method, flexible plates, shear deformations, approximation of moments, arbitrary quadrangular finite element

Abstract. An arbitrary quadrangular bending finite element based on a piecewise constant approximation of moments is proposed. The solution is based on the principles of minimum additional energy and possible displacements. The finite element allows you to consider the shear deformation, regardless of the ratio of the plate thickness to its sizes. The effect of locking in the calculation of thin plates is absent. Comparison of the results of the oblique plates calculation, annular and round plates with analytical solutions and calculation results for other programs was done. The comparison shows good accuracy in determining displacements and moments. Crushing the finite elements mesh makes the displacement values tend to exact values from above. To assess the influence of the finite element shape, the square plate calculations were made. To model the square plate, quadrangular elements were used, obtained from rectangular ones by changing the slope of one side. At the same time, the calculating accuracy of the displacements and moments decreased slightly. The proposed finite element is easy to implement. The problem solution did not require a numerical integration or the mapping of the quadrangular region to the rectangular one. The necessary expressions were obtained analytically.

Bendable plates are widely used as load-bearing structures of buildings and structures for various purposes. A lot of scientific literature has been devoted to the development of various types of finite elements for bendable plates [1-2]. To date, software systems use a wide range of finite elements for bendable plates, developed on the basis of displacement approximations and using the Lagrange's functional. Despite this, the development of new types of finite elements continues. In particular, papers [3-4] are devoted to new finite elements for the calculation of thin plates, developed on the basis of the general theory of elasticity.

Bendable plate calculation which consider shear deformations requires the use of finite elements based on various shear theories. In [5], the equations of the theory of elasticity are directly applied to allow shear deformations, and in [6-7] the Mindline-Reissner theory of plates is used. In contrast to the classical Kirchhoff's theory, in the Reissner's theory, rotation angles and vertical displacements are considered as independent variables, which allows one to lower the maximum order of derivatives in the strain energy functional and makes it possible to use first-order functions for displacements approximations. Direct use of the Reissner theory for constructing of finite elements leads to the "jamming" effect, which consists in the impossibility of using these finite elements to calculate thin plates. The "jamming" effect allows to use these finite elements only for thick plates. In [8], the new finite element is proposed, constructed on the basis of the Reissner plate theory, but suitable for calculating thin plates. To calculate a certain class of plates, there are analytical methods which use trigonometric series for the representation of the differential equation solutions of plate bend [9-10]. This approach allows us to obtain solutions with the necessary

Tyukalov, Yu.Ya. Arbitrary quadrangular finite element for plates with shear deformations. Magazine of Civil Engineering. 2021. 107(7). Article No. 10707. DOI: 10.34910/MCE.107.7

© Tyukalov, Yu.Ya., 2021. Published by Peter the Great St. Petersburg Polytechnic University. This work is licensed under a CC BY-NC 4.0

1. Introduction

accuracy for plates of a certain shape. In particular, such method is used to calculate periodic oscillations of square plates [11].

The emergence of new composite materials for plates led to the appearance of new finite elements types based on various shear theories [12-13]. For such plates the Galerkin's method is used to construct triangular and quadrangular finite elements according to the Reissner-Mindline theory of plates in [13]. In [14], the first-order shift theory was used, and displacements along the cross section vary linearly. Tasks of calculating various structures, including plates, which consider contact interaction with the base, are importance. In [15], on the basis of the finite-element contact model and the step-by-step analysis method, a numerical algorithm was developed, that allows simultaneous integration of the motion equations and the implementation of contact conditions with Coulomb friction.

For the calculation of the plates, in addition to the Lagrange functional, mixed and hybrid functionals are used [16-17]. The use of such formulations makes it possible to simplify the account of shear deformations due to the use of transverse forces and moments as unknowns, in addition to displacements. But such decisions require the observance of additional conditions to ensure convergence. In [18], a quadrangular finite element with seven degrees of freedom in each node was presented: three displacements of the median surface along the coordinate axes, two shear angles, and two normal rotation angles. This element allows to more accurately consider shear deformations when changing material properties in different directions. In [19], the finite element formulation free from the "locking" effect was presented for calculating plates with allowance for shear based on the modified Mindline theory. A separate class of problems is the plates stability analysis. To solve such problems, the finite element method is widely used, but analytical methods are also using. In [20], the analytical solution was obtained for stability analysis of thin rectangular plates. Such analytical solutions can also be used to test the finite element solutions.

In [21], the finite element model was proposed for calculating thin flexible plates, based on the moment's approximations. The solution was obtained by using the functional of the additional energy and the possible displacements principle. The resolving equations for rectangular and triangular finite elements was obtained. Such finite element models make it possible to get convergence of displacements from above. Similar solutions based on the stress' approximations were obtained for plane problems of the elasticity theory [22-23]. In [24], the method was proposed for considering shear deformations, based on the transverse forces' approximation and the additional energy functional. Shear displacements are determined independently of bend displacements. This approach is free from the locking effect and can be used to calculate any thickness plates. In [25], the finite element model also based on the moments approximations was proposed for calculating plates according of the Reissner theory.

Finite element solutions based on the moments (stresses) approximation allow one to obtain solutions that are alternative to solutions based on the Lagrange functional. This determines their relevance.

This paper aim is to build an arbitrary quadrangular finite element designed to calculate various thicknesses plates, based on the moments and transverse forces approximations. The solution will be based on papers [21, 24], in which only rectangular and triangular finite elements were considered. To evaluate the accuracy of the proposed finite element, test calculations and comparison of the results with analytical solutions and results obtained by other methods will be performed.

2. Methods

We will build the solution of the problem of bending thin plates to involve the moments approximations on the basis of the functional [1]:

1 , 12

/

2 AE • i3

nc = - j ——r M2 +MZ - 2vM M + 2(1 + v)Mz IdA, (1)

x y x y ^ ' xy

E is the elastic modulus of the plate material; t is the plate thickness; v is Poisson's ratio; M^ is bending moment directed along the X axis; M is bending moment directed along the Y axis; M is torque. We

y xy

write the functional (1) in the standard matrix form:

nc = - f MTE-1M dA. (2)

2 A s s

The following notations are introduced in expression (2):

M

M

M

y

M

xy

„-1 12 E 1 =--

E • t

1 -v 0 -v 1 0 0 0 2 (1 + v)

(3)

We introduce the notations for unknown nodal moments M ,M ,M . and for the nodal

X,l y, i xy,l

moments vector of an arbitrary quadrangular finite element

M

= \M лм лм ЛМ .M .M .M M M M .M .M . .. x, 1 y,1 xy, 1 x,2 y, 2 xy,2 x,3 y,3 xy,3 x,4 y,4 xy,4 )

(4)

In the finite element region, we will approximate the moments fields by piecewise constant functions (Fig. 1) [21, 24].

3 0.5/,

C)

Figure 1. An arbitrary quadrangular finite element: a) dashed lines go through the midpoints of

sides and divide the finite elements into regions with constant moments; b) the possible displacement of the node; c) the nodes local numbers and the coordinate system for the finite

element.

To simplify the expressions recording, we introduce auxiliary stair-step functions:

( ) i1, ( x, y)e A

W> (x,% (x, y )* A

(5)

And introduce the corresponding diagonal matrices

^ =

wг

wi

(6)

Then the moments approximations matrix in the finite element region will have simple form:

Z

к

v1 V2 V3 V4

We get

Using (8), we will express the additional strain energy of the finite element in the following form:

M = Z,M, .

s к к

пк=2 jM T ( zT E-1zk )M к ^

(7)

(8)

k

We introduce the notation for the local flexibility matrix of the finite element:

Dk = *Z kE 1z kdA-

Ak

Note that the matrix Dk is calculated analytically and has simple block-diagonal form:

(10)

D

k

A]E

-1

A2E

-1

A3E

-1

A4E

-1

(11)

Ai ^ A^ are the areas of the finite element parts in which the moments are constant (Fig. 1).

From the flexibility matrices of finite elements, the global flexibility matrix D is formed, and from s Mk the global nod

expression of functional (2):

vectors Mk the global nodes moments vector M. Using the introduced notation, we obtain the following

n c =1MT DM.

(12)

It is important that the matrix D is block-diagonal and consists of square matrices measuring 3 by 3. Therefore, the matrix D is easily analytically reversible. This circumstance greatly simplifies the construction of the problem solution.

In accordance with the principle of additional energy minimum, the moments fields must satisfy the equilibrium equations and static boundary conditions. To ensure the moments fields equilibrium, we compose the grid nodes equilibrium equations using the possible displacements principle [21, 24].

Fig. 1b shows the possible displacement of the node and the finite elements adjoined to it. To approximate vertical displacements in the finite element region, which are caused by the possible displacement of node i, we will use the following function:

Sw} (x, y) = «1,. + «2,+ «3, .y + 04 jXy,

i = 1 * 4,

(13)

x,y are point coordinates in the local coordinate system of finite element xy' (Fig. 1c). To determine the

parameters a . . included in the expression for the approximating functions (13), it is necessary to form the

J,i

following matrix:

B

1 X-

X1 y1 X1 y1

X2 y2 X2 y2

X3 y3 X3 y3

X4 y4 X4 y4

(14)

x.,y. are coordinates of node i of the finite element. If the possible displacement of node is taken equal to

i i

unity, then the parameters a . . are elements of the inverse matrix B

-1.

B

-1

a1,1 a1,2 a1,3 a1,4

a2,1 a2,2 a2,3 a2,4

a3,1 03,2 a3,3 03,4

a4,1 a4,2 a4,3 a4,4

(15)

Each column of the matrix B contains the parameters of approximating functions for possible displacement of one finite element node.

The principle of possible displacements can be written as follows:

SU. +SV. =0.

I I

(16)

SU is the internal forces work of the finite elements adjoined to the node i; SV is the external forces

i i

potential due to possible displacements. For one finite element k, the internal forces work, in the general case, is written as follows:

SUk = f i ,

f d2(Sw.) d2(Sw.) d2(Sw.) ^

M -+ M -+M i

X ,..2 y a..2

dx

From expression (12) we obtain:

8H

y dyz iswi

xy dxdy

dA.

(17)

J

d2(Sw.)

Sp .=—-—— = a„ .+ a. .y, S(p .=—-—— = o • + a. x, Sk .= 2-— = 2a. .,

^xi dx 2, i 4, f yi dy 3, i 4, i xyi dxdy 4, i

(18)

S(Pxj , S(Pvi are the possible rotation angles along the Xand Y axes, respectively; Sk . is possible

xyi

xi ' - yi torsional curvature.

Figure 2. Graphs of the rotation angles changing along the normal to the finite element sides.

Obviously, the bending curvatures

d2(dw.) d2(sw.)

2

2

are equal to zero and the finite element

dx dyz

rotates without bending. At the same time, the moments which normal to the sides of the finite element perform the work with the corresponding rotation angles. But the torques perform the work inside the finite element area. Therefore, we write the work of the finite element internal forces in the following form:

SUk = f M Sp dS + f M Sk dA

1 1 n n J n> nil

S

xy xyi

(19)

M is the moment normal to the borders (sides) of the finite element; Sp is the angle of the section

n n

rotation along the normal to the borders (sides) of the finite element; Sk are the finite element borders. In accordance with (18), torsion curvature is constant in the finite element region, therefore

f M Sk .dA = 2a . \M A k xy xyi 4,i y xy, 1 i

A"

+M A + M A + M A I. xy,2 2 xy,3 3 xy,4 4

Cosines and sines of the angles indicated in Fig. 2, are calculated by the node's coordinates:

y 2 - y1 y3- y 2 y4 - y3 y1 -y 4

cos a = —- , cosa = —- , cosa =- , cosa =

1 i 2 i 3 i 4

1

l.

X^ X _ X- X ~ X. X. X -

sina =—2, sin a =—3, sin a =—4, sin a =—1 1 l, 2 j 3 / ' 4 /

(21)

1 2 3 4

The section rotation angle along the normal to the side k of the finite element, by the possible displacement of node i, is calculated by the following formula:

dtp = 8p cosa. + d( sina. =( a ,+a .y) cosa. +(^ x) sin a • (22)

n X k y k \ 2, / 4, / J k \ 3, / 4, / J k

The points coordinates on the finite element sides indicated in Fig. 2 are determined by the following formulas:

3 X1 + X2

X / = 4 , X //

X1 + 3 X2

1

1

y1/ =

3 y1 + y2

y1// =

y1 + 3 y2

3 X2 + X3 X2 + 3 X3 3 y2 + y3

X ,=—--, X .. =—--, y .=—---, y ..

2/ 4 2// 4 2/ 4 2//

X = 3X3 + X4 X = X3 + 3X4 3/ 4 ' 3// 4

, y / =

37

3 y3 + y4

y

y2 + 3 y3

y3 + 3 y4

(23)

// 4

3 X. + X X. + 3 X 3 y. + y

= 4 1 = 4 1 = ^4 J1

X1 / = 4 , X4// = 4 , y4/ = 4 , y,//

y4 + 3 y1

4/ 4 4/ / 4 4/ 4 4/ / 4

Substituting the points coordinates in (22), we obtain the expressions for the rotation angles of the section, along the normal to the sides, for the points indicated in Fig. 2:

f \ f \

8(p

ni, j1

f

a_ .+ a. .y .

2, i 4, / ■ / J y

cos a + J

a. ,+a. .X .

3, i 4,i ■ / J

sina ., J

Sp

n, j "

A

a^ ,+a. .y , , 2, i 4, f •//

cosa + j

A

(24)

a, ,+a. .X , ,

3, i 4, i . // J

svna ., i, j = 1,2,3,4.

' j ; v

In (24), the index / denotes the node number, whose possible displacement is considered; index j denotes the side number of the finite element.

Bending moments directed along the normal to the finite element sides are calculated by the well-known formula:

2 2

M = —M sin a.— M cos a.— 2M sin a cos a.

(25)

n X l y l xy l l

In (25), moments enter with a minus sign, since the moments external to the boundary opposite to

the internal moments. Using (25), we calculate the normal moments values, that are constant for each half

of the finite element side:

Mn = -M , sin2 a-M , cos2 a - 2M , sin a cos a, 1,1 x,1 1 y,1 1 xy,1 1 1

M" =-M ,sin2a-M ,cos2a-2M ,sina cosa, 2,1 x,2 1 y,2 1 xy,2 1 1

Mn , = -M , sin2 a-M , cos2 a - 2M , sin aa cos a,

2,2 x,2 2 y,2 2 xy,2 2 2

M" = -M , sin2 a-M , cos2 a -2M , sin a cos a, (26)

3.2 x,3 2 y,3 2 xy,3 2 2

n 2 2

M , = -M ,, sin a -M ,, cos a -2M ,, sina cosa,

3.3 x,3 3 y,3 3 xy,3 3 3

M" = -M . sin2 a-M . cos2 a - 2M . sin a cos a,

4.3 x,4 3 y,4 3 xy,4 3 3

n 2 2

M. ,=-M .sin a.-M .cos a.-2M .sma. cosa,

4.4 x,4 4 y,4 4 xy,4 4 4

Mn . =-M ,sin2a-M ,cos2a.-2M sina cosa. 1,4 x,1 4 y,1 4 xy,1 4 4

The rotation angles along each side of the finite element vary linearly (Fig. 2), and the moments are constant on the half of each side, therefore

$ MnSpnds = Um" Sp i + M^Sp i + MrniSp i +Mr!Sp i + „k n n 2 1,1 ni, 1' 1 ni, 17/ 1 2 ,2 ni,2' 2 2,3 ni,2 ' 2

S (27)

M3 3S( 'l3 + M3 4S( ''l3 + M44S( 'l4 + m" 4S( ''l4).

3 , 3 ni, 3 3 3 , 4 ni, 3'' 3 4, 4 ni, 4 4 1 4 ni, 4//4 Summing up the integrals (20) and (27), we obtain the internal forces work SUk of the finite element on the possible displacement of the node i. The work can be written in vector form:

SUk = CTM , (28)

i i k

Ci is the vector whose elements are the total coefficients from expressions (19) and (27), which are

multipliers of the corresponding unknown moments. The vector Ci elements will have the following expressions:

2 2 Sp .L sin a1+Sp ..l sin a4

n/,1 1 n/,4 ' 4 4

2 2 Sp ,L cos a+Sp , ,l. cos a.

n/,1 1 n/,4' 4

2Sp ,L sina cosa + 2Sp ,,l. sina cosa- 2A,a. . n/,11 1 1 1 n/,411 4 4 4 1 4,/

2 2 Sp ,l sin a +Sp .A sin a

n/,2 2 2 n/,X ' 1 1

2 2 Sp , L cos a+Sp , A cos a

n/,21 2 2 n/,Y ' 1 1

2Sp ,L sina cosa + 2Sp . ,L sina cosa - 2A„a. . n/,2 2 2 n/,1//1 1 1 2 4,/

2 2 Sp , L sin a+Sp , A sin a

n/,3 3 3 n/,2' 2 2

2 2 Sp .Lcos a +Sp .A cos a

n/,3' 3 3 n/,27 7 2 2 2Sp ,L sina cosa + 2Sp .A. sina cosa- 2 Aa. .

n/,3 3 3 n/,2 2 2 ^4,/

2 2 Sp ,l. sin a ,+Sp , A sin a

n/,4' 4 4 n/,3'1 3 3

2 2 (29)

Sp . l cos a4 +Sp , , l3 cos a3 n/,4/ 4 4 n/,3'' 3 3

2Sp /1 sin a4 cos a4 + 2Sp . . l3 sin a3 cos a3 - 2A4a4 . n/,4' 4 4 4 n/,3''3 3 3 4 4,1

After calculating the internal forces work on the possible displacements of the all finite element nodes, we obtain

C =--i 2

5U] SU, SU.

SU

L ^M ^, L ^

4

(30)

The matrix Lk consists of four rows and twelve columns. From the finite elements matrices Lk, in

accordance with the nodes numbering, the global matrix L is formed for the entire system. The external forces potential by the possible displacement of node /

m

SV. = P .+ 2 I q Sw.dA = P ..

1 z,z k=1A z 1 z'z k

(31)

The integral in (31) is calculated over the area m of finite elements adjacent to node /'. P . is force

z, /

which is concentrated in the node; q is load which is distributed over the element area. For the case of a

z

uniformly distributed load, the integral is calculated analytically. As calculations showed, without loss of accuracy, the uniformly distributed load can be replaced by concentrated nodes. Then we get

m

SV. = P .+ Z q A j = P ., i z, i , z i, k z, i k = 1

(32)

A , is part of the finite element area adjacent to the node i (see Fig. 1). From the values P ., in

i, k z, i

accordance with the nodes numbering, the global vector of external forces P is formed. In accordance with (15), we obtain the nodes equilibrium equations system

LM + P = 0.

Using the Lagrange multipliers, equations (33) are added to the functional (11).

nc = l MT DM + wT (LM + P),

(33)

(34)

W Is the nodal displacements vector. To consider the static boundary conditions, it is necessary to compose additional equations and, using the Lagrange multipliers, add it to the functional (34). Then the vector will include these additional parameters. The necessary equations are given in [21, 24].

To obtain solving equations, we equate the derivatives (34) with respect to the vectors M and w to

zero:

(35)

DM + Lw = 0, LM + P = 0.

Expressing the vector M from the first equation and substituting it into the second equation, we

obtain:

K = LD 1LT, Kw = P,

M = D-1LTw.

(36)

When calculating thick plates, it is necessary consider additional shear deformations. In the framework of the proposed methodology, we can consider the shear state regardless for the bending state and use the following functional [24]:

nc =1 J 2 S

2k (1 + v)

E ■ t

qX+Q2y IdS,

(37)

k = 6/5 is coefficient which consider the parabolic law of the shear stresses change along the plate thickness. Q , Q are transverse forces. We introduce the notations:

Q =<

s

Q

y

E-1 = 12 (1 + v)

sh 5E ■ t

1 0 0 1

(38)

Then

nc =1J QTE—1Q dS.

yJ s sh s

(39)

We consider the case when the transverse forces are assumed constant over the finite element region. Then for one finite element we get:

n k 2 Qk Dsh, k Qk '

Ak is the finite element area; Qk =

Q

x, k

Q

y, k

is the finite element transverse forces vector. Then

Dsh, k Ak Esh.

For whole system, we can write

n=2 qT DshQ,

(41)

(42)

Q is vector of shear forces for whole system; D^ is global flexible matrix of the shear for whole system.

Using the possible displacements principle, we compose the nodes equilibrium equations under shear. To approximate possible displacements under shear, we also use expression (12). As a result of the node possible displacement, the shear deformations take place in the sections:

d(Sw , . ) d(Sw , .

c \ sh, i _ \ sh, i

oy = —-- = a_ . + a. .y, oy = —-

XZ dX 2, i 4, 'yz dy

= a. . + a. .X.

3,i 4,i

(43)

Then the work of the internal transverse forces for the finite element k, at the possible displacement of node /, can be expressed as follows:

(44)

^ = MSxzQx, k + SyzQy, k f =

Ak

= a2,/AkQx, k + a3,/AkQy, k + Qx, ka4, / J ydS + Qy, ka4, / J xdS

Ak Ak

To calculate the integrals (44), we divide the quadrangular finite element into two triangles 1 -2-3 and 1-3-4 (see Fig. 1). Next, we will perform integration over each part using triangular coordinates. Then we get the following simple expression:

a„

a

SUk = a2lAkQxk + a3lAkQyk +

/ ' ^ X, A tx y, /V

1 (A123 (y1 + y2 + y3) + A134 (y1 + y3 + y4))QX,k + r(A123 (X1 + X2 + X3) + A134 (X1 + X3 + X4))Qy,k'

(45)

^123 , ^134 are the areas of the triangles 1-2-3 and 1-3-4 (see Fig. 1). We write the work suk in matrix form:

U=csTh A'

(46)

C , is a vector whose elements are coefficients from expression (45), which are multipliers of the

sh ,l

corresponding unknown transverse forces.

c

sh, i

A ^

a_ .A, +—' 2, i k 3

, 4, i a_ .A, +—L 3, i k 3

f (A123 (y1 + y2 + y3 ) + A134 (y1 + y3 + y4 )) (A123 (X1 + X2 + X3 ) + A134 (X1 + X3 + X4 ))

(47)

After calculating the internal forces work on possible displacements of all nodes of finite element, we

obtain

SU1

su:

su:

SU

> L sh, k Qk, L sh, k

C C C C

sh,1 T

sh,2 T

sh,3 T

sh,4

(48)

The matrix L , . consists of four rows and two columns. From the finite element matrices L , , ,

sh, k sh, k

in accordance with the nodes numbering, the global matrix is formed for the whole system L^ .

The external forces potential is calculated by (32). The nodes equilibrium equations system under

shear:

LshQ + P =

(49)

Including the equilibrium equations in functional (42) with the help of Lagrange multipliers, we obtain

2

n C = 2 QTDshQ + wTh (L shQ + P ),

(50)

wsh is the nodal displacements vector, caused by the sections shear deformations. To obtain solving equations, we equate the derivatives (50) with respect to the vectors Q and w ^ to zero:

DshQ + L shw sh 0, L shQ + P = 0.

(51)

Expressing the vector Q from the first equation and substituting it into the second equation, we

obtain:

K = L D 1LT

sh sh sh sh

Ksh Wsh P Q DshLshwsh.

(52)

To consider the static boundary conditions, it is necessary to compose additional equations and, using the Lagrange multipliers, add into functional (50). Then the vector will include these additional unknowns. The necessary equations are given in [24]. Obviously, the total nodes displacements considering shear deformations are determined as the sum:

w = w + w , .

sum sh

(53)

3. Results and Discussions

As mesh refinement test, we calculated of the skew plate which is supported by hinges along the contour (Fig. 3). For such plate, the finite-difference solution is given in [26]. We used parallelogram finite elements with the largest angle of 120 degrees. The plate has the following parameters:

7 2 2

E = 10 kN/ m , ^, = 0.21, t = 0.1m, q = 10kN/m . The geometric dimensions are shown in Fig. 3.

The calculations were performed for three finite element grids: a 6 by 10 grid, a 10 by 20 grid, and a 20 by 40 grid. In addition, for comparison, the calculations of this plate were performed according to the LIRASCAD program. All results are summarized in Table 1, where SFEM is the solution according to the proposed method.

Figure 3. The skew plate with hinge supports.

Table 1. Displacements and bending moments in the plate center (Fig. 3). Mesh refinement test.

Solution (grid) w, mm Mx, kN ■ m / m My, kN ■ m / m

SFEM (6x10) 10.779 3.654 9.903

SFEM (10x20) 10.039 3.586 9.451

SFEM (20x40) 9.587 3.529 9.132

LIRA-SAPR (6x10) 9.093 3.290 8.303

LIRA-SAPR (10x20) 9.115 3.404 8.611

LIRA-SAPR (20x40) 9.149 3.449 8.731

Timoshenko [26] 9.719 - 8.712

The results in Table 1 show that the plate center displacement obtained using the proposed finite element (SFEM), when crushing the mesh tends to the exact value from above. Similar results were obtained for rectangular and triangular finite elements, also constructed on the basis of the stress approximation [22]. For the smallest grid, the obtained displacement value differs from the analytical one by 1.4 percent. The bending moments also tend to exact values from above. The value of the maximum moment, obtained on the smallest grid, is 4.8 percent more than the analytical value.

As a further example, the half-ring was calculated on the uniformly distributed load action (Fig. 6). The half-ring has a hinged support along the lines DE and DB, and clamped support along the line EA. Line AB is the symmetry axis, therefore, at the nodes lying on this line, zero torques were taken. The following data were accepted in the calculations:

E = 10000 kN/m2, /j = 0.3, t = 0.1m, q = 10kN/m2,R = 6m, r = 3m.

¡¡¡¡¡¡mi

um»»»».!:;:::

a) b)

Figure 4. Finite-element meshes of the half-ring.

The obtained results were reduced to dimensionless form:

_ 1000D -

w = w-

qR

M = x

M _ M

X Tir _ y

qR

0,M =

X y

qR

(54)

Table 2. Comparison of finite element results for the half-ring with the analytical solution.

Solution Wc My,A Mxc My ,c M x B

SFEM (Fig. 4a) 0.374 0.0134 0.0387 -0.00460 -0.0173 -0.00226

SFEM (Fig. 4b) 0.362 0.0133 0.0396 -0.00445 -0.0169 -0.00394

LIRA-SAPR (Fig. 4a) 0.350 0.0098 0.0287 -0.00422 -0.0163 -0.00352

LIRA-SAPR (Fig. 4b) 0.348 0.0108 0.0336 -0.00434 -0.0166 -0.00291

Analytical [27] 0,358 0,0118 0,0393 -0,00439 -0,0168 -0.00234

Analysis of the results shown in Table 2, demonstrates a good accuracy of the proposed finite element. The value of the radial direction moment m ., which has the prevailing values, is very close

y,a

to the analytical value. The maximum displacement obtained by the proposed method differs from the analytical value by less than one percent. The moment value in the circumferential direction is an order of magnitude smaller than the moment value in the radial direction. The proposed finite element gives a significantly greater value of this moment. It should be noted that this moment value obtained by the LIRA-SAPR program is also much greater than the analytical value. We also note that the values of another

moment M „ in the circumferential direction, obtained both by the proposed method and by the LIRA-

x, C

SAPR program, are close to the analytical value.

To assess the influence of the finite element shape (mesh distortion test), the square plate that have clamped support along the contour was calculated on the action of the uniformly distributed load. Figure 5 shows the finite element grid for the plates quarter. The calculations were performed for grids with angles from zero to 60 degrees.

The sides of the AC and AB are pinched; the sides of CD and BD are symmetry axes.

The following data were accepted for the plate:

2 2 E = 10000 kN / m , ¡j. = 0.3, t = 1m, q = 10kN / m . Table 3 shows the results of calculations of the

plate without considering shear deformations. Obviously, then larger the angle /3, the more distorted the finite element shape, and the more significantly the calculation results should differ from the values obtained for the grid of square finite elements (/ = 0).

Table 3. Calculation results for the square plate (Fig. 5). Mesh distortion test.

ß, degrees wD, mm MxD, kNm/ m M y D, kNm / m MyB, kNm / m

0 18.537 8.361 8.361 -18.251

10 18.527 8.303 8.282 -18.143

20 18.516 8.230 8.184 -18.024

30 18.503 8.140 8.072 -17.901

40 18.488 8.037 7.953 -17.774

50 18.469 7.930 7.834 -17.649

60 18.446 7.823 7.723 -17.535

First of all, we note that the displacement of the plate center depends on the angle / insignificant. At 3 = 60 degrees, the displacement value is less than the displacement for the square grid only on 0.5

percent. The value of the moment M in the clamped side, taken modulus, is less than the absolute

y B

value of this moment at B = 0 on 3.9 percent. The value of the moment M ^ , at B = 60 , decreases

x, D

by 6.4 percent. In general, we can conclude that even with a strong distortion of the finite elements shape, when the maximum internal angle becomes equal to 150 degrees, the proposed arbitrary quadrangular finite element allows us to obtain sufficiently accurate values of displacements and moments for the square-shaped plate.

As the next test example (shear locking test), the calculations of the hinged and clamped support round plates (Fig. 6) on action of a uniformly distributed load with considering shear deformations were performed. For such plates, there are Timoshenko-Mindline analytical solutions center displacements.

a) b) c)

Figure 6. Finite element mesh for a quarter of the round plate.

For a hinged round plate, this Timoshenko-Mindline decision is written as follows:

->2

qR4 (5 + ß qR w = ---— +-

62D

1 + ß

4kGt

For clamped round plate:

qR4 qR2 w = --+ -

62D 4kGt

R is the plate radius; t is the plate thickness; D = ■

Et~

12 (l V )

(55)

(56) , 5

is the plate cylindrical stiffness; k = — ;

E

G =—--- is shear modulus. The following parameters have been used for the plate:

2 (1 + m)

7 2 2

E = 10 kN / m , ¡j.= 0.3, q = 10kN / m , R = 3m . The calculation results are summarized on Table 4

and Table 5.

Table 4. The hinged round plate center displacement (Fig. 6). Shear locking test.

w, mm

t, m Fig. 6a Fig. 6b Fig. 6c Formula (55)

0.1 57.0493 56.4158 56.3331 56.4158

0.5 0.46989 0.46698 0.46435 0.46481

1 0.06401 0.06337 0.06339 0.06337

1.5 0.02157 0.02149 0.02142 0.02138

2 0.01064 0.01062 0.01060 0.01055

In formulas (55) and (56), the second terms, which are determined the value of displacements due to shear deformations, are the same. Note that when calculating according to the proposed method, displacements from shear deformations are determined independently of the bending state and will also be the same for the considered support options of the plates. The results in the Tables show the general tendency for the displacement values to approach the exact values from above, as well as good accuracy even for coarse grids.

Table 5. The clamped round plate center displacement (Fig. 6). Shear locking test.

t, m w, mm

Fig. 6a Fig. 6b Fig. 6c Formula (56)

0.1 14.3394 14.0010 13.8079 13.8908

0.5 0.12821 0.12561 0.12415 0.12461

1 0.02130 0.02102 0.02086 0.02084

1.5 0.00891 0.00885 0.00882 0.00876

2 0.00530 0.005284 0.005280 0.00524

When modeling the round plate, arbitrary quadrangular elements were used, but most of them have a square shape. But the finite elements located along the plate border have different sizes and arbitrary shapes. These elements are located near the supports, therefore, the greatest transverse forces arise in them. They will make a significant contribution to the displacements caused by the shear, and thereby determine the calculation accuracy. Thus, we note that the proposed arbitrary quadrangular finite elements allow one consider shear deformations with high accuracy and regardless of the ratio of the plate thickness to its dimensions. The shear locking effect in the calculation of thin plates is absent.

4. Conclusion

1. An arbitrary quadrangular bending finite element based on a piecewise constant approximation of moments is proposed. The solution is based on the principles of minimum additional energy and possible displacements. The finite element allows you to consider the shear deformation, regardless of the ratio of the plate thickness to its sizes. The effect of locking in the calculation of thin plates is absent.

2. There is compared the results of the oblique plates calculations, annular and round plates with analytical solutions and calculation results for other programs. The results comparison shows good accuracy in determining displacements and moments. The crushing of finite elements mesh cause the displacement values tend to exact values from above.

3. To assess the influence of the finite element shape, calculations were made for skew plate. To model the plate, quadrangular elements were used, obtained from rectangular ones by changing the slope of one side. As a result, one of the finite element internal corners varied from 90 to 150 degrees. At the same time, the calculating accuracy of the displacements and moments decreased only by 4 and 6 percent, respectively.

4. The proposed finite element is easy to implement. A numerical integration and the mapping of the quadrangular region to the rectangular one is not used to get the solution. All necessary expressions are obtained analytically.

References

1. Gallagher, R.H. Finite Element Analysis: Fundamentals. Englewood Cliffs, Prentice-Hall, 1975. 420 p.

2. Eslami, M.R. Buckling and Postbuckling of Beams, Plates and Shells. Springer. New York, 2018. 588 p.

3. Anssi, T. Karttunen, Raimo von Hertzen, J.N. Reddy, Jani Romanoff. Exact elasticity-based finite element for circular plates. Computers & Structures. 2017. Vol. 182. Pp. 219-226.

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

5. Karttunen, A.T., Hertzen, R., Reddy, J.N., Romanoff, J. Shear deformable plate elements based on exact elasticity solution. Computers & Structures. 2018. Vol. 200, Pp. 21-31.

6. Fallah, N. On the use of shape functions in the cell centered finite volume formulation for plate bending analysis based on Mindlin-Reissner plate theory// Computers & Structures. 2006. Vol. 84. Pp. 1664-1672.

7. Khezri, M., Gharib, M., Rasmussen, K.J.R. A unified approach to meshless analysis of thin to moderately thick plates based on a shear-locking-free Mindlin theory formulation. Thin-Walled Structures. 2018. Pp. 161-179.

8. Katili, I., Batoz, J.L., Maknun, I.J., Hamdouni, A., Millet, O. The development of DKMQ plate bending element for thick to thin shell analysis based on the Naghdi/Reissner/Mindlin shell theory. Finite Elements in Analysis and Design. 2015. Vol. 100. Pp. 12-27.

9. Mirsaidov, M.M., Abdikarimov, R.A., Vatin, N.I., Zhgutov, V.M., Khodzhaev, D.A., Normuminov, B.A. Nonlinear parametric oscillations of viscoelastic plate of variable thickness. Magazine of Civil Engineering. 2018. 82(6). Pp. 112-126. DOI: 10.18720/MCE.82.11

10. Sukhoterin, M.V., Baryshnikov, S.O., Knysh, T.P., Abdikarimov, R.A. Natural oscillations of a rectangular plates with two adjacent edges clamped. Magazine of Civil Engineering. 2018. 82(6). Pp. 81-94. DOI: 10.18720/MCE.82.8.

11. Abdikarimov, R., Khudayarov, B., Vatin, N.I. To Calculation of Rectangular Plates on Periodic Oscillations. MATEC Web of Conferences. 2018. Vol. 245. 01003.

12. Kumara, R., Lala, A., Singhb, B.N., Singhc, J. New transverse shear deformation theory for bending analysis of FGM plate under patch load. Composite Structures. 2019. Vol. 208. Pp. 91-100.

13. Duan, H., Ma, J. Continuous finite element methods for Reissner-Mindlin plate problem. Acta Mathematica Scientia. 2018. Vol. 38. Pp. 450-470.

14. Park, M., Choi, D.-H. A two-variable first-order shear deformation theory considering in-plane rotation for bending, buckling and free vibration analysis of isotropic plates. Applied Mathematical Modelling. 2018. Vol. 61. Pp. 49-71.

15. Lukashevich, A.A. Modelling of contact interaction of structures with the base under dynamic loading. Magazine of Civil Engineering. 2019. 89(5). Pp. 167-178. DOI: 10.18720/MCE.89.14

16. Lin, M., Junping, W., Xie, Ye. A hybridized formulation for the weak Galerkin mixed finite element method. Journal of Computational and Applied Mathematics. 2016. Vol. 307. Pp. 335-345.

17. Rodrigues, J.D., Natarajan, S., Ferreira, A.J.M., Carrera, E., Bordas, S.P.A. Analysis of composite plates through cell-based smoothed finite element and 4-noded mixed interpolation of tensorial components techniques. Computers & Structures. 2014. Vol. 135. Pp. 83-87.

18. Doa, T.V., Nguyenb, D.K., Ducc, N.D., Doanc, D.H., Buie, T.Q. Analysis of bi-directional functionally graded plates by FEM and a new third-order shear deformation plate theory. Thin-Walled Structures. 2017. Vol. 119. Pp. 687-699.

19. Senjanovic, I., Vladimir, N., Hadzic, N. Modified Mindlin plate theory and shear locking free finite element formulation. Mechanics Research Communications. 2014. Vol. 55. Pp. 95-104.

20. Ullah, S., Zhong, Y., Zhang, J. Analytical buckling solutions of rectangular thin plates by straightforward generalized integral transform method. International Journal of Mechanical Sciences. 2019. Vol. 152. Pp. 535-544.

21. Tyukalov, Yu.Ya. Finite element models in stresses for bending plates. Magazine of Civil Engineering. 2018. 82(6). Pp. 170-190. DOI: 10.18720/MCE.82.16

22. Tyukalov, Yu.Ya. Finite element models in stresses for plane elasticity problems. Magazine of Civil Engineering. 2018. 77(1). Pp. 23-37. DOI: 10.18720/MCE.77.3

23. Tyukalov, Yu.Ya. Equilibrium finite elements for plane problems of the elasticity theory. Magazine of Civil Engineering. 2019. 91(7). Pp. 80-97. DOI: 10.18720/MCE.91.8

24. Tyukalov, Yu.Ya. Calculation method of bending plates with assuming shear deformations. Magazine of Civil Engineering. 2019. 85(1). Pp. 107-122. DOI: 10.18720/MCE.85.9

25. Tyukalov, Yu.Ya. Finite element model of Reisner's plates in stresses. Magazine of Civil Engineering. 2019. 89(5). Pp. 61-78. DOI: 10.18720/MCE.89.6

26. Timoshenko, S.P., Voynovskiy-Kriger, S. Plastiny i obolochki [Plates and shells]. Moskow: Nauka, 1966. 636 p. (rus)

27. Belkin, A.Ye., Gavryushin, S.S. Raschet plastin metodom konechnykh elementov [Finite element analysis of plates]. Moskow: MGTU, 2008. 231 p. (rus)

Contacts:

Yury Tyukalov, [email protected]