Beam Elements (TYPE3)

Radioss uses a shear beam theory or Timoshenko formulation for its beam elements.

This formulation assumes that the internal virtual work rate is associated with the axial, torsional and shear strains. The other assumptions are:

No cross-section deformation in its plane.
No cross-section warping out of its plane.

With these assumptions, transverse shear is taken into account.

This formulation can degenerate into a standard Euler-Bernoulli formulation (the cross section remains normal to the beam axis). This choice is under user control.

Local Coordinate System

The properties describing a beam element are all defined in a local coordinate system.

This coordinate system can be seen in Figure 1. Nodes 1 and 2 of the element are used to define the local X axis, with the origin at node 1. The local Y axis is defined using node 3, which lies in the local XY plane, along with nodes 1 and 2. The Z axis is determined from the vector cross product of the positive X and Y axes.

The local Y direction is first defined at time $t = 0$ and its position is corrected at each cycle, taking into account the mean rotation of the X axis. The Z axis is always orthogonal to the X and Y axes.

Deformations are computed with respect to the local coordinate system displaced and rotated to take into account rigid body motion. Translational velocities

V

and angular velocities

Ω

with respect to the global reference frame are expressed in the local frame.

Beam Element Geometry

The beam geometry is user-defined by:

$A$: Cross section area
$I_{x}$: Area moment of inertia of cross section about local x axis
$I_{y}$: Area moment of inertia of cross section about local y axis
$I_{z}$: Area moment of inertia of cross section about local z axis

The area moments of inertia about the y and z axes are concerned with bending. They can be calculated using the relationships:(1)

I_{y} = \iint_{A} z^{2} d y d z

(2)

I_{z} = \iint_{A} y^{2} d y d z

The area moment of inertia about the x axis concerns torsion. This is simply the summation of the previous two moments of Ontario:(3)

I_{x} = I_{y} + I_{z}

Minimum Time Step

The minimum time step for a beam element is determined using the following relation:(4)

Δ t = \frac{a L}{c}

Where,

$c$ is the speed of sound: $\sqrt{E / ρ}$ ,

$a = \frac{1}{2} \min (\sqrt{\min (4, 1 + \frac{b}{12})} \cdot F_{1}, \sqrt{\frac{b}{3}} \cdot F_{2})$ ,

$F_{1} = \sqrt{1 + 2 d^{2}} - \sqrt{2} d$

$F_{2} = \min (F_{1}, \sqrt{1 + 2 d_{s}^{2}} - \sqrt{2} d_{s})$

$b = \frac{A L^{2}}{\max (I_{y}, I_{z})}$

$d = \max (d_{m}, d_{f})$

$d_{s} = d \cdot \max (1, \sqrt{\frac{12}{b}} \cdot \sqrt{1 + \frac{12 E}{\frac{5}{6} G b} (1 - I_{s h e a r})})$

Beam Element Behavior

Radioss beam elements behave in four individual ways:

Membrane or axial deformation
Torsion
Bending about the z axis
Bending about the y axis

Membrane Behavior

Membrane or axial behavior is the extension or compression of the beam element. The forces acting on an element are shown in Figure 2.

The force rate vector for an element is calculated using the relation:(5)

[\begin{matrix} {\dot{F}}_{x 1} \\ {\dot{F}}_{x 2} \end{matrix}] = \frac{E A}{l} [\begin{matrix} + 1 & - 1 \\ - 1 & + 1 \end{matrix}] [\begin{matrix} υ_{x 1} \\ υ_{x 2} \end{matrix}]

Where,

$E$: Elastic modulus
$l$: Beam element length
$υ_{x}$: Nodal velocity in x direction

With the force rate equation, the force vector is determined using explicit time integration:(6)

F_{x} (t + Δ t) = F_{x} (t) + {\dot{F}}_{x} Δ t

Torsion

Torsional deformation occurs when the beam is loaded with a moment M about the X axis as shown in Figure 3.

The moment rate vector is computed by:(7)

