Frame elements with lumped plasticity

Frame elements are formulated in terms of the solution variables at user-defined end nodes and extra internal degrees of freedom associated with an internal node. The three-dimensional version of the element is discussed here. The two-dimensional version is found by appropriate reduction of the three-dimensional degrees of freedom. The element has three nodes (two user-defined and one internal), 12 external degrees of freedom, and three internal degrees of freedom. Each of the two end nodes has six external degrees of freedom: three displacements and three rotations. An internal node (at the center of the element) has three displacement degrees of freedom only, as shown in Figure 1.

Figure 1. Frame element in space.

The element is formulated in a local system with the x-direction representing the axial direction and the y- and z-directions representing the directions transverse to the frame element axis. In this local coordinate system the element's degrees of freedom can be written

The elastic response of the element is governed by Euler-Bernoulli beam theory. The displacement interpolation for the deflections transverse to the frame element axis (the y- and z-directions) uses fourth-order polynomials, allowing for quadratic variation of the curvature along the beam axis. Let $\xi \in \left[-1,1\right]$ be the isoparametric coordinate along the length of the beam. Then,

The transverse displacement interpolations incorporate exact solutions to force and moment loading at the ends and constant distributed loads along the beam axis (such as gravity loading). The displacement interpolation function along the frame element axis (the x-direction) is a second-order polynomial, allowing for linear variation of the axial strain along the frame element axis:

The twist rotation degree of freedom interpolation along the beam axis (rotation about the x-axis) is linear, allowing for constant twist strain:

The generalized strains, following Euler-Bernoulli beam theory, are

where ${ϵ}_x$ is the axial strain, ${\kappa }_z$ and ${\kappa }_y$ are the beam curvatures, and ${ϵ}_t$ is the twist strain. The 15 undetermined constants in the interpolation equations for the displacements are determined by introducing the nodal degrees of freedom; that is,

The interpolations in terms of the nodal degrees of freedom are described below in the section discussing the large-displacement formulation.

The strain-displacement relationship is written in matrix form as

where $\mathbf{B}$ is a 4 $\times$ 15 matrix and

The elastic stiffness matrix is integrated numerically and used to calculate 15 nodal forces and moments—12 forces/moments (also called generalized forces) associated with the two end nodes,

and three forces associated with the internal node,

The vector of forces and moments for the frame element can be written as

The elastic stiffness is, therefore, a $15\times 15$ matrix relating the force vector, $\mathbf{F}$, and the nodal displacement vector, $\mathbf{q}$:

Material properties of frame elements can, in general, be temperature dependent. Let us define the elastic strain vector as

where $\mathit{\epsilon }$ denotes the total strain and ${\mathit{\epsilon }}^{th}$ denotes the thermal expansion strain, where only the axial strain is nonzero and is given by

where

$\alpha \left(\theta \right)$

is the thermal expansion coefficient,

$\theta$

is the current temperature at the frame element section,

${\theta }^0$

is the reference temperature for $\alpha$, and

${\theta }_I$

is the user-defined initial temperature at this point (Initial Conditions).

The temperature field is defined by the user at the element's ends and is assumed to be linear along the element axis but constant within the element cross-section. If the thermal expansion coefficient is temperature dependent, it is evaluated at the nodes. Thermal strains are calculated at the element's end nodes, and thermal strains at the integration points are interpolated from the nodal points using appropriate interpolation schemes: axial strains are interpolated linearly, curvatures are interpolated quadratically, and twist strain is constant along the frame element axis.

Initial generalized strains, ${\mathit{\epsilon }}^{initial}$, are calculated from the initial generalized forces given by the user, using the relationship

where $\mathbf{D}$ denotes the $4\times 4$ material matrix evaluated at the nodal temperature $\theta$:

and A is the cross-section area; E is Young's modulus; G is the shear modulus; $I_{11}$, $I_{12}$, and $I_{22}$ are cross-section moments of inertia; and $J_{p}G$ is the torsional stiffness. The vector of generalized initial forces includes the following components:

Initial strains, when needed, are interpolated from the nodal values to the integration points using appropriate interpolators: linear for the axial component, quadratic for the bending components, and constant for the torsional component.

Abaqus integrates the elastic stiffness matrix numerically:

where temperature-dependent material properties are evaluated at the integration points, assuming a linear variation of temperature along the element axis.

For the simplest case of a temperature-independent material and a pipe cross-section, the elastic stiffness matrix can be integrated analytically to give:

