Elastoplasticity

In this study, ideally perfect plasticity for small strains is taken into account.

In the case of small and large strain plasticity, the displacement fielducan be divided into elastic and plastic components as

u=u^e+u^p, (2.26)

whereu^eandu^p are the elastic and plastic displacements, respectively. In case of deformation, the linear strain tensorE^lcan be written as

E^l= 1

When Equation (2.26) is inserted into Equation (2.27), the strain tensor can be divided into two symmetrical tensors as

E^l=E^le+E^lp≡ 1 where E^le andE^lp are the elastic and plastic parts of the linear strain tensor.

It is noteworthy, however, that in Publication II the Green strain tensor is used in the definition of elastic forces to account for geometric nonlinearity, which makes it possible to account for the geometric stiffening effect, for example.

2.3 Constitutive models 41 For more information about the effect of geometric stiffening in different multi-body formulations, see [35, 44]. The Green strain can be presented using the displacement field (assuming that the coordinations{e1,e₂,e₃}and{e1,e₂,e₃}

It can be shown that the decomposition of the Green strain into elastic and plastic parts is more complex than in the case of the linear strain tensor E^l. Similarly to the decomposition of the linear strain tensor, the displacement field Equation (2.26) can be inserted into the Green strain which can be written as follows:

E=E^e+E^p+E^ep, (2.30)

whereE^eandE^pare the elastic and plastic parts of the Green strain tensor, and tensorE^epincludes components of both elastic and plastic parts. These tensors can be written as

As can be seen from Equation (2.33), both elastic and plastic components are combined in the tensor E^ep. Therefore, E^ep shows that it is simpler to use the displacement field decomposition for the linear strain tensor than for the nonlinear strain tensor. As a result, it is not possible to determine the elastic and plastic components of the strain tensor. This may become problematic when using this approach for large strain plasticity. Therefore, the multiplicative decomposition of the deformation gradient for defining the elastic and plastic components [9, p. 232] or, alternatively, the hypoelastic material model based on Eulerian rate theory [80] can be used. In addition, the behavior of the materials, such as metal under the influence of large strains, does not belong to the theory of

42 2 Absolute nodal coordinate formulation hyperelasticity [9, p. 231]. In Publication II, it is assumed that the combination of elastic and plastic components can be neglected, i.e.E^ep = 0. Therefore, in this study, the additive decomposition of the Green strain is used as

E=E^e+E^p. (2.34)

Numerical examples in Publication II are chosen in such a manner that the elastic and plastic strains are relatively small. The stress S depends on the elastic component of strain only, and it can therefore be defined as

S =⁴D:E^e=⁴D: (E−E^p). (2.35) In case of perfect plasticity, the Huber-von Mises yield condition is expressed as

Fy =p

3J2 −σy ≤0, (2.36)

whereJ₂ = ¹₂dev(S) : dev(S)is the second deviatoric stress invariant andσyis the yield stress. The stress deviator tensor is defined as

dev(S) =S− 1

3tr(S)I. (2.37)

The plasticity theory based on the Huber-von Mises criterion is often referred to asJ2-plasticity. The plastic strain rateE˙^p is presented with the associative flow rule as

E˙^p =γ∂Fy

∂S =γdev(S), (2.38)

whereγis the magnitude of the plastic strain and ^∂Fy

∂S is the gradient that defines the direction of the flow. If the yield condition satisfies Fy = 0, the behavior is plastic and, consequently, the increment of plastic strain can be computed from Equation (2.38). On the other hand, if the behavior is elastic, then the yield condition givesFy <0and the increment of plastic strain is zero. These relations are unified in the so-called Kuhn-Tucker complementary conditions as follows:

γ≥0, Fy ≤0 and γFy = 0. (2.39)

2.3 Constitutive models 43 In case of perfect plasticity, the elastic forces can be determined by dividing the volume integral into two parts based on total and plastic strains as follows:

δWint =δW_int^tot+δW_int^p whereδW_int^tot can be integrated by using Gaussian quadratures andδW_int^p can be defined by using plastic cells (subregions, sub volumes). The plastic cells have equal size and the plastic strain is constant over the cell. Instead of the plastic cells, one can use the plastic integration points and interpolate between them.

In order to find out the accurate distribution of plastic strains, the number of plastic cells is higher than that of the needed integration points. The number of plastic cells is defined so that the plastic strains are converged. The plastic strain distribution can be defined with an iterative procedure where the plastic strain for every cell at stepi+ 1is determined by integrating Equation (2.38) over time step [ti,ti+1] as follows: of plastic strain. At the beginning of the first time step of the iteration, the plastic strains are initialized to zero, and therefore, the generalized coordinates are solved using the total strain distributionE. At the end of the first sub step, the yield condition is checked and the plastic strains are defined at every plastic cell for the next sub step. The yield condition is verified at the end of every sub step and the computation of plastic strains is adopted in the iteration step of Newton’s iteration. As the initial value for the plastic strain distribution at the following time steps, the converged plastic strain distribution of the previous step is used.

(See more details about the algorithm used and the plastic cells [77].) 2.3.4 Other models

It is shown for continuum beam elements, where the cross-sectional deformation and rotations are defined by transverse gradient vectors, that general continuum mechanics may yield an incorrect result [57]. According to the work by Rhim and Lee, these incorrect results can be avoided by neglecting the Poisson effect

44 2 Absolute nodal coordinate formulation for continuum beam elements. This locking phenomenon is also found to be problematic in case of fully-parameterized elements based on the absolute nodal coordinate formulation [71, 29]. To neglect the Poisson effect, one can use the diagonal material stiffness matrix, which leads to similar results as obtained with conventional beam theory. A similar constitutive relation among the diagonal material stiffness matrix D = diag(E, E, E, µ, µ, µ) is also applied for beam elements based on the absolute nodal coordinate formulation [37, 19, 18]. Al-though this type of modified constitutive relation may be helpful in avoiding Poisson and thickness lockings for certain beam elements based on the absolute nodal coordinate formulation, it is not generally recommended for continuum elements. In case of plate elements, such a simplified constitutive relation will clearly lead to an incorrect solution, which can be seen from the convergence of the results in Publication I. More discussion and numerical results related to this plate element can be found in [43]. In the case of the continuum plate elements, different modified material stiffness matrices are introduced to avoid thickness locking, see for example [38]. When using the special case with Poisson’s ratio ν=0 in three-dimensional elasticity, or when using a modified material stiffness matrix, thickness locking can be avoided. The study [43] also shows that it can be used to approximate three-dimensional elasticity for thin plates modeled with fully-parameterized plate elements.

C

HAPTER

3

Summary of the findings

The main motivation of this dissertation is to propose numerical procedures in order to overcome and highlight problems associated with the finite elements based on the absolute nodal coordinate formulation. The publications in this dis-sertation focus, in part, on descriptions of the shear and transverse deformations in the finite elements based on the absolute nodal coordinate formulation. For the plate element studied in Publication I, the shear deformation is accounted for by employing bilinear shape functions. This linearization of shear deformations is used to avoid shear locking. In Publication I, a simplified material model was used in order to improve the numerical performance of the plate element. A comparison of fully-parameterized plates, where the St. Venant-Kirchhoff mate-rial model is employed, is shown in [43]. In Publication II, improvements to the stress and strain description in the higher order beam element are introduced. In the study, the shear stress distribution is defined by using a quadratic distribution in the transverse direction and linearly along the longitudinal direction of the beam. Usually, shear stress distribution is defined directly from the deformation field. In Publication III, Reissner’s nonlinear rod theory is implemented into the framework of the absolute nodal coordinate formulation using strain and kinetic energies, similar to the large rotation vector formulation by Simo and Vu-Quoc.

The transverse deformation is taken into account using the strain energy and is somewhat similar to the elastic line approach. In Publication IV, the transverse shear deformation is included in the Bernoulli-Euler type beam element based on the absolute nodal coordinate formulation by using the strain-displacement relations of the Timoshenko beam theory. This approach leads to a low amount of degrees of freedom, but suffers from slow convergence under large deformations.

45

46 3 Summary of the findings This same problem can also be seen in [70], where the Timoshenko beam theory with nonlinear strain-displacement relations are employed into the Bernoulli-Euler type beam element based on the absolute nodal coordinate formulation.

In Publication V, high frequencies due to transverse deformation are studied, and a linear continuum beam element is improved so that the Poisson phenomenon is accounted for without Poisson or thickness type locking. It is also shown that by eliminating the transverse frequencies, only minor improvements to compu-tational efficiency can be reached.

In this section, the main advantages and disadvantages of the studied elements based on the absolute nodal coordinate formulation are described. The first part of this section concentrates on studies of the beam elements and the second part on plate elements. Additionally, numerical problems due to the locking phenomena are briefly discussed in the both sections.

3.1 Studies of beam elements

Publication II

In Publication II, improvements to the description of the stress and strain distri-bution in the higher order beam element are introduced. The higher order beam element based on the absolute nodal coordinate formulation is introduced in [29], and it avoids shear locking phenomena due to the usage of higher order terms in the polynomial expansion. However, the shear deformation is directly defined from the deformation field and leads to constant transverse shear deformation.

This is due to the linear interpolation along the cross-section of the beam. It is noteworthy that the linear displacement field in the transverse direction may lead to incorrect results, provided that the shear correction factor is not employed. An accurate description for shear strains in the transverse direction is particularly important in case of nonlinear material or fatigue analysis. In the element, shear stresses are interpolated linearly along the longitudinal direction of the beam because it improves the accuracy of the shear stress. In Publication II, the improved shear stress is based on beam theory for transverse and torsional shear stresses. In order to define the shear stress due to shear force at any point in the cross-section of the beam, the quadratic distribution can be used as follows:

S_xy^S = 6Qy

3.1 Studies of beam elements 47 where W and H are width and height of the beam element, and shear force Qy =µAγxy and shear deformation γxy = r,x ·r,y. The shear stress S_xz^S is derived analogously to S_xy^S . According to the beam theory, only three com-ponents of the symmetric stress tensor make significant contributions to the beam performance. Therefore, the transverse normal stresses Syy andSzz are neglected as they are significant only in case of large surface loadings of the beam. In addition, it is assumed that the shear stressSyzis neglected in the cross-section. This is due to the assumption that the cross-section is not able to distort.

It shall be noted that a trapezoidal deformation of the cross-section does not exist in the beam elements based on the absolute nodal coordinate formulation, please see possible modes for example in [29]. The kinematic assumptions are fundamentally the same as those used in the Timoshenko theory [76]. In the improved beam element, the transverse shear stresses due to the torsion are defined by using the torsion function as follows:

S_xy^T =−2µκx

∂Φ

∂z; S_xz^T = 2µκx

∂Φ

∂y , (3.2)

whereΦ = Φ(x, y)is the torsion function andκxis the twist approximated using the definition in [59]. Due to the use of the torsion function, the axial stress due to torsion is neglected. A good approximation for the torsion function was found by using a moderate number of terms. For example, when using four terms for the polynomial expansion of torsion, the error of shear componentsS_xy^T andS_xz^T is approximately 1% at the integration points. In the improved beam element, the axial normal stress due to the axial and bending deformation is defined by using the theory of St. Venant and Kirchhoff, where the transverse normal strains are neglected and Poisson locking is therefore avoided.

Several numerical examples were used to demonstrate the accuracy of the pro-posed stress definitions. The reference solutions were obtained by using solid elements in the commercial software ABAQUS. In the case of the cantilever beam problem, it was found that maximum shear stresses (near the clamped end) differ by only about 10% from the reference solution. Accordingly, the proposed definitions for the stress and strain distribution improve the accuracy significantly when compared to the original beam element where shear stress is constant in the transverse direction. In elasto-plastic cases, the computational efficiency of higher order beam elements is highly dependent on the number of plastic cells, where comparative stress is defined as well as the order of integration for the elastic forces. In the case of the plane pendulum, each element has 32×32×8 plastic cells for height, width and length. Such discretization arrangement was

48 3 Summary of the findings needed to obtain accurate results for deformations and stresses. It was found that in this case, the computational efficiency was improved by 44.9%.

Publication III

In Publication II, the shear stress distribution was improved by using the beam theory for transverse and torsional shear stress, which is only valid for mod-erately thick beams and small strains. In Publication III, the geometrically exact beam element was implemented into the framework of the absolute nodal coordinate formulation. A geometrically exact beam theory such as Reissner’s nonlinear beam theory can be used for large strains. The implementation into this framework was done to demonstrate that the elements based on the absolute nodal coordinate formulation could be seen as a geometrically exact formulation, according to Reissner’s nonlinear beam theory. It shall be noted that in the formulation the transverse deformation can be restrained using constraints or by using additional terms for strain energy. However, in the numerical examples used for the verification of theories, the transverse deformation plays a minor role. In Publication III, two modifications of the original strain energy as intro-duced by Omar and Shabana were proposed.

In the first approach, Poisson locking was eliminated by using plane stress as-sumptions with an integration procedure. In the second approach, two different implementations based on strain energy used in the large rotation vector formula-tion by Simo and Vu-Quoc are presented. Therefore, the strains are described in the orthogonal moving frame{t1,t₂}attached to the cross-section of the beam.

In this approach, the position vector of an arbitrary point in the inertial frame {e₁,e₂}is defined as follows:

r =r₀+yt2, (3.3)

wherer₀is an arbitrary position on the elastic line (Figure 3.1) which is defined as follows:

r₀= (x+u1)e1+u2e₂, (3.4) where u1 and u2 are the components of displacement vector u. The moving vectors can be written in the inertial frame{e1,e₂}as follows:

ti =Rjiej, (3.5)

3.1 Studies of beam elements 49

Figure 3.1. Cross-section and elastic line in the beam element based on the absolute nodal coordinate formulation in the framework of the large rotation vector formulation. Points pand prefer to the same particle at different configurations after displacementu.

where values of theRjican be present in the matrix form as follows:

R=

cosθ −sinθ sinθ cosθ

, (3.6)

where θ is the time dependent angle of rotation that is defined by using the gradient vectorr,yin Publication III. It is shown in Publication III that the strains can be derived by using a pseudo-polar decomposition of the deformation gradi-ent [23, p. 110]. Unlike the case for classical polar decomposition, the rotation must be determined beforehand for pseudo-polar decomposition, and it can be defined as

F =R(I+H), (3.7)

where H is a non-symmetric strain matrix. The deformation gradient can be written by using the relationt_2,x =−t1θ,xas follows:

F = ∂r

∂x =

r_0,x−yθ,xt₁ t₂

. (3.8)

50 3 Summary of the findings The strain matrix can be written as

H=R^TF −I= t₁^T(r0,x−yθ,xt₁)−1 0 whereθ,xis the rate of rotation of the cross-section along the undeformed length of the beam, Γ1 is the strain distribution due to axial forces and the bending moment, and Γ2 describes the transverse shear deformation. The strains at the orthonormal moving base{t1,t₂}can be written as

Γ1 =t₁^Tr_,x−1−yθ,x= Γ⁰₁−yθ,x, (3.10)

Γ2=t₂^Tr_,x, (3.11)

whereΓ⁰₁ is the axial strain of the elastic line. It is noteworthy that the strain energies were described in the total Lagrangian scheme while no geometrical ap-proximations were used. In these implementations, an elastic line approach [59]

with a modified strain energy was used. Such an elastic line approach suffers from shear locking, and therefore, shear strainΓ2 was modified by using linear interpolations.

Based on Publication III, it can be concluded that beam elements based on the strain energy by Simo-Vu-Quoc produce the same convergence as the original beam element by Simo and Vu-Quoc, provided that the order of interpolations for the position and angle are identical in both approaches. In addition, compared to the classical nonlinear beam theory, it was determined that the numerical ineffi-ciency of the absolute nodal coordinate formulation was not due to the transverse eigenfrequencies, and accordingly, the description of transverse deformation will not lead to ill-conditioned finite elements.

Publication IV

In the fully-parameterized elements based on the absolute nodal coordinate for-mulation, shear deformation is usually captured through the introduction of ad-ditional gradient vectors in the element transverse direction. In Publication IV, a different approach to account for the transverse shear deformation is introduced.

The main benefit of the introduced approach is that shear deformation can be accounted for with a low number of degrees of freedom, leading to reduced computational effort. The definition of the transverse shear deformation is based

3.1 Studies of beam elements 51

Figure 3.2. Deformation of the infinitesimal section of the beam element based on the absolute nodal coordinate formulation in the framework of the Timoshenko theory.

on the Timoshenko theory [76]. In Publication IV, shear deformation is imple-mented into the two-dimensional Bernoulli-Euler type beam element based on the absolute nodal coordinate formulation using a well known approach to clas-sical beam elements [53, 50], where small displacement theory is used. In this formulation, there is no need for the use of transverse gradient vectors to capture

on the Timoshenko theory [76]. In Publication IV, shear deformation is imple-mented into the two-dimensional Bernoulli-Euler type beam element based on the absolute nodal coordinate formulation using a well known approach to clas-sical beam elements [53, 50], where small displacement theory is used. In this formulation, there is no need for the use of transverse gradient vectors to capture

In document Development of beam and plate finite elements based on the absolute nodal coordinate formulation (sivua 40-0)