[\begin{matrix} {\dot{M}}_{x 1} \\ {\dot{M}}_{x 2} \end{matrix}] = \frac{G I_{x}}{l} [\begin{matrix} + 1 & - 1 \\ - 1 & + 1 \end{matrix}] [\begin{matrix} {\dot{θ}}_{x 1} \\ {\dot{θ}}_{x 2} \end{matrix}]

Where,

$G$: Modulus of rigidity
${\dot{θ}}_{x}$: Angular rotation rate

The moment about the X axis is found by explicit time integration:(8)

M_{x} (t + Δ t) = M_{x} (t) + {\dot{M}}_{x} Δ t

Bending About Z-axis

Bending about the z axis involves a force in the y direction and a moment about the z axis as shown in Figure 4.

Two vector fields must be solved for forces and moments. The rate equations are:(9)

[\begin{matrix} {\dot{F}}_{y 1} \\ {\dot{F}}_{y 2} \end{matrix}] = \frac{E I_{z}}{l^{3} (1 + ϕ_{y})} [\begin{matrix} + 12 & 6 l & - 12 & + 6 l \\ - 12 & - 6 l & + 12 & - 6 l \end{matrix}] [\begin{matrix} v_{y 1} \\ {\dot{θ}}_{z 1} \\ v_{y 2} \\ {\dot{θ}}_{z 2} \end{matrix}]

(10)

[\begin{matrix} {\dot{M}}_{z 1} \\ {\dot{M}}_{z 2} \end{matrix}] = \frac{E I_{z}}{l^{3} (1 + ϕ_{y})} [\begin{matrix} + 6 l & (4 + ϕ_{y}) l^{2} & - 6 l & (2 - ϕ_{y}) l^{2} \\ + 6 l & (2 - ϕ_{y}) l^{2} & - 6 l & (4 + ϕ_{y}) l^{2} \end{matrix}] [\begin{matrix} v_{y 1} \\ {\dot{θ}}_{z 1} \\ v_{y 2} \\ {\dot{θ}}_{z 2} \end{matrix}]

Where,

$ϕ_{y} = \frac{144 (1 + v) I_{z}}{5 A l^{2}}$ ,

$υ$ is the Poisson's ratio.

The factor $ϕ_{y}$ takes into account transverse shear.

The time integration for both is:(11)

F_{y} (t + Δ t) = F_{y} (t) + {\dot{F}}_{y} Δ t

(12)

M_{z} (t + Δ t) = M_{z} (t) + {\dot{M}}_{z} Δ t

Bending About Y-axis

Bending about the Y axis is identical to bending about the Z axis. A force in the Y direction and a moment about the Z axis, shown in Figure 5, contribute to the elemental bending.

The rate equations are:(13)

[\begin{matrix} {\dot{F}}_{z 1} \\ {\dot{F}}_{z 2} \end{matrix}] = \frac{E l_{y}}{l^{3} (1 + Φ_{z})} [\begin{matrix} + 12 \\ - 12 \end{matrix} \begin{matrix} 6 l \\ - 6 l \end{matrix} \begin{matrix} - 12 \\ + 12 \end{matrix} \begin{matrix} + 6 l \\ - 6 l \end{matrix}] [\begin{matrix} ν_{z 1} \\ {\dot{θ}}_{y 1} \\ ν_{z 2} \\ {\dot{θ}}_{y 2} \end{matrix}]

(14)

[\begin{matrix} {\dot{M}}_{y 1} \\ {\dot{M}}_{y 2} \end{matrix}] = \frac{E l_{y}}{l^{3} (1 + Φ_{z})} [\begin{matrix} + 6 l \\ + 6 l \end{matrix} \begin{matrix} (4 + Φ_{z}) l^{2} \\ (2 - Φ_{z}) l^{2} \end{matrix} \begin{matrix} - 6 l \\ - 6 l \end{matrix} \begin{matrix} (2 - Φ_{z}) l^{2} \\ (4 + Φ_{z}) l^{2} \end{matrix}] [\begin{matrix} ν_{z 1} \\ {\dot{θ}}_{y 1} \\ ν_{z 2} \\ {\dot{θ}}_{y 2} \end{matrix}]