Twisting moments at the end nodes:

Axial forces at the nodes:

Bending moments at the end nodes and transverse forces for all three nodes:

where ${\mathbf{K}}_{\mathbf{m}}^{\mathbf{e}}$ is the bending part of the elastic stiffness matrix and takes the following form:

We assume that the displacement and rotation increments admit an additive decomposition into elastic and plastic parts. Hence,

The total forces and moments result from the elastic constitutive relation

Introduce the lumped plasticity concept such that plastic deformation can develop at the beam external (end) nodes only and develops through plastic rotations (hinges) and plastic axial displacement at one or both nodes. Further assume that the plastic deformation at the external nodes can be caused by the interaction of all three moments and the axial force. Therefore, the vector of plastic deformation at those end nodes has the following form:

The yield function $\mathrm{\Phi }$, called here the plastic interaction surface, is written in terms of the nodal forces and moments. To calculate the increment of plastic deformation, the plastic interaction surface, $\mathrm{\Phi }$, is checked at each of the two ends of the frame element during the loading history. In general, the plastic interaction surface is a function of the sectional forces, its plastic cross-sectional capacities, and the hardening parameters. The frame element is elastic if the following conditions are fulfilled at both frame ends:

The frame element is elastic-plastic if the plastic interaction surface is exceeded at one or both frame ends:

Assuming associated plasticity with the direction of the increment of plastic deformation along the outward normal to the plastic interaction surface, the following relationship holds:

where $\mathrm{\Delta }{\lambda }_I$ denotes the magnitude of the plastic deformation at end I and ${\mathbf{V}}_I$ denotes the direction of the plastic flow at that end.

The hardening model follows a nonlinear kinematic hardening rule generalized from the linear Ziegler hardening law and the relaxation term (the recall term), which introduces the nonlinearity. For details on the hardening model, see Models for metals subjected to cyclic loading. Now introduce the generalized backstress vector, ${\mathit{\alpha }}_I$, which defines the origin of the moving plastic interaction surface, and define it as

Thus, the plastic interaction surface can be expressed as a function of the generalized force ${\mathbf{S}}_I$, ${\mathrm{\Phi }}_I={\mathrm{\Phi }}_I\left({\mathbf{S}}_I\right)$ at the end I, where for each component i,

The nonlinear kinematic hardening evolution law for the increment of the backstress takes the form

where

$C_i$ and ${\gamma }_i$ are sectional parameters for the ith plastic component, which must be calibrated from the test data defining the hardening response of the cross-section. The parameters $C_i$ are the initial kinematic hardening moduli, and the parameters ${\gamma }_i$ determine the rate at which the kinematic hardening modulus decreases with increasing plastic deformation. The test data required for this implementation are the values of the sectional force and moment components as a function of generalized plastic displacements. This will be either the axial force versus plastic axial displacement or a moment versus plastic rotation of a hinge. The data can be given as pairs of values—generalized sectional force versus conjugate generalized plastic displacement—or by specifying a yield stress for the section. The curve fitting algorithm will determine the parameters $C_i$ and ${\gamma }_i$ for each component i, since the hardening law is written separately for each sectional force component, $F_i$, where the range of i depends on the number of forces and moments entering the plastic interaction surface. Integrating the hardening rule, Equation 1, the following evolution law for the backstress is obtained:

where the backstress ${\alpha }_{iI}^0$ indicates the value of the backstress at the beginning of the increment.

Plastic interaction surfaces formulated in the generalized sectional variables depend on the cross-section profile. Frame elements with lumped plasticity are valid for tubular cross-sections only, and in the simplest form the interaction surface can be expressed as an ellipsoid in a space of four sectional components: axial force and three moments. Normalized with the ultimate forces and moments for each of the sectional components, the plastic interaction surface, ${\mathrm{\Phi }}_I$, can be written for each end I of the frame element as

where $N_{x0}$, $M_{x0}$, $M_{y0}$, and $M_{z0}$ represent the cross-sectional capacities at initial yield: the axial force and three moments, respectively. Any other cross-sectional profile for which plastic interaction can be approximated well enough by the above ellipsoidal surface can be used within the lumped plasticity concept for frame elements. For two-dimensional problems modeled with frame elements the plastic interaction condition becomes an ellipsoid in the axial force and bending moment plane. By checking the plastic interaction condition at any time of the deformation at both frame element ends, it is determined that if

