close

Вход

Забыли?

вход по аккаунту

?

894

код для вставкиСкачать
INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING,
VOL.39,2337-2361 (1996)
HIGH-ACCURACY LOW-ORDER THREE-DIMENSIONAL
BRICK ELEMENTS
SHMUEL L. WEISSMAN
Symplectic Engineering Corporation, 1350 Solano Ave. #26, Albany, CA 94706, U.S.A.
SUMMARY
This paper presents a family of low-order high-accuracy three-dimensional brick elements.The elements are
formulated via a three-field variational principle. The assumed (independent)strain field is constructed from
two disjointed distributions. The first contains the lower-order distribution and its dimension is the
minimum required to satisfy stability requirements. An energy constraint, which is enforced weakly at the
element level, is used to relate the second distribution to the first. The stress field is chosen to a priori satisfy
a similar energy constraint. As a result, internal constraints (e.g. incompressibility) are automatically
satisfied by these fields, and locking behaviour is avoided. A J,-plasticity model illustrates the proposed
elements' performance in nonlinear solids. The excellent performance of the proposed elements is demonstrated with numerous challenging examples, including many that are usually modelled by shell elements.
KEY WORDS
elastoplastic; high-accuracy; locking-free;low-order; nonlinear; three-field functional
1. INTRODUCTION
This work presents a family of low-order high-accuracy eight-node brick elements that are
suitable for elastoplastic analysis (including the case of perfect plasticity). The importance of this
work stems from the following:
1. The proposed elements exhibit excellent performance even when applied to thin shell
problems; therefore, they offer a viable alternative to shell elements especially in problems
involving non-smooth shells.
2. None of the proposed elements locks at the nearly incompressible limit. This stands in sharp
contrast to other recently proposed brick elements. For example, the nine- and fifteen-parameter elements of Andelfinger et al.' and the nine-parameter element of Simo and Armero'
possess four eigenvalues that go to infinity as the incompressible limit is approached (five
when the elements are distorted). A locking-free element should exhibit only one eigenvalue
going to infinity as the nearly incompressible limit is approached.
The development of low-order high-accuracy finite elements has traditionally been one of the
main goals of finite element research, mainly because the computational effort associated with
element array processing increases linearly with the number of elements. The effort associated
with the solution of the global system of equations is proportional to the number of equations to
the third power. Therefore, a significant reduction in the computational effort can be achieved if
the number of global degrees-of-freedom is reduced, even at the expense of a somewhat larger
computational'effort at the element level.
CCC 0029-59811961142337-25
0 1996 by John Wiley & Sons, Ltd.
Received 28 April 1995
Revised 15 November 1995
2338
S. L. WEISSMAN
The design of low-order high-accuracy finite elements is based on mixed formulations. Two
approaches to mixed finite elements are commonly used: (1) a two-field (Hellinger-Reissner)
functional where the displacement and stress fields are the assumed independent variables (the
works of Pian and Tong3 and Punch and Atluri4 are representative of this approach); and (2)
a three-field (Hu-Washizu) functional, where the displacement, stress and strain fields are the
assumed independent variables. Representative of this approach are the works of Nagtegaal
et ~ l . Hughes,6
, ~
and Simo et aL7 The works of Malkus and Hughesa and Simo and Rifai’ have
shown that mixed formulations encompass, as particular cases, reduced and selective reduced
integration schemes (e.g. Reference lo), and the method of incompatible modes (e.g. Reference 11).
The design of mixed elements requires special care in the selection of the assumed fields so that
stability requirements are observed (e.g. References 12 and 13). Moreover, the use of mixed finite
elements does not guarantee high-accuracy elements, a fact that was cast into a limitation
principle by de Veubeke.I4 From a practical standpoint, it is easier to design the assumed stress
field than the assumed strain field to satisfy the above restrictions. Additionally, elements derived
from the two-field formulation are more computationally efficient. Consequently,elements based
on the Hellinger-Reissner functional are more popular. However, the Hellinger-Reissner functional is based on the complementary energy function, which for general constitutive laws is not
readily ascertainable. The Hu-Washizu functional (see e.g. Reference 15) is formulated in terms of
the strain energy function, thus making it more attractive for non-linear materials. For example,
integration schemes for elastoplasticity, such as the radial return algorithm for J2-plasticity of
Wilkins,16 are founded on the premise that the elastoplastic problem is strain driven. Consequently, in elements based on the Hu-Washizu variational principle, these integration algorithms
are formulated entirely in terms of the assumed total and plastic strain fields, and testing for
plastic loading can be performed independently at each material (integration) point. Furthermore, elements based on the Hellinger-Reissner functional are not suitable for the modelling of
perfect plasticity (see e.g. Reference 17). For these reasons, the more general Hu-Washizu
functional is used in this work.
An outline of the paper follows. Section 2 presents the strong (or local) form of the elastoplastic
boundary value problem (the Appendix contains a summary of the radial return mapping used in
the numerical simulations).Section 3 presents the weak form (includingthe construction procedure for the assumed stress and strain fields), derived from a three-field functional. Section
4 provides the assumed fields used to generate the elements proposed herein. In Section 5,
numerical results showing the excellent performance of the proposed elements are tabulated.
Finally, Section 6 presents concluding remarks.
2. THE BOUNDARY VALUE PROBLEM: STRONG FORM
This section presents the strong, or local, form of the boundary value problem addressed in this
work. The constitutive law is the rate-independent elastoplastic model, J2-plasticity with linear
isotropic and kinematic hardening, used in the examples provided in Section 5.
2.1. Kinematics
Let B c R3 denote a bounded open body whose closure has a (piece-wise)smooth boundary,
aB; where aB is assumed to have a (piece-wise)continuous with outward unit normal, n, and to
possess the following structure: a,BnaUB= 8 and a&%udvB
= aB where avBis the part of dB
on which displacements are prescribed as
UIaug= 0 (given)
(1)
HIGH-ACCURACY LOW-ORDER THREE-DIMENSIONAL BRICK ELEMENTS
2339
and 8 9 is the part of 8B on which traction is specified as
t = o*nla3 (given)
(2)
In the above equations, U(X) is the displacementfield associated with particles X c 1.
Assuming
small strains, the (total) strain field, E, is related to the displacement field, U,by
E:=
V'U
(3)
where V s is the symmetric part of the gradient operator. In component form, the total strain is
given by
Eij
= +(Vi,j
+
Vj,i)
(4)
Also, the usual assumption that
E:=
where E' and
E~
Ee
+ EP
(5)
are the elastic and plastic strain tensors, respectively, is adopted.
2.2. Balance of momentum
Let b : g --* R3and p : W + R be the (given) body force and mass density, respectively. The local
form of the balance of momentum (equilibrium)equations (neglectinginertia terms) are given by:*
in
(divo
+ pb = 0
dT = d
2.3. Constitution
A constitutive law that relates the stress tensor, a(X), to the strain tensor, E(X),is necessary to
complete the statement of the boundary value problem. To this end, a stored energy function
W (X,E(X)- E ~ ( X ) x) :Y~ + R is introduced so that the second rank stress tensor (CE 9)is
given by:'
o(X):= a, W(X,E(X) - EP(X))
(7)
and Y is the space of second rank symmetric tensors. The elastic tensor, C (a rank four tensor), is
defined by
c:=aEEw
(8)
and possesses the following symmetries: Cijkl= C k l i j = C i j l k = C j i l k . The strain energy function
W is assumed to be quadratic so that
W(Ec):=+EC:C:EC
(9)
and C is given in terms of the Lame coefficients, A and p, by
c = 118 1 + 2pI
The symmetry of the stress tensor is a result of the balance of angular momentum
If only initially isotropic and homogeneous materials are considered (as in this work), the explicit dependence of the
stress tensor on X can be dropped"
2340
S. L. WEISSMAN
where 1 and I are the rank two and rank four identity tensors, respectively. The Lame coefficients
are expressed in terms of Young’s modulus, E, and Poisson’s ratio, v, by
I=
(1
+
vE
v)(l
- 2v)
and p=-
E
2(1 v)
+
To complete the constitutive law, the evolution of the plastic strain and hardening parameters
,
x R“‘+ R be the yield function,: where rn is the dimension
must be described. Letf(a(8, E ~ ) q):Y
of the hardening parameters, vector q. An admissible state is such that
f (0, q) d 0
(12)
The surfacef(a, q) = 0 is known as the ‘yield surface’ in stress space. The evolution of the plastic
strain and hardening parameters are referred to as the ‘flow rule’ and ‘hardening law’, respectively. For the case of associative J,-plasticity, used in this work, the yield function, flow rule, and
hardening law are given by
where, following standard notations, the hardening vector, q, is taken as q:= { a , ij} where a is the
‘equivalent plastic strain’ that defines ‘isotropic hardening’ and q defines the centre or ‘kinematic
hardening’ of the von Mises yield surface (see, e.g. Reference 19); and q:= dev[a] - q, tr[q]:= 0,
and H’(a) and K ( a )are the kinematic and isotropic hardening moduli, respectively. In this work,
the following linear isotropic and kinematic hardening laws are adopted for H’(a) and K(a):
+
H’(a):= (1 - B)Ra and K(a):= ay BRa with BE[O, 13
(14)
Finally, i, is referred to as the ‘consistency parameter’, and is assumed to obey the following
Kuhn-Tacker complementary conditions:
I i, 2 O,f(a, q) d 0, and
i,f(a, q) = 0
I
In addition, f is constrained to satisfy the ‘consistency requirement’:
The Appendix provides an outline of the radial return mapping used to update the plastic strain
and hardening parameters along with the form of the algorithmic tangent. For further details
regarding these issues, see e.g. Reference 19.
3. THE WEAK FORM AND FINITE ELEMENT APPROXIMATION
This section presents a mixed finite element approximation to the elastoplastic boundary value
problem considered in Section 2. The formulation is based on a three-field (Hu-Washizu)
‘The explicit dependence of the stress on the total and plastic strains is included as a reminder that the total and plastic
strain tensors and hardening parameters are the independent variables in this development
HIGH-ACCURACY LOW-ORDER THREE-DIMENSIONAL BRICK ELEMENTS
2341
functional. For highly accurate elements in coarse meshes the approach advocated in Reference
20 and exploited in the context of two-dimensional elements in Reference 21 is adopted.
Accordingly, the assumed stress and strain fields are enlarged (over the minimum assumed fields
that satisfy stability conditions) and constrained so that the additional terms do not contribute
(weakly) to the internal energy, thus obtaining additional freedom to satisfy internal constraints.
3.1. Boundary value problem: weak form
Let the total free energy be given by the following three field functional:
I ~ ( o E,, U):=
I
[ W ( E- 9 )+ 0 : {V’U - E}] dY - nExT(U)
(17)
where the independent variables are the stress, strain and displacement fields. The plastic strain,
cP, is given by rate equations such as (13)2 (i.e. E~ is viewed as an internal variable; see e.g.
(U) is the potential energy of the external loading, which under the
Reference 22), and nExT
assumption of dead loading is given by
The weak form of the boundary value problem is obtained by looking at stationary points of the
functional (note that this is a constrained optimization problem), and is given by
Dn =
ss
+
r:[V”U - E ] ~ Y
( : [ a , W - o]dY
j9
M
kinematics
+
I
V’q:adV - D&xT.q = 0
M
constitution
(19)
V
balance of momentum
where r, ( and q are arbitrary admissible variations of the stress, strain and displacement fields,
respectively (note that q is required to satisfy the displacement boundary conditions). Since r, 5,
and q are arbitrary and independent of each other, it follows that equation (19) is composed of
three independent parts as indicated above. Furthermore, by the fundamental lemma of calculus,
the strong form of the strain-displacement,constitutive, and (after integration by parts) balance of
momentum equations are recovered.
3.2. Finite element approximation
In this work, a finite element approximation in which the assumed stress and strain fields are
approximated independently in each element, and thus are discontinuous across element boundaries, is presented.” A typical element’s domain .@:’ c a is obtained by a discretization of the
body, B, into finite dimensional regions such that a z up: .@:I, where ‘nel’ is the number of
elements in the mesh. Consequently, the total free energy is approximated by
nel
el
$By equation (19), only C-,continuity across the element boundary of the assumed stress and strain fields is strictly
required
2342
S. L. WEISSMAN
nr
where
(&, cC1, Ucl)is given by equation (17) with all fields and domain understood as element
variables and domain, respectively.The remainder of this work is carried out at the element level.
Consequently, the superscript 'el' denoting the element level can be omitted without any
confusion.
Finite element approximations can be constructed directly from equation (20) by introducing
explicit approximations to the assumed fields. Unfortunately, selecting the assumed fields to
achieve high accuracy in coarse meshes is not a trivial task in the presence of internal constraints.
To overcome this difficulty, the approach proposed by Weissman" is followed. Accordingly, the
assumed strain distribution is enlarged, over the minimum required for stability, and is constrained so that the added terms do not contribute to an energy-like quantity, consequently
providing the additional freedom required to satisfy internal constraints (e.g. incompressibility).
This constraint is appended to the free energy functional as follows:7
n(a,E, U):= J', [ W ( E-
EP)
+ cr:{VSU- + 3,:&W]dV - IIExT(U)
E}
(21)
Let the assumed stress, strain and displacement fields be approximated by
a:=
ss
E:=
E~
(22)
+ c2 = Elel + Ezez
(23)
U:= Nd
(24)
and
respectively, E~ is the added strain distribution. In the above equations, the assumed stress and
strain interpolations are constructed so that S and El contain the minimum distribution required
for stability, E l n E z = 8, and the e's and s's are the vectors of stress and strain unknowns,
respectively; the standard isoparametric interpolation is used for the displacement field where
N is the interpolation matrix and d is the vector of nodal displacements. Also let the Lagrange
multipliers be given by
A:= E'z
(25)
where E' is the interpolation matrix and z is a vector of independent parameters." The
added Euler-Lagrange equations due to 3, are used to relate Aez (the increment in ez)
to Ael by
Aezn= - (laEiTC:'gEzdV)-l[
+
jaEiT&W,,dV
EiTC:lgE1dVAel,
Recall that the functional in equation (21) represents the approximation to the true functional within the element's
domain, and all fields are finite dimensional. Consequently, the added constraint serves as an orthogonality restriction on
the acceptable distribution, rather than implying that d, W = 0. An additional constraint can be added explicitly in which
d, W is replaced by S (see Reference 20). This is not explicitly included in the present work because S and 1. are selected to
be orthogonal and thus, this added constraint is automatically satisfied
"Notethat dim(E') = dim(E,) and that the added constraint replaces a natural Euler-Lagrange equation due to the added
strain terms. This choice will be explained in Section 4 below
2343
HIGH-ACCURACY LOW-ORDER THREE-DIMENSIONAL BRICK ELEMENTS
respectively, where the subscript 'n' denotes the iteration number within the Newton iteration
used in solving the nonlinear problem, and C;lg is the algorithmic elastoplastic tangent (see the
Appendix). Let the following definitions be introduced:
STIElel,
r2,:=
+ E2ezn+ E2V, - Bd,] dV
and
r3,,:= J apbdV
+J
t dT - GTsln
(32)
where B is given by
B:= V"N.
(33)
Let S be chosen so that it is orthogonal to A,** the Lagrangian multiplier, and note that as a result
of equation (27) En too is orthogonal to A so that
The remaining (linearized) Euler-Lagrange equations can be cast into a matrix form as follows:
Eliminating the stress and strain coefficients yields the standard form:
GT(A,H;'AT)-'GAd,
L--y----J
= r3,
+ GT(A.H;lA;f)-l(H,,A~lrl.
+ rZn)
L
"
(36)
J
R"
K.
In the present work, A,, is constructed to be a fully ranked square matrix so that the stiffness
matrix, K,, and residual, R,, are given by
K,:= GTAiTH,A;'G
and R,:= r3,
+ GTA,T(rl,,+ H,,Al1r2,,)
(37)
respectively. For further details, including analysis of algorithmic stability, see Reference 20.
**This choice of S is motivated by numerical efficiency
2344
S. L. WEISSMAN
4. ASSUMED FIELDS SELECTION AND PROPOSED ELEMENTS
This section presents the assumed independent fields required for the generation of a family of
three-dimensional eight-node elements. These elements are constructed following the procedure
outlined in Section 3. Section 5 will show that these elements achieve the desired high accuracy in
coarse meshes while satisfying stability requirements.
4. I . Displacement interpolation
The displacement field is interpolated using the standard isoparametric shape functions (e.g.
Reference 23). Consequently, the displacement field is approximated by
8
where dI and NI are the vector of nodal displacement and shape function associated with node I,
respectively, and NIis given as a function of the element’s natural coordinates by
+ &O(l+ qrq)(l + <rO
NAt, q,
(39)
4.2. Stress interpolation
Following Pian and S~mihara,’~
the assumed stress field is first constructed in the element’s
natural coordinates (5, q, <) as
(40)
[9ii 9 2 2 y 3 3 9 1 2 9 2 3 9 3 i I T ’ %
where
r
1
0
0
Y=l
0
0
0
1
0
1
0
0
0
0
0 0 0
0 0
1 0 0
0 1 0
0 0 1
0 0 0
0
0
0
0
0
1
q
0
0
0
0
0
~ q < 0 0 0 0 0 0
0 0 0 ~ ~ ~ ~
0 0 0 0 0 ~ q ~
0 0 0 0 0 0 0 0 <
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0
0
0
q
0
0
0
11
1
0
0
0
0
~
q
0
0 0
0 0
0
(41)
0
is the vector of stress parameters associated with 9and thus its dimension is 18 (i.e. it contains
18 parameters).The assumed stress field is transformed from the (5, q, [) space into the (xl, x2, x3)
space by means of the following transformation:
aij =
YiIRfJXIJ
(42)
where the indices i a n d j take the values xl, x2 and x3; and the indices I and J take the values 5 q I;
and the ‘deformation gradient’ YiIis given by Yx,e
= axl/a5, etc. To preserve the ability to pass
the constant strain patch test, the transformation (42) is based on values .Fir at the element’s
center.24Consequently, the transformation (in matrix notation) is given by
S = Fox
(43)
0
2345
HIGH-ACCURACY LOW-ORDER THREE-DIMENSIONAL BRICK ELEMENTS
where
4.3. Strain interpolation
The stress interpolation in the element's natural coordinates, equation (41), is also used as the
strain g1interpolation (i.e. B1= Y ) ,while g2is given by
< < q e g o o 0 0 0
0 0 0 q q c g q o o
12
=
0 0
0 0
0 0
0 0
i
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
V
9 parameter
0
0
rcrq1:
0 0 0
0 0 0
0 0 0
J
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
-
g q o o o o
o o q ~ o
o
o o o o < q
15 parameter
The assumed strain interpolation in the (xl, x2,x3) space is obtained by using the inverse of
equation (43), and is given (in matrix form) by
E = F i T 1= FiT[&l 8 2 1
where C is the assumed strain interpolation in the element's natural coordinates.
(46)
2346
S. L. WEISSMAN
4.4. Lagrange multipliers interpolation
Analysis of equation (26) shows that stability is obtained for
E’= E2
(47)
(see Reference 20). This choice is also in agreement with the requirement that the enhancement to
the assumed fields does not contribute to the internal energy. Moreover, it naturally gives rise to
an Euler-Lagrange equation identical to the one due to the constraint introduced in equation
(21). However, this choice fails to satisfy equation (34), for elements with a non-constant
transformation Jacobian J = det(9). For this reason, the following modified interpolation is
used for the Lagrange multipliers:tt
4.5. Proposed elements
The assumed fields presented above are used to generate three-dimensional eight-node brick
elements. Three elements are presented that differ by the enhancement’s dimension. The three
options are 9, 15 and 24 parameters, as shown in equation (45), and accordingly are denoted by
B8-9P (8-node brick element with 9 parameters), B8-15P and B8-24P, respectively. The motivation for these selections is as follows: To overcome locking behaviour only the normal stresses/
strains must be enhanced; thus the selection of nine-parameters’ enhancement, which provides the
same approximation in the normal strains in all directions. In this case, however, the shear terms
are approximated by an incomplete linear distribution (in the element’s natural coordinates).
A better approximation is obtained if a complete linear distribution is adopted (hence the
15-parameter element). The third element is constructed to have the same distributions for both
the normal and shear terms, and thus the choice of the 24-parameter element. Finally, the stress
field can be used to formulate an element based on the Hellinger-Reissner formulation. This
element is identical to that proposed by Pian and Tong.3
4.6. Remarks on implementation
The procedure presented above is best suited for conveying the proposed ideas. From an
implementational standpoint, however, it is best to operate using the sparse matrices in the
element’s natural co-ordinates. Thus, for efficiency, the assumed field should not be pushed
forward to the physical space; rather, the tangent matrix and B-matrix should be pulled
back to the element’s natural co-ordinates (e.g B = FtB).
As will be shown in the next section, the last enhancement (i.e. extension to 24 parameters)
does not significantly improve the results over the B8-15P elements and therefore, should
not be used in practice. It is included precisely to demonstrate that adding modes does not
necessarily improve element performance.
++Notethat in the case of constant Jacobian elements the Lagrangian multiplier’s interpolation provided by equation (48)
is equivalent to that of equation (47)
HIGH-ACCURACY LOW-ORDER THREE-DIMENSIONAL BRICK ELEMENTS
2347
5. NUMERICAL EXAMPLES
The performance of the proposed elements is now evaluated with a number of challenging
problems selected from the literature. First, in Section 5.1, the elements are shown to pass the
constant strain patch test. Second, in Section 5.2, the eigenvalues of a single element's stiffness
matrix are examined to demonstrate that the proposed elements do not lock at the nearly
incompressible limit. Third, the higher-order patch test is used to evaluate the proposed elements'
performance in bending-dominated problems. The next three examples, in Sections 5.4, 5.5 and
5.6 are used to illustrate the excellent performance of the proposed elements in elastoplastic
problems (including the case of perfect plasticity). The following three examples, Sections 5.7,5.8
and 5.9, along with Section 5.6, are used to present the excellent performance of the proposed
three-dimensional brick elements in shell-like problems. Finally, Section 5.10 illustrates the
performance of the proposed elements in the context of finite strain elasticity.
5.1. Constant stress/strain patch test
A common objective in element design is that regardless of how skewed the element is, it should
model a state of constant strain exactly. This characteristic is demonstrated with the mesh
presented in Reference 25. The displacements at the vertices are specified to coincide with
a homogeneous state of strain (see Reference 25 for specified values), while the interior nodes'
displacements are left free to be determined. The material constants used are: Young's modulus,
E = 1.OE + 6, and Poisson's ratio, v = 0.25. All three proposed elements reproduce the exact
nodal displacements as specified in MacNeal and Harderz5 and thus, pass this test. The same
mesh is also used with only a minimum number of constrained displacements (to prevent rigid
body modes) and applied external loads that yield a constant state of stresslstrain. The proposed
elements pass this version of the patch test, too.
5.2. Eigenvalue analysis of a single element
To test the elements' performance at the nearly incompressible limit, and to show that they are
free of spurious zero energy modes, an eigenvalue analysis of a single element mesh is performed.
First, the eigenvalues of an undistorted element of side length L = 1 are presented in Table I.$$
Second, the eigenvalues of a rhombic-shaped element of side length L = 1 and distortion angle
@ = 45" are given in Table 11. A perspective view of the elements' geometry is shown in Figure 1.
The material properties used in both cases are: E = 1 and v = 0.4999.
It is evident from the tables that all elements possess the correct non-zero eigenvalues and that
all three are free of spurious locking behaviour at the nearly incompressible limit. The B8-9P
element, as expected, exhibits a higher sensitivity to mesh distortion, as manifested by the higher
eigenvalues (2-18). This result stands in sharp contrast to the nine-parameter elements of
Andelfinger et al.' and Simo and Armero.2 These latter elements possess four large eigenvalues,
three of which are undesired, indicating possible volumetric locking in some cases. Moreover, the
~
15-parameter element of Andelfinger et al.' and the 12-parameter element of Simo et ~ 1 . ' possess
four and three large eigenvalues, respectively. Locking behaviour of the Simo et a1.26element can
be avoided, however, through a modification of the isoparametric shape function derivatives (see
"The eigenvalues 19-24 are zero for both cases considered (i.e. the six rigid body modes)
2348
S. L. WEISSMAN
Table I. Eigenvalues for an undistorted element of side length L = 1
Eigenvalue
All elements
1
2500E
+3
2
3.333E - 1
3
3.333E - 1
4
3.333E - 1
5
3.333E - 1
6.
3.333E - 1
Eigenvalue
All elements
7
3.333E - 1
8
3.333E - 1
9
3.333E - 1
10
2.222E - 1
11
1.111E - 1
12
1.111E - 1
Eigenvalue
All elements
13
l.111E - 1
14
5.555E - 2
15
5.5558 - 2
16
5-555E - 2
17
5-555E - 2
18
5555E - 2
Table 11. Eigenvalues for a rhombic - shaped element of side length L = 1 and
4 = 45"
Eigenvalue
B8-24/15P
B8-9P
1
3.333E + 3
3.333E 3
2
6.372E - 1
7.690E - 1
3
5695E - 1
6.724E - 1
4
5WOE - 1
6.372E - 1
5
6
5'OOOE - 1 4303E - 1
5'000E - 1 5.000E - 1
Eigenvalue
B8-24/15P
B8-9P
7
1.962E - 1
3.934E - 1
. 8
1.897E - 1
3.545E - 1
1-603E - 1
3'150E - 1
9
10
1.3278 - 1
2'385E - 1
11
8.334E - 2
1'962E - 1
12
8.012E - 2
1.766E - 1
Eigenvalue
B8-24/15P
B8-9P
13
5.9658 - 2
8.334E - 2
14
15
2W0E - 2
6.647E - 2
1'069E - 2
4.755E - 2
16
9.849E - 3
2.508 - 2
17
4.277E - 3
2.387E - 2
18
3.736E - 3
1.607E - 2
+
Figure 1. Perspective views of an undistorted and a rhombic-shaped element
the eigenvalue analysis in References 1 and 26). A less expected result is that the BS-24P and
B8-15P elements yield identical results even for the distorted case. Finally, it is noted that the
eigenvalues obtained for the Pian and Tong3 element are identical to those obtained for the
B8-24P and B8-15P elements for both cases.
5.3. Bending of a cantilever beam: Higher order patch test
This is a standard benchmark problem used to assess element sensitivity to mesh distortion. In
this problem, a cantilever beam of length L = 10, height H = 2, and width W = 1 is loaded by an
end unit moment. To highlight the proposed elements' performance in the presence of internal
2349
HIGH-ACCURACY LOW-ORDER THREE-DIMENSIONALBRICK ELEMENTS
10.0
I,
0.5
0.5
Figure 2. Higher-order patch test mesh, plan view
c
5 1.0
g
0
9 0.8
0
a
BN
._
0.6
‘;ii 0.4
E
b 0.2
I
0.0 I
0
I
1
.
I
I
I
1
2
3
4
5
A
Figure 3. Higher order patch test, v = 0.49999; vertical displacement of the bottom edge, normalized with beam solution
constraints, the lateral displacements are fixed (i.e. a state of plane strain is simulated) and the
material constants are E = 1 and v = 0.4999 (i.e. nearly incompressible material). The beam is
modelled by two elements, as shown in Figure 2. Figure 3 shows the lower edge vertical
displacement, normalized with the exact beam theory solution of 56.2565, as A is increased from
0 to 5 (see Figure 2). To further emphasize that the deterioration in the elements’performance is
due to geometrical distortion, the same problem is repeated for the case of v = 0. The lower edge
vertical displacement (normalized with the exact beam solution of 75.0)is presented in Figure 4.
It is evident from comparing Figures 3 and 4 that the deterioration in performance is indeed
due to the geometrical distortion. In fact, the performance of the B8-9P element is slightly better
for the nearly incompressiblecase. The B8-15PYB8-24P and the Pian and Tong3 elements yield
identical results, and exhibit excellent insensitivity to mesh distortion. This identity is, however,
restricted to linear elastic constitutive law. Finally, the B9-9P element reveals significant sensitivity to excessive mesh distortion and thus, as expected, it is shown to be inferior relative to the
other proposed elements.
5.4. Conjined plastic flow
This problem was suggested by Andelfinger et al.’ as a challenging problem for elements’
performance under confined plastic flow. In this problem, a ‘regular brick’ (termed in Andelfinger
2350
S. L. WEISSMAN
I
0.04
0
.
I
I
1
2
I
.
3
4
5 A
Figure 4. Higher order patch test, v = 0; vertical displacement of the bottom edge, normalized with beam solution
Figure 5. A regular mesh of 5 x 5 x 5 elementsis used to model (taking advantage of symmetry conditions)a regular brick
of side length L = 100 and height h = 50, loaded by a uniform pressure load of 4 = 250 acting on its center area of 20 x 20
et al. as ‘regular block‘) of side length L = 100 and height h = 50 is fixed at its bottom, and
a uniform pressure of q = 250 is applied to the center 20 x 20 area of its top surface. Using
symmetry conditions, only one quadrant is discretized by a uniform 5 x 5 x 5 mesh, as shown in
Figure 5. The material properties are: E = 210000-0, v = 0-3, c,,= 250-0, B = 1.0 (i.e. isotropic
hardening), and R = 1/210.
Figure 6 presents the load-deflection curves obtained for the proposed elements and for the
displacement model. While the three proposed elements yield slightly different curves, they all
have the same slope at the end (this difference diminishes under mesh refinement).This is in sharp
contrast to load-deflection curve for the displacement model, which exhibits a higher slope. The
failure of the displacement model to achieve the correct slope is in agreement with the findings
of Nagtegaal et al.’ for the double-edge-notched tensile specimen, where this behaviour
HIGH-ACCURACY LOW-ORDER THREE-DIMENSIONAL BRICK ELEMENTS
2351
L
0
CI
0
rz
I.............
DlSP
Figure 6. Load-deflection curves for simulations of the regular brick problem using various elements
I
48
1
I
Figure 7. Cook's membrane problem (plan view)
was observed for very fine meshes. It should also be noted that the nine-parameter element of
Andelfinger et al.' also overestimates the load carrying capacity of the brick. However, the
twenty-one- and thirty-mode elements of Andelfinger et al.' yield results similar to the B8-15P
and B8-24P elements.
5.5. Cook's membrane problem: elastoplastic
A tapered panel (of unit thickness) is clamped at one end and loaded by an in-plane bending
load at the other end, as shown (plan view) in Figure 7. This benchmark problem demonstrates
the proposed elements' superior coarse mesh accuracy in bending-dominated elastoplastic problems. The material properties are: E = 70, v = 3, o,, = 0.243, = 0.1,and ff = 0-15. A plane strain
simulation is imposed by fixing the out-of-plane displacements. Figure 8 presents the convergence
(of the centre node on the loaded edge) under uniform mesh refinement. For comparison, the
results obtained for the plane strain implementation of the mean dilation (B-bar) element of
2352
S. L. WEISSMAN
2
4
8
16
32 NELlSide
Figure 8. Cook‘s membrane problem+lastoplastic plane strain simulation.Convergencefor the proposed elements and
the B-bar (plane strain) element
Nagtegaal et aL5 are provided.@The superiority of the proposed elements is evident from the
figures. It is interesting to note that even the B8-9P element exhibits superior performance over
the mean dilation element.
5.6. Clamped square plate: elastoplastic
This example, typically modelled by plate and shell elements, is introduced to demonstrate the
excellent performance of the proposed elements even when the element’s geometrical aspect ratio
is large. The plate’s characteristics are: side length L = 10 and thickness t = 0 1 (a plan view is
shown in Figure 9); the material properties are E = 10920.0, v = 0.3 and by = 1000.0 (i.e. an
elastic-perfectly plastic material); and the plate is loaded by a uniformly distributed pressure.
First, the convergence under uniform mesh refinement (with only one element through the
thickness) for the elastic case is considered. The transverse displacement of the mid-surface at the
plate’s centre is presented in Figure 10 (normalized with the Kirchhoff solution of w = 12.6). The
results obtained for the T1 plate element” are given for comparison. It is evident that the
three-dimensional brick elements compare favourably with the plate element results. Second, the
load-deflection curves for the mesh of the 16 x 16 (in-plane)elements with two and four elements
through the thickness are presented in Figure 11(stress control simulation). Comparison is made
with the upper Kirchhoff limit load,22given by
Fuppcr
= (45-05/&
h2ay
(49)
where h is the plate thickness. The results for the mesh with only one element through the
thickness are not reported since it does not have enough integration points through the thickness.
As demonstrated by Figure 10, the fact that there are only two integration points through the
thickness does not affect the accuracy in the (linear)elastic case. However, during the elastoplastic
#Recently a number of elements have been proposed that outperform the mean dilation element. However, the mean
dilation element remains the standard benchmark in the literature. Therefore, the mean dilation element is used here to
provide a measuring stick for the performance of the proposed elements
HIGH-ACCURACY LOW-ORDER THREE-DIMENSIONAL BRICK ELEMENTS
2353
Figure 9. Clamped square plate of side length L = 10, thickness h = 0.1. Due to symmetry only one quadrant is
discretized (plan view)
0.950
4
256
64
16
1024
NEL
Figure 10. Uniformlyloaded square clamped plate; convergenceof the centre transverse.displacementwith uniform mesh
refinement
B8-15P14
----o---.
om+ 0
1
5
.
I
.
.
I
.
a
- - - . - .
1
I
I
1
-
I
-1
10 15 20 25 30 35 40 45 50 55 60 65
W
Figure 1 1. Load-displacement curves for the clamped square plate-dastoplastic
2354
S.L. WEISSMAN
(a)
(b)
Figure 12. Circular plate meshes-plan view
I3
.--B
Q
1.00-
0.75
-
U
-
1 0.50 1
Y
L
-
U
41
.-
0.25
//
I____-----
-
Q
E
c
p
1
o.m!
1
I
I
100
Number of elements in the mesh
10
lo00
Figure 13. Clamped thin circular plate--convergence of the transverse displacement at R
refinement
=0
with uniform mesh
phase, more integration points through the thickness are required. In particular, when only one
element is used through the thickness, the integration points are relatively far from the upper and
lower surfaces; thus, the model predicts elastic, and consequently stiffer, response even when the
lower and upper surfaces have plastified.
5.7. Clamped thin circular plate
In this example a thin clamped circular plate of radius R = 5 and thickness t = 0.1 is subject to
a uniformly distributed unit load. This classical benchmark problem is used to demonstrate that
plate and shell elements are free from shear locking when irregular shaped elements are used.
Taking advantage of symmetry conditions, the plate is meshed using 3, 12,48 and 768 elements
(only one element through the thickness). Two mesh topology options, illustrated in Figure 12,
are examined. The linear elastic material properties are E = 10920.0 and v = 03.
The convergence of the transverse displacement at the plate’s centre (R = 0) with uniform mesh
refinement is shown in Figure 13 (normalized with the Kirchhoff solution of w = 9.78348).The
HIGH-ACCURACY LOW-ORDER THREE-DIMENSIONAL BRICK ELEMENTS
2355
t’
I F = 1.0
Figure 14. Pinched cylinder with end diaphragms, perspective view
results obtained for the T1 plate element2’ are given for comparison. The plate element
outperforms the brick elements. This result is expected as plate elements are specifically designed
for this class of problems. However, the proposed brick elements produce very good results. With
only twelve elements, 92 and 78 per cent of the Kirchhoff solution is obtained when using mesh
a and mesh b, respectively. Note that all the proposed elements yield virtually identical results.
Furthermore, observe the significant improvement in course meshes obtained when using meshes
of type a over meshes of type b. Finally, it should be reiterated that for this linear elastic problem,
the element of Pian and Tong3 yields identical results to these reported for the proposed elements.
5.8. Pinched cylinder with end diaphragms
A cylinder of length L = 600, radius R = 300, and thickness t = 3, with two end rigid
diaphragms, is loaded by pinching vertical forces at the middle section. The material properties
are E = 3.0 E 6 and v = 0.3. The objectives of this problem are: (1) to test the proposed
elements’ performance in a curved geometry, and (2) to test the importance of the geometrical
approximation. To this end, using symmetry conditions, only one octant of the cylinder is
discretized by uniform meshes, as shown in Figure 14.
First, the convergence under uniform mesh refinement is examined (only one element is used
through the thickness). Figure 15 presents the displacement under the load (normalized with an
analytical solution of 1432488 E - 5). The results obtained for the displacement model eight-node
brick element and for the shell element of Simo et aL2* (the mixed approximation) are provided
. ~ outperforms
~
the proposed brick
for comparison. The shell element of Simo et ~ 1 clearly
elements. However, roughly the same number of elements is required to achieve a low level of
error (e.g. 1 per cent). It should also be noted that the rate of convergence with mesh refinement is
about the same for the shell element and the proposed brick elements. Thus, the proposed
elements provide a viable alternative to shell elements. This alternative is particularly attractive in
problems involving non-smooth geometry, such as the one described below.
Second, the geometry is fixed to that obtained by the 4 x 4 mesh, and the convergence under
uniform mesh refinement of the displacement under the load is presented in Figure 16
+
2356
S. L. WEISSMAN
16 NEWSIDE
8
4
Figure 15. Pinched cylinder with end diaphragms. Convergence of the displacement under the load with uniform mesh
refinement
- -- - - - B8-24/15/9P
-
'
0.0
BE-DISP
CONVERGED
I
I
I
I
4
8
16
32
NEWSIDE
Figure 16. Pinched cylinder with end diaphragms-geometry fixed to that obtained by the 4 x 4 mesh. Convergence of the
displacement under the load with uniform mesh refinement
(normalized with a converged solution of 8-95 E - 6). The converged displacement in this case is
only 49 per cent of the analytical solution obtained for the smooth cylinder. This result re-affirms
the known sensitivity of curved shell-like structures to geometrical perturbations. Moreover, it
demonstrates the need for the proposed elements since such geometry is outside the scope of
classical shell theory (see e.g. Reference 29).
5.9. Pinched hemispherical shell with 18" hole
The pinched hemispherical shell with an 18" hole, proposed in Reference 25, is used to evaluate
the performance of the proposed elements in doubly curved structures. The shell's dimensions are:
radius R = 10 and thickness t = 0.04.The material properties are E = 6.875 E + 7 and v = 0-3.
The shell is subjected to alternating inward and outward forces (F = 2.0) 90" apart. Using
symmetry boundary conditions, only one quadrant of the shell is modelled by uniformed meshes,
as shown in Figure 17.
HIGH-ACCURACY LOW-ORDER THREE-DIMENSIONAL BRICK ELEMENTS
2357
t=
Figure 17. Hemispherical shell with 18" hole, perspective view
-g
1.0-
-.-EB 0.8n 0.6
FN
t
-
0.4 -
zb 0.2 4
16
32
NEL/SIDE
Figure 18. Hemispherical shell with 18" hole. Convergence of the displacement under the load with uniform mesh
refinement
Figure 18 presents the convergence of the displacement under the loads (normalizedwith 0.093,
a value quoted in Reference 28 as obtained by asymptotic expansiony for the proposed elements.
The results reported in Reference 28 for the mixed shell approximation and those obtained for the
displacement model eight-node brick elements are also provided for comparison. Once again the
shell element exhibits superior performance. However, as in the former example, roughly the same
number of elements is required to achieve a low level of error.
5.10. Finite strain elasticity: Cook's problem
The extension of the formulation presented herein to finite strain problems will be addressed in
a subsequent paper. However, in the context of finite strain elasticity, formulated in the reference
w MacNeal and Harderz5propose the value of 0.094 for the displacement under the load
2358
S.L. WEISSMAN
configuration, this extension is straightforward. This example is intended to illustrate the
performance of the proposed elements in this context. The changes required to transform the
small strain formulation presented herein into a finite strain one is to add the geometrical
stiffness,and to interpret the stress as the second Piola-Kirchhoff stress and the strain as the
Lagrangian strain.
To illustrate the performance of the proposed elements, consider once again Cook's problem of
Example 5.5, now with a compressible neo-Hook model characterized by the following strain
energy function:
K
w =Z
( t ( -~ 1)2 - In(J)) + P (tr(C) - 3)
Figure 19. Finite deformation analysis of Cook's problem; deformed mesh and outline of the undeformed mesh (dashed
line)
4
2
4
8
16
32
NEL./Side
Figure 20. Finite deformation analysis of Cook's problem; Convergence of the finite element solution
HIGH-ACCURACY LOW-ORDER THREE-DIMENSIONAL BRICK ELEMENTS
2359
where C is the Cauchy-Green strain tensor, J’ is the determinant of C, and K and p are the bulk
and shear moduli, respectively. As in Example 5.5, a plane strain condition is imposed by fixing
the out-of-plane displacements. The material properties are: K = 400942.0 and p = 80.1938 (i.e.
a nearly incompressiblematerial). The simulation is carried out using force control with steps of
AF = 5.0 and F E [0.0,250*0]. The deformed shape for the 64-element mesh is shown in Figure 19,
and the convergence of the top-corner vertical displacement with mesh refinement is given in
Figure 20.
6. CONCLUDING REMARKS
A family of three-dimensional eight-node brick elements that exhibit high accuracy in coarse
meshes and are suitable for the modelling of elastoplastic (including elastic-perfectly plastic)
problems has been presented. The elements are formulated based on a three-field functional of the
Hu-Washizu type; consequently, the classical strain-driven algorithms of plasticity are directly
applicable. This is in contrast with other highly accurate elements that are based on the
Hellinger-Reissner functional and require the development of special algorithms (see, e.g.
Reference 17). The latter are not suitable for the modelling of elastic-perfectly plastic problems.
Also, the elements derived via the Hellinger-Reissner functional are based on the complementary
energy, which is not always readily ascertainable.
The elements presented pass the constant stress/strain patch test, and exhibit good distortion
insensitivity. In addition, these elements exhibit good performance even when the element’s
geometrical aspect ratio is very large, as shown by numerous problems commonly used to
evaluate the performance of shell elements. Moreover, unlike the recently proposed elements of
Andelfinger et al.’ and Simo and Armero’ the nine-parameter element proposed herein does not
lock at the nearly incompressible limit.
Finally, the B8-15P and B8-24P elements yield almost identical results. Therefore, since no
advantage is to be gained by adding the extra nine modes, it is recommended that the B8-15P be
chosen for practical computations. (The B8-24P element was presented here to demonstrate that
adding the extra terms has only a marginal effect.)
APPENDIX: RADIAL RETURN MAPPING AND THE ALGORITHMIC TANGENT
The radial return mapping used to update the plastic strain and hardening parameters is
summarized in this Appendix. Also provided is the algorithmic elastoplastic tangent. For
a detailed discussion of integration algorithms for elastoplasticity and the algorithmic tangent,
see Simo and Hughes.’’
A.I. Return mapping
Step 1: Given L-:, a,, qn and
evaluate the trial deviator stress and q
trial
-
devCaln+1- 2pdev[En+1 - &:I
and
q::;
Step 2. Evaluate the yield function
= dev[a:?;
- q,,]
2360
S. L. WEISSMAN
Step 3. Check for plastic flow
< 0)O then
if(f:i
set (0). + = (o):?:
go to step 5
else
go to step 4
end if
Step 4. Compute Ay.+
= y.+
- y.
and update plastic strain and hardening parameters
where
(58)
Step 5. Compute the stress and algorithmic elastoplastic tangent
Q,+~ = k trace(&,+
c;y1
where k
=
{
+ dev[o]fp\
(59)
- 2pn::\
k 1 @ 1 + 2 , ~ ( 1 - +1@1)
k l @ l+2~16~(1-31@1)-2p6,n:Y\@n::~
if:::f
iff;?\
<0
>O
(60)
+ $,u is the bulk modulus, and the constants b1 and d2 are given by
=1
REFERENCES
1. U. Andelfinger, E. Ramm and D. Roehl, ‘2D- and 3D-enhanced assumed strain elements and their application in
plasticity’, Proc. COMPLAST 111, Barcelona, Spain, 1992.
2. J. C. Simo and F. Armero, ‘Geometrically nonlinear enhanced strain mixed methods and the method of incompatible
modes’, Int. j. numer. methods eng., 33, 1413-1449 (1992).
3. T. H. H. Pian and P. Tong, ‘Relations between incompatible displacement model and hybrid stress model’, Int. j.
numer. methods eng., 22, 173-181 (1986).
4. E. F. Punch and S. N. Atluri, ‘Development and testing of stable, invariant, isoparametric curvilinear 2- and 3-D
hybrid-stress elements’, Comp. Methods Appl. Mech. Eng., 47, 331-356 (1984).
5. J. C. Nagtegaal, D. M. Parks and J. C. Rice, ‘On numerically accurate finite element solutions in the fully plastic
range’, Comp. Methods Appl. Mech. Eng., 4, 153-177 (1974).
6. T. J. R. Hughes, ‘Generalization of selective integration procedures to anisotropic and nonlinear media’, Int. j. numer.
methods eng., 15, 1413-1418 (1980).
7. J. C. Simo, R. L. Taylor and K. S. Pister, ‘Variational and projection methods for volume constraint in finite
deformation elastoplasticity’, Comp. Methods Appl. Mech. Eng., 51, 177-208 (1985).
8. D. S. Malkus and T. J. R. Hughes, ‘Mixed finite element methods-reduced and selective integration techniques:
A unification of concepts’, Comp. Methods Appl. Mech. Eng., 15, 63-81 (1978).
9. J. C. Simo and M. S. Rifai, ‘A class of mixed assumed strain methods and the method of incompatible modes’, Int. j .
numer. methods eng., 29, 1595-1638 (1990).
HIGH-ACCURACY LOW-ORDER THREE-DIMENSIONAL BRICK ELEMENTS
2361
10. N. Bicanic and E. Hinton, ‘Spurious modes in two-dimensional isoparametric elements’, Int. j. numer. methods eng., 14,
1545-1557 (1979).
11. E. L. Wilson, R. L. F. W. P. Doherty and J. Ghaboussi, ‘Incompatible displacement models, in S. J. Fenves,
N. Perrone, A. R. Robinson and W. C. Schnobrich (eds.), Numerical and Computer Models in Structural Mechanics,
Academic Press, New York, 1973.
12. I. Babushka, ‘Error bounds for finite element method‘, Numer. Math., 16, 322-333 (1971).
13. F. Brezzi, ‘On the existence, uniqueness and approximation of saddle-point problems arising from lagrange multipliers, Rev. Francaise d’dutomati. Informat. Recherche Operationnelle, &R2, 129-1 51 (1974).
14. B. F. de Veubeke, ‘Displacement and equilibrium models in the finite element method’, in 0.C. Zienkiewicz and G. S.
Holister (eds.), Stress Analysis, John Wiley, London, 1965.
15. K. Washizu, Variational Methods in Elasticity and Plasticity, 3rd edn., Pergamon, Elmsford, NY, 1982.
16. M. L. Wilkins, ‘Calculation of Elastic-Plastic Flow’, in B. Alder, et al. (eds.), Methods of Computational Physics 3,
Academic Press, New York, 1964.
17. J. C. Simo, J. G. Kennedy and R. L. Taylor, ‘Complementary mixed finite element formulations of elastoplasticity’,
Comp. Methods Appl. Mech., Eng., 74, 177-206 (1989).
18. J. E. Marsden and T. J. R. Hughes, Mathematical Foundations of Elasticity, Englewood Cliffs, Prentice-Hall,
N.J., 1983.
19. J. C. Simo and T. J. R. Hughes, Elastoplasticity and Viscoplasticity: Computational Aspects, Springer, Berlin, to appear.
20. S. L. Weissman, ‘A mixed formulation of nonlinear-elastic problems’, Comput. Struct., 44, 813-822 (1992).
21. S. L. Weissman and M. Jamjian, ‘Two-dimensional elastoplasticity: Approximation by mixed finite elements’, Int.
j . numer. methods eng., 36, 3703-3727 (1993).
22. J. Lubliner, Plasticity Theory, Macmillan, New York, 1990.
23. T. J. R. Hughes, The Finite Element Method, Englewood Cliffs, Prentice-Hall, N.J., 1983.
24. T. H. H.Pian and K. Sumihara, ‘Rational approach for assumed stress finite elements’, Int. j . numer. methods eng., 20,
1685-1695 (1984).
25. R. H. MacNeal and R. L. Harder, ‘A proposed standard set of problems to test finite element accuracy’, Finite
Elements Anal. Des., 1, 3-20 (1985).
26. J. C . Simo, F. Armero and R. L. Taylor, ‘Improved versions of assumed enhanced strain tri-linear elements for
3D-finite deformation problems’, Comp. Methods Appl. Mech. Eng., 110, 359-386 (1993).
27. T. J. R. Hughes and T. E. Tezduyar, ‘Finite elements based upon mindlin plate theory with particular reference to the
four node bilinear isoparametric element’, J . Appl. Mech., 46, 587-596 (1981).
28. J. C. Simo, D. D. Fox and M. S. Rifai, ‘On a stress resultant geometrically exact shell model. Part 11: The linear theory;
computational aspects’, Comp. Methods Appl. Mech. Eng., 73, 53-92 (1989).
29. P. M. Naghdi, ‘Theory of shells’, in C. Truesdell (ed.), Handbuck der Physik, Vol. Via/2, Mechanics of Solids 11,
Springer, Berlin, 1983.
Документ
Категория
Без категории
Просмотров
4
Размер файла
1 164 Кб
Теги
894
1/--страниц
Пожаловаться на содержимое документа