Where, $Φ_{z} = \frac{144 (1 + ν) I_{y}}{5 A l^{2}}$ .

Like bending about the Z axis, the factor $Φ_{z}$ introduces transverse shear.

With the time integration, the expression is:(15)

F_{z} (t + Δ t) = F_{z} (t) + {\dot{F}}_{z} Δ t

(16)

M_{y} (t + Δ t) = M_{y} (t) + {\dot{M}}_{y} Δ t

Material Properties

A beam element may have two different types of material property:

Elastic
Elasto-plastic

Elastic Behavior

The elastic beam is defined using material LAW1 which is a simple linear material law.

The cross-section of a beam is defined by its area $A$ and three area moments of inertia $I_{x}$ , $I_{y}$ and $I_{z}$ .

An elastic beam can be defined with these four parameters. For accuracy and stability, the following limitations should be respected:(17)

L > \sqrt{A}

(18)

0.01 A^{2} < I_{y} < 100 A^{2}

(19)

0.01 A^{2} < I_{z} < 100 A^{2}

(20)

\begin{array}{l} 0.1 (I_{y} + I_{z}) < I_{x} < 10 (I_{y} + I_{z}) \end{array}

Elasto-plastic Behavior

A global plasticity model is used.

The main assumption is that the beam cross section is full and rectangular. Optimal relations between sections and section inertia are:(21)

12 I_{y} I_{z} = A^{4}

(22)

I_{x} = I_{y} + I_{z}

However, this model also gives good results for the circular or ellipsoidal cross-section. For tubular or H cross-sections, plasticity will be approximated.

Recommendations:(23)

L > \sqrt{A}

(24)

0.1 A^{4} < 12 I_{y} I_{z} < 10 A^{4}

(25)

0.01 < I_{y} / I_{z} < 100

(26)

0.5 (I_{y} + I_{z}) < I_{x} < 2 (I_{y} + I_{z})

Global Beam Plasticity

The elasto-plastic beam element is defined using material LAW2:(27)

σ_{y} = (A + B ε_{p}^{n}) (1 + C \ln \frac{\dot{ε}}{{\dot{ε}}_{0}})

The increment of plastic strain is:(28)

Δ ε_{p} = \frac{Δ W_{p l a s t i c}}{σ_{y}}

The equivalent strain rate is derived from the total energy rate:(29)

\dot{ε} = \frac{Δ W_{t o t a l}}{σ_{e q} Δ t}

Yield stress:(30)

σ_{e q} = \sqrt{\frac{F_{x}^{2}}{A^{2}} + \frac{3}{A} (\frac{M_{x}^{2}}{I_{x x}} + \frac{M_{y}^{2}}{I_{y y}} + \frac{M_{z}^{2}}{I_{z z}})}

σ_{e q} > σ_{y}

, you perform a radial return on the yield surface:(31)

F_{x}^{p a} = F_{x} \frac{σ_{y}}{σ_{e q}}

and for

i

= x, y, z:(32)

M_{i}^{p a} = M_{i} \frac{σ_{y}}{σ_{e q}}

Inertia Computation

The computational method of inertia for some kinds of elements as beam is particular as the inertia has to be transferred to the extremities of the beam. The nodal inertias are computed in function of the material density

ρ

, the cross-section area

S

, the element length

L

and the moments of inertia

I_{x x}, I_{y y}, I_{z z}

:(33)

M A X (\begin{matrix} (\frac{ρ S L}{2}) (\frac{L^{2}}{12}) + (\frac{ρ L}{2}) • M A X (\begin{matrix} I_{y y} & ; & I_{z z} \end{matrix}) & ; \begin{matrix} (\frac{ρ L}{2}) • I_{x x} \end{matrix} \end{matrix})