CC BY
9Magazine of Civil Engineering. 2022. 110(2). Article No. 11004
Magazine of Civil Engineering
journal homepage: http://engstroy.spbstu.ru/
ISSN 2712-8172
DOI: 10.34910/MCE.110.4
Finite element models based on the approximation of discontinuous stress fields
A.A. Lukashevich
St. Petersburg State University of Architecture and Civil Engineering, St. Petersburg, Russia E-mail: aaluk@bk.ru
Keywords: finite element method, plane elasticity problem, discontinuous stress approximation, penalty function method, functional of additional energy, Reissner functional
Abstract. The paper develops the finite element method (FEM) in the form of the force method and in the mixed form for the calculation of structures. At present, the displacement-based finite element method is mainly used for engineering calculations. Stress-based and mixed finite element formulations are not so widely spread, but in some cases these formulations can be more effective in particular with respect to calculating stresses and also obtaining a two-sided estimate of the exact solution of the problem. The finite element models based on the approximation of discontinuous stress fields and the use of the penalty function method to satisfy the equilibrium equations are considered. It is shown that the continuity of both normal and tangential stresses only on the adjacent sides of the finite elements contributes to the expansion of the class of statically admissible stress fields. At the same time, the consistent approximation is provided, both of the main part of the functional of additional energy, and its penalty part. The necessary matrix relations for rectangular and triangular finite elements are obtained. The effectiveness of the developed models is illustrated by numerical studies. The calculation results were compared with the solution on the FEM in displacements, as well as with the results obtained using other schemes of approximating the stresses in the finite element. It is shown that the model of discontinuous stress approximations gives the bottom convergence of the solution, both in stresses and in displacements. At the same time, the accuracy on the stresses here is much higher than in the displacement-based FEM or when using conventional stress approximation schemes.
One of the important steps in the design process of building structures and constructions is to improve the method of their calculation. This problem is very relevant, since it helps to give valid results, which ultimately reduces the efforts and costs in the construction process. Nowadays, the calculation of structures is carried out, as a rule, using the finite element method (FEM). Many fundamental works have been devoted to the development of the theoretical principles of the FEM and questions of its application in structural mechanics [1-3]. They consider the basic variational principles and the corresponding FEM formulations, on the basis of which the finite element models can be constructed to solve different problems.
There are three main forms of FEM, each form is an analogue of one of the three classical methods of structural mechanics of rod systems - the displacement method, the force method, and the mixed method. The most widely used engineering calculations method is the FEM in the displacement method form. The mathematical formulation of the problem here is based on the variational Lagrange principle, i.e., the minimum principle of the total potential energy of the system. The main unknowns in this case are the displacements of the nodal points of the discrete model of the structure. The stresses are secondary and can be calculated by numerically differentiating the displacements. The finite element method in displacements is widely used to solve geometrically and physically nonlinear problems [4-6], constructively
Lukashevich, A.A. Finite element models based on the approximation of discontinuous stress fields. Magazine of Civil Engineering. 2022. 110(2). Article No. 11004. DOI: 10.34910/MCE.110.4
© Lukashevich, A.A., 2022. Published by Peter the Great St. Petersburg Polytechnic University. This work is licensed under a CC BY-NC 4.0
1. Introduction
nonlinear problems with unilateral constraints [7-10], problems of stability and dynamics of structures [11-14].
The advantages and disadvantages of the displacement method are well known. The huge advantage of this FEM form is its complete formalization (and, accordingly, the ease of implementation to software), as well as good stability of the solution with guaranteed convergence to the lower boundary. However, the accuracy of stresses (forces) determination is much lower than the displacements, although it is stress values which are more important in the structural strength analysis. In this connection, special refinement algorithms are often used to calculate stresses, e.g. [15]. In addition, since the approximate solution in displacements corresponds to the lower boundary, the values of both displacements and stresses are underestimated relative to the exact values.
Attempts to overcome these shortcomings, based on the use of FEM schemes directly in stresses (forces) or in a mixed form, have been made repeatedly, but this problem is still far from complete and remains one of the most important problems of the finite element method application in structural mechanics.
Castiliano's minimum principle of additional energy and the related FEM schemes in the form of the force method, as well as the Reissner variational principle (mixed method), have not received such a wide application. This situation is caused by a number of circumstances, in particular, the need to satisfy the equilibrium equations in the force method or an increase in the number of unknowns in the mixed method. However, in some cases, these approaches can be more effective, especially with regard to the accuracy of calculating stresses. In addition, performing dual calculations based on the alternative forms of FEM allows, as a rule, to obtain a two-side estimate of the exact solution of the problem [1, 2].
The FEM schemes, in which the search of the solution is based on the approximation of stresses (forces) in the finite element region, were considered in [1, 2, 16-22]. In [19], the combination of principles of possible displacements and possible changes in the stress state is used to find the solution. The solution to the problems of the elasticity theory in stresses on the basis of the functional of additional energy is considered in [21, 22]. In this case, the principle of possible displacements and Lagrange multipliers are used to satisfy the equilibrium equations.
The main advantage of the FEM in the form of the force method is that the main unknowns here are stresses. And if there were no certain difficulties in implementing the force method [1, 3], the stress values could be obtained with the same degree of accuracy as the displacements in the FEM scheme of the same name. In addition, the use of Castiliano's principle gives the upper boundary of the approximate solution (i.e., the stresses are overestimated), which, in principle, is better for strength calculations than the underestimated estimate. On the other hand, solving the problem in stresses can supplement the usual calculation of FEM in displacements, including from the point of view of the two-sided estimation of the solution of the problem on energy. Nevertheless, there are no algorithms that are equally as simple and stable, with guaranteed convergence in an extensive class of problems, similar to the FEM in the form of the displacement method.
The variational formulation of the mixed method is based on the principle of the stationarity of various forms of the Reissner functional. Some variants of the FEM in a mixed formulation were considered in [2, 3, 23-31]. With this approach, displacements and stresses within each finite element are approximated simultaneously, therefore, there is no need to overestimate the requirements for the continuity of the desired functions and their derivatives. On the contrary, it is possible to set the necessary approximations, and since the mixed variational principles lead to a mixed form of relations between stresses and displacements for the finite element, more reliable solution can be obtained.
On the other hand, the Reissner functional is not convex; its surface at the point of stationary has a prominent saddle-like behavior. The system of resolving equations corresponding to the formulation of a mixed type is not a symmetric and positive definite. Therefore, the approximate solution obtained by the FEM in the mixed form is characterized by some imbalance in the fulfillment of the equilibrium conditions and the conditions for the compatibility of deformations and, when crushing the mesh, it can approach the exact solution both from below and above. These circumstances to some extent make the direct use of the Reissner functional in the FEM more difficult [1, 3].
This paper is devoted to the development of alternative forms of the FEM based on the force and mixed methods. It is proposed to use the approximation of discontinuous stress fields in the finite element region to construct a solution to the plane problem of the elasticity theory. Such an approach may have advantages in terms of more accurately fulfilling the equilibrium conditions, as well as in solving contact problems and some problems related to stress concentration.
2. Methods
Let us consider the following version of the FEM in the force method form, which allows us to efficiently solve a wide class of problems directly in stresses. It is based on the approximation of discontinuous stress fields and the use of the penalty method to satisfy equilibrium equations. The corresponding variational formulation for the plane problem of the elasticity theory has the form of a modified functional of additional energy [32, 33]:
n (p) = -1 j>}r [D]-1{ p}dQ + J {u} [Ls ]T {p}dS + a j( [ A]T{p} + {p]j ([ A]T{p} + {p}) dQ (1)
•a
Q
under the additional condition at the region boundary: [Ls]T{a} = {gs} e S
g '
Here { p} = { px py Txy }T is the stress vector; [D] is the stiffness coefficients matrix; [A] is the differentiation operations matrix in the equilibrium equations; [ Ls ] is the direction cosine matrix of the external normal to the boundary S = Sg + Su of the plane region Q ; {gs} is the surface force vector at
the boundary Sg ; {us } is the vector of given displacements at the boundary Su ; {p} = { px py }T is the volumetric force vector.
When solving problems in stresses, the main restrictions on the smoothness of the desired functions are imposed by the equilibrium equations - it is necessary to ensure the existence of piecewise continuous derivatives of the stress components. Therefore, according to them, for the plane problem of the elasticity
theory, the differentiability of the normal stresses px is only on x, py is only on y, and tangential stresses Txy is both on x and on y is required. Thus, normal stresses can undergo discontinuities at sites
perpendicular to the element boundaries. In comparison with the continuous approximation, the use of such an assumption contributes to the expansion of the class of statically allowable stress fields, among which the solution is sought. This allows to minimize the additional energy functional to a greater extent - as a result, the solution of the problem is, on average, closer to the exact one. The use of discontinuous approximations leads to the need to apply a special class of finite elements. The location of the nodal points here should ensure that the conditions of discontinuity of normal stresses are satisfied.
Figure 1 shows one of the simplest elements of this kind, which is a rectangle, the four nodes at the vertices of which correspond to tangential stresses, two nodes on each of the sides correspond to normal
stresses px, two nodes on the upper and lower sides correspond to stresses py .
Figure 1. Eight-node finite element with discontinuous stress approximation.
Here the stress distribution is given by the following approximating polynomials:
px = ax + a2 x; py = a3 + a4y; TXy = a5 + a6 x + a7y + a8 xy. Substituting the approximation data in the equilibrium equations, we obtain:
a2 + a7 + a8x + px = 0; a4 + a6 + a8y + py = 0.
(2)
(3)
Obviously, for the constant volumetric element forces px, Py the conditions (3) can be satisfied only if the following condition is identically satisfied:
a8 x = a8 y = 0. (4)
Since this is not provided by the exact integration of the penalty term of functional in Eq. (1), a method for artificial lowering the accuracy of calculating the integrals was proposed for the reduced finite element, which facilitates the zeroing of the expression in the integral part of the penalty term [33]. Despite the fact that this finite element gives the acceptable accuracy and convergence of the results of solving problems in stresses, there are certain disadvantages. First, to calculate the coefficients of the deformability matrix related to the penalty term, it is necessary to apply the procedure of numerical integration. Secondly, for the main part of functional (1) and its penalty term, different stress approximations are used - for the penalty term, they are respectively of a lower order than for the main part of the functional. The model proposed here for approximating discontinuous stress fields for a finite element of a plane problem of the elasticity theory is free from these disadvantages.
Consider the above condition (4). For the finite element (Figure 1) and approximation (2) described above, exact execution of (4) is possible only when using the single-point integration scheme of the penalty term, i.e. when the integration order is lowered. However, locating nodes of the finite element on the x, y
axes (Figure 2(a)), the condition (4) is reduced to the equality a8x = a8y.
Substituting this equality into approximation (2), we obtain the following variants of polynomials for
xy ■
2 2 t = a5 + a6 x + a7 y + a8 x ; t = a5 + a6 x + a7 y + a8 y .
(5)
For each of the polynomials (5), we define the form functions for nodes and average them. Then the nodal stress vector for the finite element
{= {°1x T1xy °2y t2 xy °3x T3xy ° 4y t4 xy f
the following form functions will be a respond:
1 - x 1 _ x x
N13x = —+—; N13 =—+—+—-1,3x 2 2a 'xy 4 2a 4a2
y2 at 1 - y •
N2,4y=2+^'
N
4b
2,4 xy
1_ y
= — + — -
(6)
22 xy
+ . (7)
4 ' 2b 4a2 4b
2
The approximation of the tangential stresses within a given finite element, in particular, for the shape function Nlxy, is shown in Figure 2(b).
Figure 2. Four-node finite element with discontinuous stress approximation: (a) location of the nodes; (b) approximation of the tangential stresses.
Expressing the stresses using the form functions (7) through the nodal values (6) and substituting them into stationarity condition of the functional (1) for finite element, we obtain the matrix equation for the
vector {a)e:
[Ka]e {= {U}e
(8)
The deformability matrix of the finite element can be represented as the sum of two matrixes [Ka]e = [K0]e + [Ka]e . The first term [K0]e = J[N.]T [D]-1 [N.]dQ corresponds the main part of
thefuncona, (1), the aecond [Kar = -a j([* [N..[* [N.] dQ „ .s ^ pa. The righ Part
is the vector of given nodal displacements
[ Na ] =
a
ts {U}e = J[Na]T [Ls]-1K}dS + a j([^]T [Na])T{p}
da.
Su
Nix 0 0 0 Nsx 0 0 0
0 0 N2 y 0 0 0 N4 y 0
0 N1xy 0 N2 xy 0 N3xy 0 N4 xy
a
is the matrix of form functions for
stresses in finite element.
Perfoming the matrix operations and integrating the obtained expressions, we obtain:
[ Ko]e = k
[Ka ]e = a
"113^ 0 37 n 0 -7 n 0 37 n 0 "
0 240 0- -180v 0 120 0 -180v
37 n 0 113^ 0 37 n 0 -7 n 0
0 -180v 0 240 0- -180v 0 120
-7 n 0 37 n 0 113^ 0 37 n 0
0 120 0- -180v 0 240 0 -180v
37 n 0 -7 n 0 37 n 0 113 n 0
_ 0 -180v 0 120 0- -180v 0 240 _
4m + n 0 - m - n -1 n - 2m 0 - m - n 1
0 3m -1 0 0 - 3m 1 0
- m - n -1 m + 4n 0 - m - n 1 m - 2n 0
-1 0 0 3n 1 0 0 - 3n
n - 2 m 0 - m - n 1 4m + n 0 - m - n -1
0 - 3m 1 0 0 3m -1 0
- m - n 1 m - 2n 0 - m - n -1 m + 4n 0
1 0 0 - 3n -1 0 0 3n
ab „ N a b
Here: k =-; u = 2(1 - v); m = — ; n =—, where E, v are the elastic modulus and
180E 3b 3a
Poisson's ratio; a, b are the size of the rectangular finite element.
The penalty parameter a has the physical meaning of the deformability coefficient of the weakest possible elastic base. Its value should be large enough and limited from above only by the accuracy of calculations. The recommendations for choosing the values a for different problems are given in [21, 32,
34]. As applied to the plane problem of the elasticity theory, a = (107^ 109)k can be accepted.
Similarly, the approximation of discontinuous stress fields can be applied to construct finite elements in stresses of other geometric shapes, as well as extended to finite elements of the mixed form. In particular, for a finite element in the form of a right triangle (Figure 3(a)), the stress distribution can be represented by the following polynomials:
ox = ai + x; o y = a3 + a4 y; r^y = a5 + a6 x + a7 y.
(9)
It is clear that the fulfillment of the condition similar to (4) is not required here. The following form functions will respond to the vector of nodal stresses of a finite element
{°} - {°1x T1xy °2y T2xy °3x °3y T3xy} :
Nlx = --; Nlxy = -X; Nly = -y; Nlxy = -{; = 1 + X N3y = 1 + {; N^ = 1 + X + £. (10) a a b b a b a b
We replace the nodal stress components o3y, o3x, t^^ with the tangent and normal components to the inclined face of the element o^n, T^n (Figure 3(b)). The nodal stress vector will then be
{o}e = {o1x T1xy o2y T2xy o3n T3n Y .
Figure 3. Triangular finite element with discontinuous stress approximation.
After substituting the form functions (10) and converting the stress components during turn of the axes, we obtain the following terms that make up the deformability matrix of a triangular finite element:
[^ole -
X 0 0 0 0 0
0 1 vlm -vl2 0 0
0 vlm X(l2 -m2) lm(2 X+2+v) 0 -lm(2+v)
0 -vl2 lm(2 X+2+v) 2(m2 -vl2) 0 -vm2
0 0 0 0 X 0
0 0 -lm(2+v) -vm2 0 1
[Ka ]e - a
n
0
lm-nx
,2
0 n
-n-lm-x
-l -2n-lm -2lm-n-m2 1 0 -1 0
lm-nx
-n- lm- x
x-P
X+n -lm+rn+2 P -lm -lm- m H-lm-x
-l 2 -2n lm -2lm-n m2
x+n-lm+rn+2 P-lm -lm-nx
-1
0
V-lm-x
n - m2 +4lm -m2 -2 ¡a. - lm -2lm-¡a -12
-m2-2¡a-lm -2lm-^-l
A A
- A A
Here: k —
2ab
2-
; X = 2(1 + v) ; q = a; /u = b; ft = n + n ; x = (l2 - , where E, v are the 3E b a
elastic modulus and Poisson's ratio; a, b are the dimensions of the legs of the triangular finite element; l,
m are cosine and sine to the normal of the inclined face.
Let us consider a mixed finite element model that allows to solve effectively the problems of the elasticity theory directly both in displacements and stresses. It is also based on the approximation of discontinuous stress fields and using the penalty method to satisfy equilibrium equations. The variational formulation of the mixed problem corresponds to the principle of stationarity of various forms of the Reissner functional, which directly includes the components of both displacements and stresses as well. As already indicated, this functional is not convex; its surface at the stationary point has a prominent saddle-like shape. This circumstance significantly complicates its use in the finite element method (the matrix of coefficients in this case is not symmetric and positive definite).
K
A convex mixed functional can be obtained by subtracting the Lagrange functional from the first form of the Reissner functional [33]:
n(u,o) = nRi (u,o) - nL(u) = J{o}T[A]{u}dQ - 1 J{o}T [D]-1{o}dQ -
a 2 a
2 a s,
, (11)
- - J([A]{u])T [D][A]{u}da - J({u}-{us}J [Ls]T{a}dS.
u
T ,
under the additional conditions: [A]T { a} + {p} = 0 e a , [Ls ]T { a} = {gs } e Sg .
The convexity of the functional (11) is obtained by moving the equilibrium conditions from natural (for the functional nr (u, a) to additional ones. We use the penalty method in order to satisfy these conditions:
n(u,a) = J{a}T[A]{u}da - - J{a}T [D]-1{a}da - 1 J([A]{u})T [D][A]{u}da +
a 2 a 2 a (12)
f J ([A]T { a} + { p})T ([A]T { a} + {p})da - J ({ u} - { us}f [Ls ]T { a} dS.
+ a J \^A] { o} + { p}) ^[A]{ o} + ^ isy, ^s
a Su
with the additional condition [Ls]T{o} = {^} e Sg .
It is easy to see that the reduced functional is equivalent to the variational statement of the problem in a least squares form. This circumstance provides the convexity condition for the functional (12), which makes it more convenient for applying the finite element method. Thus, setting an approximating expression
for displacements {u} = [Nu]{u}e and, accordingly, for stresses {o} = [No]{o}e , and varying the
functional (12) in a discrete form, we obtain the following matrix equilibrium equation for a finite element with a symmetric, positive definite coefficient matrix:
[Ku ] [Kua ] [Kau ]e [Ka ]'
xJ{u}e[ = J{R}Ê >, (13)
!{ a}e J [{^ }ef
where [Ku ]e = J( [A] [Nu ] ) T [D][A][Nu ] da and [Ka ]e = J [N. ]T [D]-1[N. ] da -
a a
- a j([A]T [No] )T[A]T [No] dQ are the stiffness and deformabilify mat^xes, respectively;
Q
[Kuo]e =-J([A][Nu])T[No]dQ and [Km]e =([Kuo]e)T are the mixed matrixes; ], [No] are
a
the form function matrixes for displacements and stresses, respectively; {u}e, {o}e are the vectors of
nodal generalized displacements and stresses, respectively; {R}e, {U}e are the vectors of generalized reactions and displacements in the nodes of the finite element, respectively.
The approximation of displacements within a finite element (Figure 4) can be similar to that usually used in the FEM in displacement method form. For stresses, however, it is more efficient to use the discontinuous approximations, such, as those discussed above, e.g. (2).
Figure 4. Eight-node finite element with discontinuous stress approximation.
The FEM solutions considered above using discontinuous stress approximations are implemented in a computer program. With its help, numerical studies and a comparative analysis of the finite element solutions based on the formulation of the displacement method, the force method, and the mixed method were performed.
3. Results and Discussion
Let us consider the problem of unilateral contact of a three-layer array with a rigid base (Figure 5). The upper and lower layers are concrete (elasticity modulus Ei = 1.9-104 MPa, v = 0.18), the middle layer is brickwork (E2 = 0.27-104 MPa, v = 0.16). Under the action of a distributed load applied at the left end of the calculated region, the zone of separation of the lower surface of the array from the rigid base is formed. The dimensions of the contact and separation zones are calculated using the standard iterative algorithm for unilateral constraints.
Figure 5. Unilateral contact of the three-layer array with the rigid base.
The solution of the problem based on discontinuous stress approximations was compared with the traditional FEM scheme in the displacement method form, as well as with the relatively accurate solution obtained with a sufficiently dense finite element mesh using the LIRA-SAPR software package [35]. Comparison of the proposed approach with the traditional FEM was carried out under the same discretization of the computational region, i.e. 18 x 6 elements.
The stresses ax along the vertical cross-section at a distance of 0.4 m from the left face of the array are compared in Figure 6(a) and the stresses ay along the lower face of the contact zone of the array with the rigid base are compared in Figure 6(b). The solid line shows the exact solution, the shaded circles -the FEM solution in the displacement method form, the not shaded circles - in the force method form.
-100 0 100 200 300 crx, kN/m2 0 0.2 0.4 0.6 0.8 x, m
Figure 6. Stress distribution in characteristic sections of the three-layer array.
The graphs show that the stresses obtained by the FEM in the form of the force method are closer to the exact values, than the form of the displacement method gives. This is especially manifested at the region edges under consideration, as well as near the interfaces of different-modular materials.
In order to compare various forms of the finite element method, the problem of determining the stressstrain state of a cantilever plate was solved (Figure 7(a)). The Poisson's ratio was assumed to be 0.3, the remaining values: the elastic modulus E, the plate thickness h, and the load p were expressed in general form (in computer calculations they were set equal to one). Meshes of different densities were used, e.g. Figure 7(b) shows the 2*3 mesh when calculating the plate by the mixed method.
Figure 7. The problem of the cantilever plate.
The stresses and displacements obtained by the considered approaches were compared with the results of the FEM calculation in the displacement method form. Table 1 shows (in the terms of E, h, p) the vertical displacements v of the middle of the free edge (x = 12 m, y = 0) and the stress ax at the upper point of the fixed edge of the plate (x = 0, y = 3 m) for different finite element meshes. The stresses at the extreme points here are obtained using linear interpolation.
Table 1. Comparison of the plate calculation results with different FEM forms.
Displacement method
Force method
Mixed method
Mesh Number of JP 1 Number of rxfP 1 Number of JP1
unknowns 1 Eh J 1 h J unknowns 1 h J unknowns 1 Eh J 1 h J
2x2 12 152.5 6.67 24 16.63 42 156.2 15.97
3x2 18 156.2 6.83 34 16.20 58 183.8 15.74
3x3 24 184.4 8.94 48 14.27 80 189.4 14.39
4x3 40 200.8 10.34 62 14.42 102 197.6 14.52
5x4 70 213.1 11.70 98 14.56 158 211.7 14.72
6x4 96 216.6 12.19 116 14.90 186 215.9 14.81
Exact solution 225.0 15.10 15.10 225.0 15.10
The calculation results for the given variants of the finite element meshes show the following. The FEM in the force method form provides a monotonic convergence of the solution in stresses to exact values from below. The stresses obtained by the mixed method approach the exact solution at first (for coarse meshes) from above, then (for fine meshes) also from below. Moreover, the use of discontinuous approximations of stress fields allows to achieve better convergence and accuracy of the solution in comparison with the displacement method with approximately the same order of the resolving equations system.
As it can be seen, on coarse meshes, the stresses obtained by the FEM in the force method form are significantly closer to the exact values than those obtained by the FEM in displacements. Under crushing the mesh (with the equivalent dimension of the problem), the accuracy of determining the stress by the force method and the mixed method also remains higher than by the displacement method. So, with the number of unknowns of the order of 100, the values of the stresses obtained by the force method differ from the exact values by 3.6 %, the mixed method by 3.8 %, and the displacement method by 19.3 %. In turn, for the identical meshes, the mixed FEM allows one to obtain stresses with the same accuracy as FEM in stresses, and displacements - as FEM in displacements.
Figure 8 illustrates the convergence of solutions obtained using various forms of the FEM. Obviously, according to the accuracy of calculating the stresses, the approaches considered here are more effective in comparison with the traditional FEM scheme in displacements.
■te)
200 150 100
Exact solution
: method
v Displace nent methc d
Number of unknowns
15 -\ >
12 ■■■
9 ■■■
Exact solution
^ Force t method ^ lixed methc
Dlspla cement me thod
Number of unknowns
40
80
120
160
40
80
120
160
Figure 8. The convergence of solutions obtained by the different FEM forms.
The problem of bending a cantilever plate under the action of a vertical distributed load applied to its free end is presented below in Figure 9(a). The plate thickness is 1 m, the elastic modulus is 10000 kN/m2, the Poisson's ratio is 0.25. This example was given in [21], where the application of the finite element models in stresses at constant and piecewise-constant approximations of stresses in the finite element region was considered. Below, we compare the solutions constructed on the basis of approximations of discontinuous stress fields with the results obtained in [21] using conventional schemes of stress approximation in the finite element.
Table 2 shows the stress values in the upper fibers of the section at the clamp (x = 0, y = 1 m), obtained using different stress approximation types in the finite element. The calculation results are presented for four variants of the finite element meshes. For comparison, the lower row of the table shows the stress values obtained by FEM in the LIRA-SAPR software package for a relatively dense mesh.
Table 2. Comparison of the plate calculation results for different types of stress approximation.
Constant stresses Piecewise constant stresses Discontinuous stresses
Mesh Ox, kN/m2 kN/m2 тíy, kN/m2 Ox, kN/m2 kN/m2 тíy, kN/m2 Ox, kN/m2 kN/m2 тXy, kN/m2
12x4 13.96 0.95 1.11 17.56 3.12 1.75 22.18 3.42 2.75
24x8 16.53 1.43 1.45 19.55 3.52 2.40 23.25 3.67 3.08
48x16 18.91 1.78 1.93 22.06 4.03 3.22 23.62 3.80 3.44
96x32 21.45 2.10 2.49 25.16 4.65 4.06 24.17 3.97 3.69
192x64
(LIRA- 24.02 3.46 3.62 24.02 3.46 3.62 24.02 3.46 3.62
SAPR)
The Table shows that even on the coarse meshes, the stress values obtained on the basis of discontinuous approximations are more accurate compared to using the constant and piecewise constant
stress approximations. So, for the mesh 12*4, the stress values ox differ upward from the constant and
piecewise-constant stresses by 37 % and 21 %, respectively. The stresses oy and zxy in our case are
also much more accurate than with ordinary schemes of approximations.
When the mesh is crushed, the model of discontinuous approximations gives a monotonous and stable convergence of stresses to the lower boundary. In this case, the numerical stress values are in the gap between the stress values obtained with the constant and piecewise-constant approximations. For
example, for mesh 96*32, the stress ox is 11% higher than at constant stresses in the element, but 4 %
less than for the piecewise-constant stresses variant. Figures 9(b-d) illustrate the convergence of solutions in stresses for different types of stress approximations.
0 1000 2000 3000 0 1000 2000 3000
Figure 9. The convergence of solutions in stresses for different types of approximation.
4. Conclusions
Summarizing the results of numerical researches, we can draw the following conclusions:
1. The finite element method in the form of the force method in the general case gives the upper convergence of the solution in stresses. This allows us to obtain a two-sided estimate of the exact solution of the corresponding problem. At the same time, the finite-element solutions built on the stress approximation are much closer to exact values than those obtained in the displacement method formulation.
2. The use of discontinuous stress approximations contributes to the expansion of the class of statically admissible stress fields. In this case, the same approximation is provided for both the main part of the additional energy functional and its penalty part. This allows us to obtain more accurate stress values as compared to the continuous stress approximation.
3. The mixed finite-element models allow to obtain alternative solutions to the considered problems, thereby ensuring their greater validity and reliability. Under identical meshes, the mixed finite element method allows us to obtain stresses with the same accuracy as the FEM in stresses, and displacements -as the FEM in the displacement method form.
4. Discontinuous stress approximations give a fast and stable convergence of the resulting solution to the lower boundary. At the same time, the accuracy of the stress in this case is significantly higher than when using conventional stress approximation schemes, in particular, constant and piecewise constant approximations.
References
1. Gallagher, R.H. Finite element analysis: Fundamentals. Englewood Cliffs: Prentice-Hall, 1975. 416 p.
2. Zienkiewicz, O.C., Taylor, R.L., Zhu, J.Z. The Finite Element Method: Its Basis and Fundamentals. Butterworth-Heinemann, 2013. 756 p.
3. Bathe, K.J. Finite element procedures. Englewood Cliffs: Prentice-Hall, 2016. 1043 p.
4. Lalin, V.V., Rybakov, V.A., Morozov, S.A. The finite elements research for calculation of thin-walled bar systems. Magazine of Civil Engineering. 2012. No. 27(1). Pp. 53-73. (rus). DOI: 10.5862/MCE.27.7.
5. Paola, M., Scimemi, G. Finite element method on fractional visco-elastic frames. Computers & Structures. 2016. Vol. 164. Pp. 1522. DOI: 10.1016/j.com pstruc.2015.10.008.
6. Lukashevich, A.A., Lukashevich, N.K., Timohina, E.I. Modelling of contact interaction with allowance for nonlinear compliance in unilateral constraints. IOP Conference Series: Materials Science and Engineering. 2018. Vol. 463. 042054. DOI: 10.1088/1757-899X/463/4/042054.
7. Wriggers, P., Schroder, J., Schwarz, A. A finite element method for contact using a third medium. Computational Mechanics. 2013. 52(4). Pp. 837-847. DOI: 10.1007/s00466-013-0848-5.
8. Lukashevich, A.A. Computational solution of contact problems with unilateral constraints and friction by step-by-step analysis. Advanced Materials Research. 2014. Vol. 941-944. Pp. 2264-2267. DOI: 10.4028/www.scientific.net/AMR.941-944.2264.
9. Wriggers, P., Rust, W.T., Reddy, B.D. A virtual element method for contact. Computational Mechanics. 2016. 58(6). Pp. 1039-1050. DOI: 10.1007/s00466-016-1331-x.
10. Lukashevich, A.A. Computational modelling of stiffness and strength properties of the contact seam. Magazine of Civil Engineering. 2018. 81(5). Pp. 149-159. DOI: 10.18720/MCE.81.15.
11. Rybakov, V.A., Lalin, V.V., Ivanov, S.S., Azarov, A.A. Coordinate functions quadratic approximation in Slivker's semi-shear stability theory. Magazine of Civil Engineering. 2019. 89(5). Pp. 115-128. DOI: 10.18720/MCE.89.10.
12. 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.
13. Alshurafa, S., Hanan, H., Polyzois, D. Finite element method for the static and dynamic analysis of FRP guyed tower. Journal of Computational Design and Engineering. 2019. Vol. 6. No. 3. Pp. 436-446. DOI: 10.1016/j.jcde.2018.08.004.
14. Lukashevich, A.A., Lukashevich, N.K. Numerical solution of contact problems with friction under dynamic loads. IOP Conference Series: Materials Science and Engineering. 2020. Vol. 753. 022058. DOI: 10.1088/1757-899X/753/2/022058.
15. 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. DOI: 10.1016/0045-7825(77)90005-6.
16. Muller, B., Starke, G. Stress-based finite element methods in linear and nonlinear solid mechanics. Advanced Finite Element Technologies. 2016. Pp. 69-104. DOI: 10.1007/978-3-319-31925-4_4.
17. Kuss, F., Lebon, F. Stress based finite element methods for solving contact problems: comparisons between various solution methods. Advances in Engineering Software. 2009. Vol. 40. 8. Pp. 697-706. DOI: 10.1016/j.advengsoft.2008.11.013.
18. Kuo, Y.-L. Stress-based finite element analysis of sliding beams. Appl. Math. Inf. Sci. 2015. Vol. 9. 2. Pp. 609-616. DOI: 10.12785/amis/092L36.
19. Negrozov, O.A., Akimov, P.A. Combined application of finite element method and discrete-continual finite element method. AIP Conference Proceedings 1800. 2017. Vol. 1800. 040014. DOI: 10.1063/1.4973055.
20. Peng, Y., Bai, Y., Guo, Q. Analysis of plane frame structure using base force element method. Structural Engineering and Mechanics. 2017. 62(1). Pp. 11-20. DOI: 10.12989/sem.2017.62.1.011.
21. 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.
22. 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.
23. Wriggers, P., Carstensen, C. Mixed finite element technologies. Springer verlag GMBH. 2009. 206 p. DOI: 10.1007/978-3-211-99094-0.
24. Qiu, W., Demkowicz, L. Mixed hp-finite element method for linear elasticity with weakly imposed symmetry: stability analysis // SIAM J. Numer. Anal. (SIAM Journal on Numerical Analysis) 2011. Vol. 49. 2. Pp. 619-641. DOI: 10.1137/100797539.
25. Stupishin, L.U., Nikitin, K.E. Mixed finite element for geometrically nonlinear orthotropic shallow shells of revolution // Advanced Materials Research. 2014. Vol. 919-921. Pp. 1299-1302. DOI: 10.4028/www.scientific.net/AMR.919-921.1299.
26. 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. DOI: 10.1016/j.cam.2016.01.004.
27. 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. 1. Pp. 26-54. DOI 10.1002/nme.5201.
28. Ignatyev, A.V., Ignatyev, V.A. On the efficiency of the finite element method in the form of the classical mixed method. Procedia Engineering. 2016. Vol. 150. Pp. 1760-1765. DOI: 10.1016/j.proeng.2016.07.167.
29. Ignatyev, A.V., Ignatyev, V.A., Gamzatova, E.A. Analysis of thin plates with excluding the displacements of the finite element as an absolutely rigid body by the fem in the form of a classical mixed method. Izvestiya vuzov. Stroitelstvo. 2018. 711(3). Pp. 5-13. (rus)
30. Gaur, H. A New Stress Based Approach for Nonlinear Finite Element Analysis. Journal of Applied and Computational Mechanics. 2019. Vol. 5. 3. Pp. 563-576. DOI: 10.22055/JACM.2019.29045.1548.
31. Lalin, V.V., Rybakov, V.A., Ivanov, S.S., Azarov, A.A. Mixed finite-element method in Slivker's semi-shear thin-walled bar theory. Magazine of Civil Engineering. 2019. 89(5). Pp. 79-93. DOI: 10.18720/MCE.89.7.
32. Zienkiewicz, O.C. Constrained variational principles and penalty function methods in finite element analysis. Conference on the Numerical Solution of Differential Equations. Lecture Notes in Mathematics, 1974. Vol. 363. Pp. 207-214. DOI: 10.1007/bfb0069138.
33. Lukashevich, A.A. Sovremennyye chislennyye metody stroitelnoy mekhaniki [Modern numerical methods of structural mechanics]. Izd-vo KhGTU. Khabarovsk, 2003. 135 p. (rus)
34. Taylor, R.L., Zienkiewicz, O.C. Complementary energy with penalty functions in finite element analysis. Energy methods in finite element analysis. Chichester: Wiley, 1979. Pp. 153-173.
35. Vodopyanov, R.Yu., Titok, V.P., Artamonova, A.E., Romashkina, M.A. LIRA-SAPR software package. User manual. Educational examples. Edited by Gorodetsky, A.S. Electronic edition, 2017. 535 p. https://www.liraland.ru/files/lira2017/format-pdf/.
Contacts:
Anatoliy Lukashevich, aaluk@bk.ru
Received 31.05.2020. Approved after reviewing 05.04.2021. Accepted 05.04.2021.