1. ${\mathrm{\Phi }}_1<0$ and ${\mathrm{\Phi }}_2<0$, the frame element stays elastic.

2. ${\mathrm{\Phi }}_I<0$ and ${\mathrm{\Phi }}_J\ge 0(I\ne J)$, the frame element is elastic-plastic. If the plastic condition at the end J is exceeded, an iterative procedure is needed to find the final deformation state at the end of the increment. Either end I stays elastic and end J becomes plastic, or both ends become plastic.

3. ${\mathrm{\Phi }}_I\ge 0$ and ${\mathrm{\Phi }}_J\ge 0$, the frame element is elastic-plastic. If the plastic condition is exceeded at one or both nodes, an iterative procedure is needed to find the final deformation state at the end of the increment. Depending upon the ratio of plastic deformation at both ends, one or both ends will become plastic.

The integration of the plasticity model for frame elements follows the same general rule as described in Integration of plasticity models.

To solve for the value of the deformation and sectional forces at the end of the increment for an arbitrary load increment, an iterative process is required. To set up an appropriate Newton loop, the following relationships are used, with some of them linearized:

Elastic equilibrium equation:

where F stands for the generalized force at the end of the increment.

The associated flow rule:

The hardening evolution law:

The backstress definition:

The plastic interaction surfaces at both frame ends:

The frame elements admit large overall displacements and rotations; however, it is assumed that the strains are small. Accordingly, the nonlinear geometric formulation corresponds to Euler-Bernoulli beam theory superposed on a rotating reference frame. An Euler-Bernoulli displacement field, ${\mathbf{q}}^{eb}$, is defined relative to this rotation reference configuration, which causes straining. The Euler-Bernoulli displacement field is defined as follows.

Let $\overline{\mathbf{x}}$ and $\overline{\mathit{\varphi }}$ be the average position and average rotation of the frame element in the deformed configuration:

The motion of the rotating reference system is defined by the rigid body motion, where the translational part of the motion is the displacement vector $\left(\overline{\mathbf{x}}-{\mathbf{X}}_3\right)$, since ${\mathbf{X}}_3$ initially corresponds to the element centroid. The rotation part of the motion is the rotation matrix created from the average rotation vector:

The average rotation defines the rotation of the element's local directions from the reference values $\left(\mathbf{T},{\mathbf{N}}_1,{\mathbf{N}}_2\right)$ to the current values $\left(\mathbf{t},{\mathbf{n}}_1,{\mathbf{n}}_2\right)$ through

To define the strain-inducing rotation contributions, we multiplicatively decompose the rotation at each node,

where ${\mathbf{r}}_I$ is the strain-inducing part of the nodal rotation. Since the strains are assumed small, define the Euler-Bernoulli rotations at node I as the axial vector ${\widehat{\mathit{\varphi }}}_I^{eb}$, such that

Using Equation 8 in Equation 7, we can solve for ${\widehat{\mathit{\varphi }}}_I^{eb}$ by quaternion extraction. In components relative to the reference element coordinate directions

The Euler-Bernoulli displacement field is the difference between the position of the node relative to the element center and the reference position of the node relative to the reference center rotated to the current configuration by the average rotation:

or in components relative to the rotated local element coordinate directions

Once the equivalent Euler-Bernoulli displacements and rotations are determined from the nonlinear displacements and rotations, standard expressions following Euler-Bernoulli beam theory are used. The element interpolations are

The strain increments, following Euler-Bernoulli beam theory, are

where ${ϵ}_x$ is the axial strain, ${\kappa }_z$ and ${\kappa }_y$ are the bending strains, and ${ϵ}_t$ is the twist strain.

The user can supply hardening data for sectional forces in the form of pairs of values relating the axial nodal force to the plastic extension and the bending and twisting moments to the plastic rotations. If both the hardening data and the yield stress value for the frame section are omitted, Abaqus assumes that the frame element will remain elastic.

The curve fitting algorithm is used for the evolution hardening equation. At least three pairs of data are required for Abaqus to fit the curve and solve for the constants $C_i$ and ${\gamma }_i$ for each generalized force component, i, as shown in Figure 2.

The parameters ${\gamma }_i$ are scaled according to the following formula:

where $V_i^0$ denotes the component i of the plastic direction at the beginning of plastic deformation. For the plastic interaction surface in the form of Equation 2, the scaling factor is equal to $V_i^0=2/S_i^0$.

Figure 2. Hardening model for a frame element.