Beam element formulation

Beam element in undeformed configuration

In the undeformed configuration, we use capital letters for all quantities. We assume that the undeformed state has no warping, so the position of a material point is given by

\hat{X} (S, S^{α}) = X (S) + S^{α} N_{α} (S)

and the curvatures and twist are defined by

\begin{matrix} B_{α} & ​ = ϵ_{α}^{β} T \cdot \frac{d N_{β}}{d S}, \\ B & = N_{2} \cdot \frac{d N_{1}}{d S} = - N_{1} \cdot \frac{d N_{2}}{d S}, \end{matrix}

where $T$ is the unit vector orthogonal to $N_{1}$ and $N_{2}$ ; i.e.,

T = N_{1} \times N_{2} .

We assume that the section normal $T$ coincides with the beam tangent $S$ .

In the element the position of a point on the axis is interpolated from nodal positions $X^{N}$ with standard interpolation functions as

X = g_{N} (ξ) X^{N},

where $ξ$ is a parametric coordinate, typically varying between $- 1$ and 1 along the element. The beam axis tangent is readily computed as

S = \frac{d X}{d ξ} / \frac{d S}{d ξ} = \frac{d g_{N}}{d ξ} X^{N} / | \frac{d g_{M}}{d ξ} X^{M} | .

The section normals are interpolated from user-defined nodal normals $N^{N}$ . However, we cannot use a simple interpolation since this would not create integration point normals orthogonal to the beam tangent $S$ . Hence, we use a two-step approach. First, we create approximate normals by the interpolation

{\bar{N}}_{α} = g_{N} (ξ) N_{α}^{N} .

Then, we orthonormalize these vectors with respect to $S$ by

{\tilde{N}}_{α} = ({\bar{N}}_{α} - S S \cdot {\bar{N}}_{α}) / | ({\bar{N}}_{α} - S S \cdot {\bar{N}}_{α}) | (no summation)

and subsequently, with respect to each other by

N_{α} = ({\tilde{N}}_{α} + ϵ_{α}^{β} {\tilde{N}}_{β} \times S) / | ({\tilde{N}}_{α} + ϵ_{α}^{β} {\tilde{N}}_{β} \times S) | (no sum on α),

where it has been assumed that ${\bar{N}}_{α}$ and $S$ form a right-handed system. This provides $T = S$ .

The curvature and the twist in the initial configuration are calculated directly from ${\bar{N}}_{α}$ as

\begin{matrix} B_{α} & = ϵ_{β}^{α} T \cdot \frac{d {\bar{N}}_{β}}{d S} = ϵ_{α}^{β} T \cdot \frac{d g_{N}}{d ξ} N_{α}^{N} / \frac{d S}{d ξ}, \\ B & = \frac{1}{2} ϵ_{α}^{β} {\bar{N}}_{β} \cdot \frac{d {\bar{N}}_{α}}{d S} = \frac{1}{2} ϵ_{α}^{β} {\bar{N}}_{β} \cdot \frac{d g_{N}}{d ξ} N^{N} / \frac{d S}{d ξ} . \end{matrix}

The “average” twist is taken, since, in general,

{\bar{N}}_{2} \cdot \frac{d {\bar{N}}_{1}}{d S} \neq - {\bar{N}}_{1} \cdot \frac{d {\bar{N}}_{2}}{d S} .

The gradient of the normal vectors is then obtained as

\frac{d N_{α}}{d S} = ϵ_{α}^{β} (- B_{β} T + B N_{β});

and, therefore, the curvature and twist are also equal to

\begin{matrix} B_{α} & ​ = ϵ_{α}^{β} T \cdot \frac{d N_{β}}{d S}, \\ B & = N_{2} \cdot \frac{d N_{1}}{d S} = - N_{1} \cdot \frac{d N_{2}}{d S} . \end{matrix}

The procedure followed above to derive $N_{α}$ and $\frac{d N_{α}}{d S}$ is not unique but provides values that satisfy the proper orthonormality conditions. This procedure is followed only for the undeformed configuration. For subsequent configurations $n_{α}$ and $\frac{d n_{α}}{d S}$ are obtained individually by forward integration of the kinematic equations.

Change of position, warping, and normal direction

We assume that the position of the beam axis and the orientation of the normals can undergo (independent) changes. The change in position of the axis is described by the velocity vector $v = v (ξ)$ , which can be obtained from nodal velocities $v^{N}$ with the standard interpolation functions as

v (ξ) = g_{N} (ξ) v^{N} .

The change in orientation of the normals is described by the spin vector $ω = ω (ξ)$ , which is obtained from nodal spin vectors $ω^{N}$ with the same interpolation functions $g_{N} (ξ)$ , as

{\dot{n}}_{α} = ω (ξ) \times n_{α}, ω (ξ) = g_{N} (ξ) ω^{N} .

Rigid body motion is included since the original position $X (ξ)$ is obtained by the same interpolation as $v (ξ)$ . The rate of change of warping is also defined in terms of the rate of change of nodal warping $w^{N}$ with the standard interpolation functions as

\dot{w} (ξ) = g_{N} (ξ) {\dot{w}}^{N} .

The velocity and spin describe the rate of change of the position and orientation. Finite changes in position are obtained by integration of the velocities over a finite time increment as

Δ x (ξ) = \int_{t}^{t + Δ t} v (ξ) d t = g_{N} (ξ) \int_{t}^{t + Δ t} v^{N} d t = g_{N} (ξ) Δ X^{N} .

Similarly, for the warping,

Δ w (ξ) = \int_{t}^{t + Δ t} \dot{w} (ξ) d t = g_{N} (ξ) \int_{t}^{t + Δ t} {\dot{w}}^{N} d t = g_{N} (ξ) Δ w^{N} .

The spin is related to the rate of change of the rotation quaternion $q$ by

2 q^{†} \dot{q} = ω, q = (\cos \frac{θ}{2}, \sin \frac{θ}{2} n), q^{†} = (\cos \frac{θ}{2}, - \sin \frac{θ}{2} n),

where $θ = θ n$ is the total or Euler rotation.

The relation between spin and quaternion can be integrated exactly if it is assumed that the spin is constant over the time increment $Δ t$ . We then define

Δ ϕ = ω Δ t, Δ ϕ = | Δ ϕ |, n = Δ ϕ / Δ ϕ,

and the incremental rotation quaternion follows from

Δ q = (\cos \frac{Δ ϕ}{2}, \sin \frac{Δ ϕ}{2} n) .

This allows us to update the normal directions at the nodes and at stress points with

n_{α} = Δ q N_{α} Δ q^{†} .

In this equation we use the notation that $N_{α}$ is $n_{α}$ at the beginning of the current increment; i.e.,

N_{α} \overset{def}{=} n_{α} (t) .

For the entire motion the new normal directions can be formally expressed as

n_{α} = q N_{α} q^{†},

where $q$ is defined by the product rule

q = \prod_{i = 1}^{n} Δ q^{(i)} .

Here i is an increment and $Δ q^{(i)}$ is the rotation quaternion for that increment. The section normal $t$ is updated the same way.

Change of curvature and twist

The curvature and the twist involve the derivative of the normal vector with respect to S. From the update rule for $n_{α}$ ,

\frac{d n_{α}}{d S} = \frac{d Δ q}{d S} N_{α} Δ q^{†} + Δ q N_{α} \frac{d Δ q^{†}}{d S} + Δ q \frac{d N_{α}}{d S} Δ q^{†} .

The second term can be written as

Δ q N_{α} \frac{d Δ q^{†}}{d S} = {[\frac{d Δ q}{d S} N_{α}^{†} Δ q^{†}]}^{†} = - {[\frac{d Δ q}{d S} N_{α} Δ q^{†}]}^{†} .

Hence, the scalar parts of the first two terms cancel each other, and the vector parts are the same. Therefore,

\frac{d n_{α}}{d S} = 2 V (\frac{d Δ q}{d S} N_{α} Δ q^{†}) + Δ q \frac{d N_{α}}{d S} Δ q^{†} .

Because $Δ q$ is a rotation quaternion, its inverse is equal to its conjugate ( $Δ q^{- 1} = Δ q^{†}$ ), and, hence, we can write

\begin{matrix} \frac{d n_{α}}{d S} & ​ = 2 V (\frac{d Δ q}{d S} Δ q^{†} Δ q N_{α} Δ q^{†}) + Δ q \frac{d N_{α}}{d S} Δ q^{†} \\ ​ = V (2 \frac{d Δ q}{d S} Δ q^{†} n_{α}) + Δ q \frac{d N_{α}}{d S} Δ q^{†}, \end{matrix}

where $p = V (p)$ denotes the vector part of a quaternion. For the first term we use the relation

\frac{d Δ q}{d S} = (- \frac{1}{2} \sin \frac{Δ ϕ}{2} \frac{d Δ ϕ}{d S}, \frac{1}{2} \cos \frac{Δ ϕ}{2} n \frac{d Δ ϕ}{d S} + \sin \frac{Δ ϕ}{2} \frac{d n}{d S}) .

This leads to

2 \frac{d Δ q}{d S} Δ q^{†} = n \frac{d Δ ϕ}{d S} + \sin Δ ϕ \frac{d n}{d S} + (1 - \cos Δ ϕ) n \times \frac{d n}{d S} = Δ ϱ,

which is a vector. From the definitions of $Δ ϕ$ and $n$ it follows that

\begin{matrix} \frac{d Δ ϕ}{d S} & = Δ ϕ^{- 1} Δ ϕ \cdot \frac{d Δ ϕ}{d S} = n \cdot \frac{d Δ ϕ}{d S}, \\ \frac{d n}{d S} & = Δ ϕ^{- 1} \frac{d Δ ϕ}{d S} - Δ ϕ^{- 2} Δ ϕ n \cdot \frac{d Δ ϕ}{d S} = Δ ϕ^{- 1} [\frac{d Δ ϕ}{d S} - n n \cdot \frac{d Δ ϕ}{d S}] . \end{matrix}

These results provide

\begin{matrix} Δ ϱ & ​ = n n \cdot \frac{d Δ ϕ}{d S} + \frac{\sin Δ ϕ}{Δ ϕ} \frac{d Δ ϕ}{d S} - \frac{\sin Δ ϕ}{Δ ϕ} n n \cdot \frac{d Δ ϕ}{d S} + \frac{(1 - \cos Δ ϕ)}{Δ ϕ} n \times \frac{d Δ ϕ}{d S} \\ ​ = \frac{\sin Δ ϕ}{Δ ϕ} \frac{d Δ ϕ}{d S} + \frac{(1 - \cos Δ ϕ)}{Δ ϕ} n \times \frac{d Δ ϕ}{d S} + (1 - \frac{\sin Δ ϕ}{Δ ϕ}) n n \cdot \frac{d Δ ϕ}{d S} . \end{matrix}

We see that $Δ ϱ \Rightarrow \frac{d Δ ϕ}{d S}$ if $Δ ϕ \to 0$ .

For the second term we express $\frac{d N_{α}}{d S}$ in terms of the curvature and twist at the beginning of the increment, which yields

Δ q \frac{d N_{α}}{d S} Δ q^{†} = Δ q ϵ_{α}^{β} (- B_{β} T + B N_{β}) Δ q^{†} = ϵ_{α}^{β} (- B_{β} t + B n_{β}) .

Combining these terms,

\frac{d n_{α}}{d S} = Δ ϱ \times n_{α} + ϵ_{α}^{β} (- B_{β} t + B n_{β}) .

The current curvature and twist are, therefore,

b_{α} = ϵ_{α}^{β} t \cdot \frac{d n_{β}}{d S} = ϵ_{α}^{β} Δ ϱ \cdot (n_{β} \times t) - ϵ_{α}^{β} ϵ_{β}^{γ} B_{γ} = Δ ϱ \cdot n_{α} + B_{α}

and

b = n_{2} \cdot \frac{d n_{1}}{d S} = Δ ϱ \cdot (n_{1} \times n_{2}) + ϵ_{1}^{2} B = Δ ϱ \cdot t + B .

Hence, the current curvature and twist are updated by a summation over all increments as given by

\begin{matrix} b_{α}^{(n)} & ​ = \overset{n}{\sum_{i = 1}} Δ ϱ^{(i)} \cdot n_{α}^{(i)} + B_{α}, \\ b^{(n)} & ​ = \overset{n}{\sum_{i = 1}} Δ ϱ^{(i)} \cdot t_{α}^{(i)} + B . \end{matrix}

First variations

The first variations of the geometric quantities are readily obtained. Recall that

\begin{matrix} λ & ​ = | \frac{d x}{d S} |, \\ γ_{α} & = s \cdot n_{α} (with s = λ^{- 1} \frac{d x}{d S}), \\ b_{α} & ​ = ϵ_{α}^{β} t \cdot \frac{d n_{β}}{d S}, \\ b & ​ = n_{2} \cdot \frac{d n_{1}}{d S}, \\ χ & ​ = \frac{d w}{d S} . \end{matrix}

It follows that

\begin{matrix} δ λ & ​ = s \cdot \frac{d δ x}{d S}, \\ δ γ_{α} & = δ s \cdot n_{α} + s \cdot δ n_{α} = - λ^{- 1} γ_{α} δ λ + λ^{- 1} \frac{d δ x}{d S} \cdot n_{α} + s \cdot (δ ϕ \times n_{α}) \\ ​ \approx λ^{- 1} n_{α} \cdot \frac{d δ x}{d S} - ϵ_{α}^{β} n_{β} \cdot δ ϕ, \\ δ χ & ​ = \frac{d δ w}{d S}, \end{matrix}

where in the expression for $δ γ_{α}$ we have assumed that $γ_{α} ≃ 0$ and $s ≃ t$ . For the variations in the curvature and twist, we note that $δ t = δ ϕ \times t$ . Hence, it follows that

\begin{matrix} δ b_{α} & = ϵ_{α}^{β} (δ t \cdot \frac{d n_{β}}{d S} + t \cdot \frac{d δ n_{β}}{d S}) = ϵ_{α}^{β} [(δ ϕ \times t) \cdot \frac{d n_{β}}{d S} + t \cdot (δ ϕ \times \frac{d n_{β}}{d S} + \frac{d δ ϕ}{d S} \times n_{β})] \\ = ϵ_{α}^{β} \frac{d δ ϕ}{d S} \cdot (n_{β} \times t) = n_{α} \cdot \frac{d δ ϕ}{d S}, \\ δ b & = δ n_{2} \cdot \frac{d n_{1}}{d S} + n_{2} \cdot \frac{d δ n_{1}}{d S} = (δ ϕ \times n_{2}) \cdot \frac{d n_{1}}{d S} + n_{2} \cdot (δ ϕ \times \frac{d n_{1}}{d S} + \frac{d δ ϕ}{d S} \times n_{1}) \\ = \frac{d δ ϕ}{d S} \cdot (n_{1} \times n_{2}) = t \cdot \frac{d δ ϕ}{d S} . \end{matrix}

These terms will be used in the virtual work equation that will be discussed later. In using the above expressions the rotational quantities are obtained by interpolation from nodal variational quantities that are assumed to be valid for the velocity fields as

\begin{matrix} δ x (ξ) & ​ = g_{N} (ξ) δ x^{N}, \\ δ w (ξ) & ​ = g_{N} (ξ) δ w^{N}, \\ δ ϕ (ξ) & ​ = g_{N} (ξ) δ ϕ^{N} . \end{matrix}

Solution with Newton's algorithm

Newton's algorithm involves linearization of the incremental equations. The equations must be linearized around the current (latest) state. At the integration points these equations simply take the form

\begin{matrix} d λ & ​ = s \cdot \frac{d d x}{d S}, \\ d γ_{α} & ​ \approx λ^{- 1} n_{α} \cdot \frac{d d x}{d S} - ϵ_{α}^{β} n_{β} \cdot d ϕ, \\ d b_{α} & ​ = \frac{d d ϕ}{d S} \cdot n_{α}, \\ d b & ​ = \frac{d d ϕ}{d S} \cdot t, \\ d χ & ​ = \frac{d d w}{d S} . \end{matrix}

The corrections $d x$ and $d w$ are readily obtained from the nodal corrections with the interpolation functions $g^{N} (ξ)$ as

d x = g^{N} (ξ) d x^{N}, d w = g^{N} (ξ) d w^{N} .

The corrections $d ϕ$ are defined with respect to the end of the increment: we call such corrections “compound” corrections. However, it has been assumed that the total incremental rotation $Δ ϕ$ would be interpolated with $g^{N} (ξ)$ as

Δ ϕ (ξ) = g^{N} (ξ) Δ ϕ^{N} .

Linearization of this equation yields

d Δ ϕ (ξ) = g^{N} (ξ) d Δ ϕ^{N},

where $d Δ ϕ$ are corrections in $Δ ϕ$ in an additive sense such that

Δ ϕ^{n e w} = Δ ϕ^{o l d} + d Δ ϕ .

To relate the additive correction $d Δ ϕ$ to the compound correction $d ϕ$ , we use the formula obtained for $Δ ϕ$ to find

d ϕ = Δ q^{†} d q = \frac{\sin Δ ϕ}{Δ ϕ} d Δ ϕ + \frac{(1 - \cos Δ ϕ)}{Δ ϕ} n \times d Δ ϕ + (1 - \frac{\sin Δ ϕ}{Δ ϕ}) n n \cdot d Δ ϕ .

Similarly, at the nodes we find

\begin{matrix} d ϕ^{N} & ​ = \frac{\sin Δ ϕ^{N}}{Δ ϕ^{N}} d Δ ϕ^{N} + \frac{(1 - \cos Δ ϕ^{N})}{Δ ϕ^{N}} n^{N} \times d Δ ϕ^{N} \\ + (1 - \frac{\sin Δ ϕ^{N}}{Δ ϕ^{N}}) n^{N} n^{N} \cdot d Δ ϕ^{N} (no summation) . \end{matrix}

Now we assume that the difference in incremental rotation along a beam element is small; i.e.,

| Δ ϕ (ξ) - Δ ϕ^{N} | ≪ 1,

which implies that either

| n (ξ) - n^{N} | ≪ 1 and | Δ ϕ (ξ) - Δ ϕ^{N} | ≪ 1

| Δ ϕ (ξ) | ≪ 1 and | Δ ϕ^{N} | ≪ 1.

In the first case we can use the approximations $n (ξ) \approx n^{N} \approx \bar{n}$ and $Δ ϕ (ξ) \approx Δ ϕ^{N} \approx Δ \bar{ϕ}$ which, with the use of the interpolation functions for $d Δ ϕ$ , gives

\begin{matrix} d ϕ & ​ \approx g_{N} (ξ) [\frac{\sin Δ \bar{ϕ}}{Δ \bar{ϕ}} d Δ ϕ^{N} + \frac{(1 - \cos Δ \bar{ϕ})}{Δ \bar{ϕ}} \bar{n} \times d Δ ϕ^{N} + (1 - \frac{\sin Δ \bar{ϕ}}{Δ \bar{ϕ}}) \bar{n} \bar{n} \cdot d Δ ϕ^{N}], \\ d ϕ^{N} & ​ \approx \frac{\sin Δ \bar{ϕ}}{Δ \bar{ϕ}} d Δ ϕ^{N} + \frac{(1 - \cos Δ \bar{ϕ})}{Δ \bar{ϕ}} \bar{n} \times d Δ ϕ^{N} + (1 - \frac{\sin Δ \bar{ϕ}}{Δ \bar{ϕ}}) \bar{n} \bar{n} \cdot d Δ ϕ^{N} . \end{matrix}

In the second case it follows directly that

d ϕ \approx g_{N} (ξ) d Δ ϕ^{N}, d ϕ^{N} \approx d Δ ϕ^{N} .

In either case

d ϕ (ξ) \approx g_{N} (ξ) d ϕ^{N} .

This approximate relationship is used only in the creation of the Jacobian, so the approximation will at worst result in reduced speed of convergence.

Once a nodal update vector $d ϕ^{N}$ is obtained, an exact update procedure is followed. This is achieved by a transformation into quaternions, use of exact quaternion update formula, and transformation of the results back into an incremental Euler rotation vector:

d q^{N} = (d q^{N}, d q^{N}) = (\cos \frac{d ϕ^{N}}{2}, \frac{d ϕ^{N}}{d ϕ^{N}} \sin \frac{d ϕ^{N}}{2}), (no summation)

Δ q^{N} = (Δ q^{N}, Δ q^{N}) = (\cos \frac{Δ ϕ^{N}}{2}, \frac{Δ ϕ^{N}}{Δ ϕ^{N}} \sin \frac{Δ ϕ^{N}}{2}), (no summation)

\begin{matrix} Δ q_{n e w}^{N} & ​ = d q^{N} Δ q^{N} \\ ​ = (d q^{N} Δ q^{N} - d q^{N} \cdot Δ q^{N}, d q^{N} Δ q^{N} + d q^{N} Δ q^{N} + d q^{N} \times Δ q^{N}), (no summation) \end{matrix}

Δ ϕ_{n e w}^{N} = 2 \frac{q_{n e w}^{N}}{| q_{n e w}^{N} |} \tan^{- 1} (| q_{n e w}^{N} |, q_{n e w}^{N}) (no summation) .

The incremental rotation vector at the integration point is obtained by interpolation. Subsequently, we can calculate the updated integration point normals and the incremental curvature and twist.

Second variations

For the calculation of the Jacobian, we also need the second variation in the generalized quantities. These follow from the first variations:

\begin{matrix} d δ λ = & d s \cdot \frac{d δ x}{d S} = λ^{- 1} \frac{d d x}{d S} \cdot (I - ss) \cdot \frac{d δ x}{d S}, \\ d δ γ_{α} = ​ & λ^{- 2} γ_{α} d λ δ λ - λ^{- 1} γ_{α} d δ λ - λ^{- 1} δ λ d γ_{α} - λ^{- 2} d λ \frac{d δ x}{d S} \cdot n_{α} \\ + λ^{- 1} \frac{d δ x}{d S} \cdot (d ϕ \times n_{α}) - λ^{- 1} d λ s \cdot (δ ϕ \times n_{α}) + λ^{- 1} \frac{d d x}{d S} \cdot (δ ϕ \times n_{α}) \\ + \frac{1}{2} s \cdot [δ ϕ \times (d ϕ \times n_{α}) + d ϕ \times (δ ϕ \times n_{α})] \\ \approx & - λ^{- 1} (δ λ d γ_{α} + d λ δ γ_{α}) + \frac{1}{2} δ ϕ \cdot (n_{α} t + t n_{α}) \cdot d ϕ \\ - λ^{- 1} ϵ_{α}^{β} (δ λ d ϕ \cdot n_{β} + d λ δ ϕ \cdot n_{β}) + λ^{- 1} ϵ_{α}^{β} (\frac{d δ x}{d S} \cdot n_{β} t \cdot d ϕ + \frac{d d x}{d S} \cdot n_{β} t \cdot δ ϕ) \\ \approx & - λ^{- 2} \frac{d δ x}{d S} \cdot (n_{α} t + t n_{α}) \cdot \frac{d d x}{d S} + λ^{- 1} ϵ_{α}^{β} (\frac{d δ x}{d S} \cdot n_{β} t \cdot d ϕ + δ ϕ \cdot t n_{β} \cdot \frac{d d x}{d S}) \\ + \frac{1}{2} δ ϕ \cdot (n_{α} t + t n_{α}) \cdot \underset{\cdot}{ϕ}, \end{matrix}

where we have again used $γ_{α} \approx 0$ .

For the second variation of the curvature we find

\begin{matrix} d δ b_{α} = ​ & ϵ_{α}^{β} [d δ t \cdot \frac{d n_{β}}{d S} + δ t \cdot \frac{d d n_{β}}{d S} + d t \cdot \frac{d δ n_{β}}{d S} + t \cdot \frac{d d δ n_{β}}{d S}] \\ = & ϵ_{α}^{β} [\frac{1}{2} δ ϕ \cdot (t \frac{d n_{β}}{d S} + \frac{d n_{β}}{d S} t) \cdot d ϕ - (t \cdot \frac{d n_{β}}{d S}) (δ ϕ \cdot d ϕ) \\ + (δ ϕ \times t) \cdot (\frac{d d ϕ}{d S} \times n_{β} + d ϕ \times \frac{d n_{β}}{d S}) + (d ϕ \times t) \cdot (\frac{d δ ϕ}{d S} \times n_{β} + δ ϕ \times \frac{d n_{β}}{d S}) \\ + t \cdot (\frac{1}{2} \frac{d d ϕ}{d S} δ ϕ \cdot n_{β} + \frac{1}{2} d ϕ \frac{d δ ϕ}{d S} \cdot n_{β} + \frac{1}{2} d ϕ δ ϕ \cdot \frac{d n_{β}}{d S}) \\ + t \cdot (\frac{1}{2} \frac{d δ ϕ}{d S} d ϕ \cdot n_{β} + \frac{1}{2} δ ϕ \frac{d d ϕ}{d S} \cdot n_{β} + \frac{1}{2} δ ϕ d ϕ \cdot \frac{d n_{β}}{d S}) \\ - t \cdot (\frac{d δ ϕ}{d S} \cdot d ϕ n_{β} + δ ϕ \cdot \frac{d d ϕ}{d S} n_{β} + δ ϕ \cdot d ϕ \frac{d n_{β}}{d S})] . \end{matrix}

Rewriting the second and third terms and combining the others yields

\begin{matrix} d δ b_{α} = ​ & ϵ_{α}^{β} [δ ϕ \cdot (t \frac{d n_{β}}{d S} + \frac{d n_{β}}{d S} t) \cdot d ϕ - 2 (t \cdot \frac{d n_{β}}{d S}) (δ ϕ \cdot d ϕ) \\ + \frac{1}{2} \frac{d δ ϕ}{d S} \cdot (n_{β} t + t n_{β}) \cdot d ϕ + \frac{1}{2} δ ϕ \cdot (n_{β} t + t n_{β}) \cdot \frac{d d ϕ}{d S} \\ - δ ϕ \cdot n_{β} t \cdot \frac{d d ϕ}{d S} + (t \cdot \frac{d n_{β}}{d S}) (δ ϕ \cdot d ϕ) - δ ϕ \cdot \frac{d n_{β}}{d S} t \cdot d ϕ \\ - d ϕ \cdot n_{β} t \cdot \frac{d δ ϕ}{d S} + (t \cdot \frac{d n_{β}}{d S}) (d ϕ \cdot δ ϕ) - d ϕ \cdot \frac{d n_{β}}{d S} t \cdot δ ϕ] \\ = & \frac{1}{2} ϵ_{α}^{β} [\frac{d δ ϕ}{d S} \cdot (n_{β} t - t n_{β}) \cdot d ϕ - δ ϕ \cdot (n_{β} t - t n_{β}) \cdot \frac{d d ϕ}{d S}] . \end{matrix}

The second variation of twist is

\begin{matrix} d δ b = ​ & d δ n_{2} \cdot \frac{d n_{1}}{d S} + δ n_{2} \cdot \frac{d d n_{1}}{d S} + d n_{2} \cdot \frac{d δ n_{1}}{d S} + n_{2} \cdot \frac{d d δ n_{1}}{d S} \\ = & \frac{1}{2} δ ϕ \cdot (n_{2} \frac{d n_{1}}{d S} + \frac{d n_{1}}{d S} n_{2}) \cdot d ϕ - (n_{2} \cdot \frac{d n_{1}}{d S}) (δ ϕ \cdot d ϕ) \\ + (δ ϕ \times n_{2}) \cdot (\frac{d d ϕ}{d S} \times n_{1} + d ϕ \times \frac{d n_{1}}{d S}) + (d ϕ \times n_{2}) \cdot (\frac{d δ ϕ}{d S} \times n_{1} + δ ϕ \times \frac{d n_{1}}{d S}) \\ + n_{2} \cdot (\frac{1}{2} \frac{d d ϕ}{d S} δ ϕ \cdot n_{1} + \frac{1}{2} d ϕ \frac{d δ ϕ}{d S} \cdot n_{1} + \frac{1}{2} d ϕ δ ϕ \cdot \frac{d n_{1}}{d S}) \\ + n_{2} \cdot (\frac{1}{2} \frac{d d ϕ}{d S} d ϕ \cdot n_{1} + \frac{1}{2} δ ϕ \frac{d d ϕ}{d S} \cdot n_{1} + \frac{1}{2} δ ϕ d ϕ \cdot \frac{d n_{1}}{d S}) \\ - n_{2} \cdot (\frac{d δ ϕ}{d S} \cdot d ϕ n_{1} + δ ϕ \cdot \frac{d d ϕ}{d S} n_{1} + δ ϕ \cdot d ϕ \frac{d n_{1}}{d S}) . \end{matrix}

Then, following the same procedure as for $d δ b_{α}$ ,

\begin{matrix} d δ b = ​ & δ ϕ \cdot (n_{2} \frac{d n_{1}}{d S} + \frac{d n_{1}}{d S} n_{2}) \cdot δ ϕ - 2 (n_{2} \cdot \frac{d n_{1}}{d S}) (δ ϕ \cdot d ϕ) \\ + \frac{1}{2} \frac{d δ ϕ}{d S} \cdot (n_{1} n_{2} + n_{2} n_{1}) \cdot d ϕ + \frac{1}{2} δ ϕ \cdot (n_{1} n_{2} + n_{2} n_{1}) \cdot \frac{d d ϕ}{d S} \\ - δ ϕ \cdot n_{1} n_{2} \cdot \frac{d d ϕ}{d S} + (n_{2} \cdot \frac{d n_{1}}{d S}) (δ ϕ \cdot d ϕ) - δ ϕ \cdot \frac{d n_{1}}{d S} n_{2} \cdot d ϕ \\ - d ϕ \cdot n_{1} n_{2} \cdot \frac{d d ϕ}{d S} + (n_{2} \cdot \frac{d n_{1}}{d S}) (d ϕ \cdot δ ϕ) - d ϕ \cdot \frac{d n_{1}}{d S} n_{2} \cdot δ ϕ \\ ​ = \frac{1}{2} [\frac{d δ ϕ}{d S} \cdot (n_{1} n_{2} - n_{2} n_{1}) \cdot d ϕ - δ ϕ \cdot (n_{1} n_{2} - n_{2} n_{1}) \cdot \frac{d d ϕ}{d S}] . \end{matrix}

Finally, we note

d δ χ = 0.

The resulting second variations are summarized as

\begin{matrix} d δ λ & ​ = λ^{- 1} \frac{d d x}{d S} \cdot (I - ss) \cdot \frac{d δ x}{d S}, \\ d δ γ_{α} & ​ = - λ^{- 2} \frac{d δ x}{d S} \cdot (n_{α} s + s n_{α}) \cdot \frac{d d x}{d S} + λ^{- 1} ϵ_{α}^{β} (\frac{d δ x}{d S} \cdot n_{β} t \cdot d ϕ + δ ϕ \cdot n_{β} t \cdot \frac{d d x}{d S}) \\ + \frac{1}{2} δ ϕ \cdot (n_{α} t + t n_{α}) \cdot \underset{\cdot}{ϕ}, \\ d δ b_{α} & ​ = \frac{1}{2} ϵ_{α}^{β} [\frac{d δ ϕ}{d S} \cdot (n_{β} t - t n_{β}) \cdot d ϕ - δ ϕ \cdot (n_{β} t - t n_{β}) \cdot \frac{d d ϕ}{d S}], \\ d δ b & = \frac{1}{2} [\frac{d δ ϕ}{d S} \cdot (n_{1} n_{2} - n_{2} n_{1}) \cdot d ϕ - δ ϕ \cdot (n_{1} n_{2} - n_{2} n_{1}) \cdot \frac{d d ϕ}{d S}], \\ d δ χ & ​ = 0. \end{matrix}

Strains

The gradient of the current position of a point in the section with respect to the coordinate S is

\frac{\partial \hat{x}}{\partial S} = \frac{d x}{d S} + \frac{d f}{d S} S^{α} n_{α} + f S^{α} \frac{d n_{α}}{d S} + \frac{d w}{d S} ψ t + w ψ \frac{d t}{d S} .

We will keep only terms up to order $\frac{h}{ℓ}$ . We assume that $\frac{d f}{d S} ≪ \frac{1}{ℓ}$ , and—hence—the second term can be neglected. We also assume that $w ψ = O (\frac{h^{2}}{ℓ}),$ and—since $\frac{d t}{d S} = O (\frac{1}{ℓ})$ —we can neglect the last term as well. However, the warping function may vary rapidly near the ends where warping is constrained. Hence, $\frac{d w}{d S} ψ = O (\frac{h}{ℓ})$ and should be preserved. With those approximations,

\frac{\partial \hat{x}}{\partial S} \approx \frac{d x}{d S} + f S^{α} \frac{d n_{α}}{d S} + \frac{d w}{d S} ψ t = λ s + f S^{α} ϵ_{α}^{β} (- b_{β} t + b n_{β}) + χ ψ t .

The gradients with respect to $S^{α}$ are

\frac{\partial \hat{x}}{\partial S^{α}} = f n_{α} + w \frac{\partial ψ}{\partial S^{α}} t .

Correspondingly, in the original configuration

\begin{matrix} \frac{\partial \hat{X}}{\partial S} = ​ & \frac{d X}{d S} + S^{α} \frac{d N_{α}}{d S} = T + S^{α} ϵ_{α}^{β} (- B_{β} T + B N_{β}), \\ \frac{\partial \hat{X}}{\partial S^{α}} & ​ = N_{α} . \end{matrix}

The above relations are readily inverted to yield

\begin{matrix} \frac{\partial S}{\partial \hat{X}} & ​ = {(1 - S^{γ} ϵ_{γ}^{δ} B_{δ})}^{- 1} T, \\ \frac{\partial S^{α}}{\partial \hat{X}} & ​ = N_{α} - {(1 - S^{γ} ϵ_{γ}^{δ} B_{δ})}^{- 1} S^{β} ϵ_{β}^{α} B T . \end{matrix}

The deformation gradient then becomes

\begin{matrix} F = ​ & \frac{\partial \hat{x}}{\partial S} \frac{\partial S}{\partial \hat{X}} + \frac{\partial \hat{x}}{\partial S^{α}} \frac{\partial S^{α}}{\partial \hat{X}} \\ = & {(1 - S^{γ} ϵ_{γ}^{δ} B_{δ})}^{- 1} [λ s T + (- f S^{α} ϵ_{α}^{β} b_{β} + χ ψ) t T + f S^{α} ϵ_{α}^{β} b n_{β} T \\ + (1 - S^{α} ϵ_{α}^{β} B_{β}) f n_{α} N_{α} + (1 - S^{α} ϵ_{α}^{β} B_{β}) w \frac{d ψ}{d S^{η}} t N_{η} \\ - S^{α} ϵ_{α}^{β} B f n_{β} T - S^{α} ϵ_{α}^{β} B w \frac{\partial ψ}{\partial S^{β}} t T] . \end{matrix}

We define the initial length ratio R as

R = 1 - S^{α} ϵ_{α}^{β} B_{β} .

In the expression enclosed within square brackets above, terms of order $h^{2} / ℓ^{2}$ can again be neglected. However, it is assumed that the section may have low resistance to torsion and that, hence, the warping and twist may be large. This is particularly true for thin walled open sections. Hence, we obtain:

\begin{matrix} F = R^{- 1} { & λ sT + [- S^{α} ϵ_{α}^{β} (f b_{β} + B w \frac{\partial ψ}{\partial S^{β}}) + χ ψ] tT + \\ f S^{α} ϵ_{α}^{β} (b - B) n_{β} T + w \frac{\partial ψ}{\partial S^{β}} {tN}_{β} + R f n_{β} N_{β}} . \end{matrix}

We calculate the components of $F$ in a corotational system with the approximation $t \cdot s = 1$ . This provides

\begin{matrix} F_{11} & = t \cdot F \cdot T = R^{- 1} [λ - S^{β} ϵ_{β}^{α} (f b_{α} + B w \frac{\partial ψ}{\partial S^{α}}) + χ ψ], \\ F_{(α + 1) 1} & = n_{α} \cdot F \cdot T = R^{- 1} (λ γ_{α} + ϵ_{β}^{α} f S^{β} (b - B)), \\ F_{1 (β + 1)} & = t \cdot F \cdot N_{β} = R^{- 1} w \frac{d ψ}{d S^{β}}, \\ F_{(α + 1) (β + 1)} & = n_{α} \cdot F \cdot N_{β} = f δ_{α β} . \end{matrix}

We again neglect all terms of order $h^{n} / ℓ^{n}$ for $n \geq 2$ , except the term involving $B w$ . The equations then simplify to

\begin{matrix} F_{11} & ​ = λ - S^{β} ϵ_{β}^{α} (f b_{α} + B w \frac{\partial ψ}{\partial S^{α}} - λ B_{α}) + χ ψ, \\ F_{(α + 1) 1} & ​ = λ γ_{α} + ϵ_{β}^{α} f S^{β} (b - B), \\ F_{1 (β + 1)} & ​ = w \frac{d ψ}{d S^{β}}, \\ F_{(α + 1) (β + 1)} & ​ = f δ_{α β} . \end{matrix}

Consistent with traditional shell and beam theories, we slightly adapt the term involving the initial curvature $B_{α}$ —instead of multiplying it with $λ$ , we multiply it with f. Such a change does not significantly increase the error in the curvature calculation, since we do not properly account for the initial curvature in the volume integration anyway. Hence, we find for $F_{11}$ :

F_{11} = λ - S^{β} ϵ_{β}^{α} (f (b_{α} - B_{α}) + B w \frac{\partial ψ}{\partial S^{α}}) + χ ψ .

We now make a multiplicative decomposition of $F$ into a “stretch” part $F^{s}$ and a “distortion” part $F^{d}$ such that $F = F^{d} \cdot F^{s}$ .

For $F^{s}$ we choose

F_{11}^{s} = λ, F_{(α + 1) 1}^{s} = F_{1 (β + 1)}^{s} = 0, and F_{(α + 1) (β + 1)}^{s} = f δ_{α β};

and, hence, $F^{d}$ is

\begin{matrix} F_{11}^{d} & = λ^{- 1} F_{11} = 1 - S^{β} ϵ_{β}^{α} λ^{- 1} (f (b_{α} - B_{α}) + B w \frac{\partial ψ}{\partial S^{α}}) + λ^{- 1} χ ψ, \\ F_{(α + 1) 1}^{d} & = λ^{- 1} F_{(α + 1) 1} = γ_{α} + ϵ_{β}^{α} f λ^{- 1} S^{β} (b - B), \\ F_{1 (β + 1)}^{d} & = f^{- 1} F_{1 (β + 1)} = f^{- 1} w \frac{\partial ψ}{\partial S^{β}}, \\ F_{(α + 1) (β + 1)}^{d} & = f^{- 1} F_{(α + 1) (β + 1)} = δ_{α β} . \end{matrix}

Since $F^{s}$ is a diagonal tensor, the logarithmic “stretch” strains are immediately available as

e_{11}^{s} = \ln λ, γ_{(α + 1) 1}^{s} = 0, γ_{1 (α + 1)}^{s} = 0, e_{(α + 1) (β + 1)}^{s} = δ_{α β} \ln (f) .

Since the “distortional” strains are small, we obtain them from $F^{d}$ with the Green-Lagrange formula:

E^{d} = \frac{1}{2} (F^{d}^{T} \cdot F^{d} - I) .

For the components this yields

\begin{matrix} E_{11}^{d} = ​ & - S^{β} ϵ_{α}^{β} λ^{- 1} (f (b_{α} - B_{α}) + B w \frac{\partial ψ}{\partial S^{α}}) + λ^{- 1} χ ψ \\ - S^{β} ϵ_{α}^{β} (b_{α} - B_{α} + B w \frac{\partial ψ}{\partial S^{α}}) λ^{- 1} χ ψ + γ_{α} ϵ_{α}^{β} f λ^{- 1} S^{β} (b - B) \\ + \frac{1}{2} S^{γ} ϵ_{γ}^{α} S^{β} ϵ_{β}^{δ} (b_{α} - B_{α} + B w \frac{\partial ψ}{\partial S^{α}}) (b_{δ} - B_{δ} + B w \frac{\partial ψ}{\partial S^{δ}}) \\ + \frac{1}{2} λ^{- 2} χ^{2} ψ^{2} + \frac{1}{2} γ_{α} γ_{α} + \frac{1}{2} f^{2} λ^{- 2} S^{α} S^{α} {(b - B)}^{2}, \end{matrix}

\begin{matrix} E_{(α + 1) 1} = ​ & \frac{1}{2} [f^{- 1} w \frac{\partial ψ}{\partial S^{α}} + γ_{α} + ϵ_{α}^{β} f λ^{- 1} S^{β} (b - B) \\ - λ^{- 1} w \frac{\partial ψ}{\partial S^{α}} ϵ_{β}^{γ} (b_{γ} - B_{γ} + B w \frac{\partial ψ}{\partial S^{γ}}) + f^{- 1} w \frac{\partial ψ}{\partial S^{α}} λ^{- 1} χ ψ], \\ E_{(α + 1) (β + 1)} = ​ & \frac{1}{2} f^{- 2} w^{2} \frac{\partial ψ}{\partial S^{α}} \frac{\partial ψ}{\partial S^{β}} . \end{matrix}

Note that these strains are small. Since the various terms have different dependence on the cross-sectional coordinates, this leads to the conditions

\begin{matrix} S^{β} ϵ_{α}^{β} f λ^{- 1} (b_{α} - B_{α}) & ​ ≪ 1, \\ λ^{- 1} χ ψ & ​ ≪ 1, \\ γ_{α} & ​ ≪ 1, \\ f^{- 1} w \frac{\partial ψ}{\partial S^{α}} + ϵ_{α}^{β} f λ^{- 1} S^{β} (b - B) & ​ ≪ 1. \end{matrix}

The last condition will generally require that both w and $b - B$ are small, since $\partial ψ / \partial S^{α}$ will not be proportional to $ϵ_{α}^{β} S^{β}$ . However, for thin walled open section beams the proportionality is approximately satisfied and, therefore, we obtain in that case

w \frac{\partial ψ}{\partial S^{α}} \approx - f^{2} λ^{- 1} ϵ_{α}^{β} S^{β} (b - B) .

Note that within the desired accuracy this equation even holds for other sections: in that case both the right-hand and left-hand side are very small. Substitution in the expression for the Green-Lagrange strain yields the (small) distortional strains:

\begin{matrix} e_{11}^{d} & ​ = - f λ^{- 1} S^{β} ϵ_{β}^{α} (b_{α} - B_{α}) + λ^{- 1} χ ψ + \frac{1}{2} f^{2} λ^{- 2} S^{α} S^{α} (b^{2} - B^{2}), \\ γ_{(α + 1) 1}^{d} & ​ = γ_{α} + ϵ_{β}^{α} f λ^{- 1} S^{β} (b - B) + f^{- 1} w \frac{\partial ψ}{\partial S^{α}}, \\ γ_{1 (α + 1)}^{d} & ​ = γ_{(α + 1) 1}^{d}, \\ e_{(α + 1) (β + 1)}^{d} & ​ = + \frac{1}{2} f^{2} λ^{- 2} S^{α} S^{β} {(b - B)}^{2} . \end{matrix}

The total strains are obtained simply by addition as

\begin{matrix} e_{11} & ​ = \ln λ - f λ^{- 1} S^{β} ϵ_{β}^{α} (b_{α} - B_{α}) + λ^{- 1} χ ψ + \frac{1}{2} f^{2} λ^{- 2} S^{α} S^{α} (b^{2} - B^{2}), \\ γ_{(α + 1) 1} & ​ = γ_{α} + ϵ_{β}^{α} f λ^{- 1} S^{β} (b - B) + f^{- 1} w \frac{\partial ψ}{\partial S^{α}}, \\ γ_{1 (α + 1)} & ​ = γ_{(α + 1) 1}, \\ e_{(α + 1) (β + 1)} & ​ = δ_{α β} \ln (f) + \frac{1}{2} f^{2} λ^{- 2} S^{α} S^{β} {(b - B)}^{2} . \end{matrix}

We assume that there are no stresses in the $(α + 1) (β + 1)$ directions. Hence, the strains in these directions do not contribute to virtual work and do not need to be considered any further.

It is useful to split the total warping w in two parts: a part due to “free” warping $w_{f}$ minus a part due to warping prevention $w_{p}$ : $w = w_{f} - w_{p}$ . We assume that the warping function $ψ$ is chosen such that the free warping is related to the twist with the relation

w_{f} = f^{2} λ^{- 1} (b - B) and, hence, w_{p} = f^{2} λ^{- 1} (b - B) - w .

This makes it possible to write the expression for $γ_{(α + 1) 1}$ as

γ_{(α + 1) 1} = γ_{α} + f λ^{- 1} (ϵ_{β}^{α} S^{β} + \frac{\partial ψ}{\partial S^{α}}) (b - B) - f^{- 1} w_{p} \frac{\partial ψ}{\partial S^{α}} .

It is desirable to choose the cross-sectional resultants such that they are completely uncoupled. In addition, we assume that the axial strain variation across the section due to the second-order term in the twist is not significant. Therefore, we consider only the average axial strain due to the second-order twist term. Hence, we introduce the average axial strain

\bar{e} = \ln λ - f λ^{- 1} S_{c}^{β} ϵ_{β}^{α} (b_{α} - B_{α}) + λ^{- 1} χ \bar{ψ} + \frac{1}{2} (\frac{I_{p}}{A} + S_{c}^{α} S_{c}^{α}) f^{2} λ^{- 2} (b^{2} - B^{2}),

where $S_{c}^{β}$ are the coordinates of the centroid, $I_{p}$ is the polar moment of inertia, and $\bar{ψ}$ is the average value of the warping function:

\bar{ψ} = \frac{1}{A} \int_{A} ψ d A .

Similarly, we introduce the average shear strain

{\bar{γ}}_{α} = γ_{α} + f λ^{- 1} (ϵ_{β}^{α} S_{c}^{β} + \frac{1}{A} \int_{A} \frac{\partial ψ}{\partial S^{α}} d A) (b - B) - f^{- 1} w_{p} \frac{1}{A} \int_{A} \frac{\partial ψ}{\partial S^{α}} d A .

This last expression can be simplified by the introduction of the shear center coordinates $S_{s}^{β}$ , which are related to the warping function by

ϵ_{β}^{α} S_{s}^{β} = ϵ_{β}^{α} S_{c}^{β} + \frac{1}{A} \int_{A} \frac{\partial ψ}{\partial S^{α}} d A, or S_{s}^{β} = S_{c}^{β} + \frac{1}{A} \int_{A} ϵ_{β}^{α} \frac{\partial ψ}{\partial S^{α}} d A .

This yields

{\bar{γ}}_{α} = γ_{α} + f λ^{- 1} ϵ_{β}^{α} S_{s}^{β} (b - B) + f^{- 1} w_{p} ϵ_{β}^{α} (S_{c}^{β} - S_{s}^{β}) .

Note that the average value is in fact the value at the shear center if warping prevention is absent ( $w_{p} = 0$ ). However, for full warping prevention ( $w_{p} = f^{2} λ^{- 1} (b - B)$ ) the average value corresponds to the value obtained at the centroid.

Instead of the original warping function $ψ (S^{α})$ we now introduce a modified warping function $Ω (S^{α})$ related to $ψ$ by

Ω (S^{α}) \overset{def}{=} ψ (S^{α}) - \bar{ψ} .

This function in fact represents the classical definition of a warping function with an area weighted average of zero. The average value $\bar{ψ}$ can be obtained from the classical warping function with the condition that $ψ (0) = 0$ :

\bar{ψ} = ψ (S^{α}) - Ω (S^{α}) = ψ (0) - Ω (0) = - Ω (0) = - Ω_{\circ} .

The expression for the average axial strain then becomes

\bar{e} = \ln λ - f λ^{- 1} S_{c}^{β} ϵ_{β}^{α} (b_{α} - B_{α}) - λ^{- 1} χ Ω_{\circ} + \frac{1}{2} (\frac{I_{p}}{A} + S_{c}^{α} S_{c}^{α}) f^{2} λ^{- 2} (b^{2} - B^{2}),

and the location of the shear center is

S_{s}^{β} = S_{c}^{β} + \frac{1}{A} \int_{A} ϵ_{β}^{α} \frac{\partial Ω}{\partial S^{α}} d A .

The strains can, thus, be written in the form

\begin{matrix} e_{11} & ​ = \bar{e} - f λ^{- 1} (S^{β} - S_{c}^{β}) ϵ_{β}^{α} (b_{α} - B_{α}) + λ^{- 1} χ Ω, \\ γ_{(α + 1) 1} & ​ = {\bar{γ}}_{α} + f λ^{- 1} (ϵ_{β}^{α} (S^{β} - S_{s}^{β}) + \frac{\partial Ω}{\partial S^{α}}) (b - B) + f^{- 1} w_{p} (ϵ_{β}^{α} (S_{s}^{β} - S_{c}^{β}) - \frac{\partial Ω}{\partial S^{α}}) . \end{matrix}

The second term in the expression for $γ_{(α + 1) 1}$ is proportional to the shear strain field for pure elastic torsion:

γ_{α}^{t} (S^{γ}) \overset{def}{=} ϵ_{β}^{α} (S^{β} - S_{s}^{β}) + \frac{\partial Ω}{\partial S^{α}} .

We use this definition to eliminate the gradient of $Ω$ from the expression for the shear strain, which yields as final expression for the strains

\begin{matrix} e_{11} & ​ = \bar{e} - f λ^{- 1} (S^{β} - S_{c}^{β}) ϵ_{β}^{α} (b_{α} - B_{α}) + λ^{- 1} χ Ω, \\ γ_{(α + 1) 1} & ​ = {\bar{γ}}_{α} + f λ^{- 1} γ_{α}^{t} (b - B) + f^{- 1} w_{p} (ϵ_{β}^{α} (S^{β} - S_{c}^{β}) - γ_{α}^{t}) . \end{matrix}

Virtual work

Since it was assumed that there are no stresses in the $(α + 1) (β + 1)$ directions, the virtual work contribution is

\begin{matrix} δ Π & ​ = \int_{V} (σ_{11} δ e_{11} + τ_{(α + 1) 1} δ γ_{(α + 1) 1}) d V \\ ​ = \int_{ℓ} d ℓ \int_{A} (σ_{11} δ e_{11} + τ_{(α + 1) 1} δ γ_{(α + 1) 1}) d A . \end{matrix}

The strain variations are obtained by linearization of the expressions for the strains

\begin{matrix} δ e_{11} & ​ = δ \bar{e} - f λ^{- 1} (S^{β} - S_{c}^{β}) ϵ_{β}^{α} δ b_{α} + λ^{- 1} Ω δ χ, \\ δ γ_{(α + 1) 1} & ​ = δ {\bar{γ}}_{α} + f λ^{- 1} γ_{α}^{t} δ b + f^{- 1} (ϵ_{β}^{α} (S^{β} - S_{c}^{β}) - γ_{α}^{t}) δ w_{p}, \end{matrix}

where all terms of the order of the “distortional” strains have been neglected. From the expressions for the average axial strain and average shear strain, we obtain the average axial and shear strain variations as

\begin{matrix} δ \bar{e} & ​ = λ^{- 1} δ λ - λ^{- 1} f S_{c}^{β} ϵ_{β}^{α} δ b_{α} - λ^{- 1} Ω_{\circ} δ χ + f^{2} λ^{- 2} (\frac{I_{p}}{A} + S_{c}^{α} S_{c}^{α}) b δ b, \\ δ {\bar{γ}}_{α} & ​ = δ γ_{α} + λ^{- 1} f S_{s}^{β} ϵ_{β}^{α} δ b + f^{- 1} ϵ_{β}^{α} (S_{c}^{β} - S_{s}^{β}) δ w_{p}, \end{matrix}

where $δ w_{p} = f^{2} λ^{- 1} δ b - δ w$ and again all terms of the order of the “distortional” strains have been neglected.

We now introduce the generalized section forces as defined below:


Axial force	$F = \int_{A} σ_{11} d A$
Shear forces	$F_{α} = \int_{A} τ_{(α + 1) 1} d A$
Bending moments	$M^{α} = \int_{A} - f (S^{β} - S^{β}_{c}) ϵ_{β}^{α} σ_{11} d A$
Twisting moment	$M_{t} = \int_{A} f γ_{α}^{t} τ_{(α + 1) 1} d A$
Warping moment	$M_{w} = \int_{A} f (ϵ_{β}^{α} (S^{β} - S^{β}_{c}) - γ_{α}^{t}) τ_{(α + 1) 1} d A$
Bimoment	$W = \int_{A} Ω σ_{11} d A$

which transforms the virtual work contribution into

δ Π = \int_{ℓ} (F δ \bar{e} + F_{α} δ {\bar{γ}}_{α} + M^{α} λ^{- 1} δ b_{α} + M_{t} λ^{- 1} δ b + M_{w} f^{- 2} δ w_{p} + W λ^{- 1} δ χ) d ℓ .

Observe that the total torque T relative to the centroid of the section is the sum of the twisting moment and the warping moment:

T = \int_{A} f ϵ_{β}^{α} (S^{β} - S_{c}^{β}) τ_{(α + 1) 1} d A = M_{t} + M_{w} .

The rate of change of virtual work

To obtain the rate of change of virtual work, we first transform the integrations in the virtual work equation to the original volume such that

δ Π = \int_{ℓ^{\circ}} d ℓ^{\circ} \int_{A^{\circ}} (σ_{11} f^{2} λ δ e_{11} + τ_{(α + 1) 1} f^{2} λ δ γ_{(α + 1) 1}) d A^{\circ} .

The strain variations relative to the original state are then

\begin{matrix} f^{2} λ δ e_{11} & ​ = f^{2} λ δ \bar{e} - f^{3} (S^{β} - S_{c}^{β}) ϵ_{β}^{α} δ b_{α} + f^{2} Ω δ χ, \\ f^{2} λ δ γ_{(α + 1) 1} & ​ = f^{2} λ δ {\bar{γ}}_{α} + f^{3} γ_{α}^{t} δ b + f λ (ϵ_{β}^{α} (S^{β} - S_{c}^{β}) - γ_{α}^{t}) δ w_{p} . \end{matrix}

The strain variations relative to the original state are

\begin{matrix} f^{2} λ δ \bar{e} & ​ = f^{2} δ λ - f^{3} S_{c}^{β} ϵ_{β}^{α} δ b_{α} - f^{2} Ω_{\circ} δ χ + f^{4} λ^{- 1} (\frac{I_{p}}{A} + S_{c}^{α} S_{c}^{α}) b δ b, \\ f^{2} λ δ {\bar{γ}}_{α} & ​ = f^{2} λ δ γ_{α} + f^{3} S_{s}^{β} ϵ_{β}^{α} δ b + f λ ϵ_{β}^{α} (S_{c}^{β} - S_{s}^{β}) δ w_{p} . \end{matrix}

The rate of change of virtual work is then

\begin{matrix} d δ Π & = \int_{ℓ^{\circ}} d ℓ^{\circ} \int_{A^{\circ}} [f^{2} λ (d σ_{11} δ e_{11} + d τ_{(α + 1) 1} δ γ_{(α + 1) 1} \\ + σ_{11} (2 f d f δ λ - 3 f^{2} S^{β} ϵ_{β}^{α} d f δ b_{α} + 2 f ψ d f δ χ + f^{4} λ^{- 1} (\frac{I_{p}}{A} + S_{c}^{α} S_{c}^{α}) d b δ b \\ + f^{2} d δ λ - f^{3} S^{β} ϵ_{β}^{α} d δ b_{α} + f^{2} ψ d δ χ) \\ + τ_{(α + 1) 1} (f^{2} d λ δ γ_{α} + 2 λ f d f δ γ_{α} + 3 f^{2} S^{β} ϵ_{β}^{α} δ b d f + f \frac{\partial ψ}{\partial S^{α}} d λ δ w \\ + λ \frac{\partial ψ}{\partial S^{α}} d f δ w + λ f^{2} d δ γ_{α} + f^{3} S^{β} ϵ_{β}^{α} d δ b + λ f \frac{\partial ψ}{\partial S^{α}} d δ w)] d A ​^{\circ}, \end{matrix}

where we have neglected terms of order b. The changes in stress follow from the constitutive law

{\begin{matrix} d σ_{11} \\ d τ_{21} \\ d τ_{31} \end{matrix}} = [\begin{matrix} D_{11} & D_{12} & D_{13} \\ D_{21} & D_{22} & D_{23} \\ D_{31} & D_{32} & D_{33} \end{matrix}] {\begin{matrix} d e_{11} \\ d γ_{21} \\ d γ_{31} \end{matrix}} .

We approximate $d e_{11}$ and $d γ_{11}$ with the same relations as used for virtual work:

\begin{matrix} d e_{11} & ​ = d \bar{e} - f λ^{- 1} (S^{β} - S_{c}^{β}) ϵ_{β}^{α} d b_{α} + λ^{- 1} Ω d χ, \\ d γ_{(α + 1) 1} & ​ = d {\bar{γ}}_{α} + f λ^{- 1} γ_{α}^{t} d b + f^{- 1} (ϵ_{β}^{α} (S^{β} - S_{c}^{β}) - γ_{α}^{t}) d w_{p}, \end{matrix}

and for the average strain rates

\begin{matrix} d \bar{e} & ​ = λ^{- 1} d λ - λ^{- 1} f S_{c}^{β} ϵ_{β}^{α} d b_{α} - λ^{- 1} Ω_{\circ} d χ + f^{2} λ^{- 2} (\frac{I_{p}}{A} + S_{c}^{α} S_{c}^{α}) b d b, \\ d {\bar{γ}}_{α} & ​ = d γ_{α} + λ^{- 1} f S_{s}^{β} ϵ_{β}^{α} d b + f^{- 1} ϵ_{β}^{α} (S_{c}^{β} - S_{s}^{β}) d w_{p} . \end{matrix}

We now neglect all terms that involve the product of a stress tensor; a variation in curvature, twist, or warping; and a change in axial strain of the cross-section. Using the previously obtained expressions for the second variations in $λ, b_{α}, b, γ_{α}$ , and $χ$ and transforming the results into the current state, this provides

\begin{matrix} d δ Π & = \int_{ℓ} d ℓ \int_{A} [d σ_{11} δ e_{11} + d τ_{(α + 1) 1} δ γ_{(α + 1) 1} \\ + σ_{11} (2 λ^{- 1} δ λ f^{- 1} d f + λ^{- 1} d δ λ - λ^{- 1} f S^{β} ϵ_{β}^{α} d δ b_{α} + f^{2} λ^{- 2} (\frac{I_{p}}{A} + S_{c}^{α} S_{c}^{α}) d b δ b) \\ + τ_{(α + 1) 1} (d δ γ_{α} + λ^{- 1} f S^{β} ϵ_{β}^{α} d δ b + f^{- 1} \frac{\partial ψ}{\partial S^{α}} d δ w)] d A . \end{matrix}

The incremental moments, forces, etc. are defined as

\begin{matrix} d F & ​ = \int_{A} d σ_{11} d A, \\ d F_{α} & ​ = \int_{A} d τ_{(α + 1) 1} d A, \\ d M^{α} & ​ = \int_{A} - f (S^{β} - S_{c}^{β}) ϵ_{β}^{α} d σ_{11} d A, \\ d M_{t} & ​ = \int_{A} f γ_{α}^{t} d τ_{(α + 1) 1} d A, \\ d M_{w} & ​ = \int_{A} f (ϵ_{β}^{α} (S^{β} - S_{c}^{β}) - γ_{α}^{t}) d τ_{(α + 1) 1} d A, \\ d W & ​ = \int_{A} Ω d σ_{11} d A . \end{matrix}

For the determination of the initial stress stiffness, we assume that the second variations of the warping function and its derivative vanish: $d δ w = d δ χ = 0$ . Consequently,

d δ w_{p} = f^{2} λ^{- 1} d δ b,

and, therefore, only the torque relative to the centroid plays a role in the initial stress contribution to the rate of change of virtual work:

\begin{matrix} d δ Π = \int_{ℓ} & (d F δ \bar{e} + d F_{α} δ {\bar{γ}}_{α} + d M^{α} λ^{- 1} δ b_{α} + d M_{t} λ^{- 1} δ b + d M_{w} f^{- 2} δ w_{p} + d W λ^{- 1} δ χ \\ + F d δ \bar{e} + F_{α} d δ {\bar{γ}}_{α} + M^{α} λ^{- 1} d δ b_{α} + T λ^{- 1} d δ b) d ℓ . \end{matrix}

The second variations of $\bar{e}$ and ${\bar{γ}}_{α}$ contain contributions of the second variation in curvature and twist. These can be separated out by defining the bending moments and the torque relative to the origin of the cross-sectional coordinate system:

\begin{matrix} {\bar{M}}^{α} & = M^{α} + \int_{A} f S_{c}^{β} ϵ_{β}^{α} σ_{11} d A = M^{α} + f S_{c}^{β} ϵ_{β}^{α} F, \\ \bar{T} & = T + \int_{A} f S_{c}^{β} ϵ_{β}^{α} τ_{(α + 1) 1} d A = T + f S_{c}^{β} ϵ_{β}^{α} F_{α} . \end{matrix}

The expression for the rate of change of virtual work then takes the form

\begin{matrix} d δ Π = \int_{ℓ} & (d F δ \bar{e} + d F_{α} δ {\bar{γ}}_{α} + d M^{α} λ^{- 1} δ b_{α} + d M_{t} λ^{- 1} δ b + d M_{w} f^{- 2} δ w_{p} + d W λ^{- 1} δ χ + \\ F (λ^{- 1} d δ λ + 2 λ^{- 1} δ λ f^{- 1} d f + f^{2} λ^{- 2} (\frac{I_{p}}{A} + S_{c}^{α} S_{c}^{α}) d b δ b) + \\ F_{α} d δ γ_{α} + {\bar{M}}^{α} λ^{- 1} d δ b_{α} + \bar{T} λ^{- 1} d δ b) d ℓ . \end{matrix}

We propose to make the “cross-section size” a function of the stretch $λ$ . For the thermal stretch of the section we use isotropic expansion $(f_{θ} = λ_{θ}),$ and for the “mechanical” stretch of the section we assume an effective Poisson's ratio $ν (f_{m} = λ_{m}^{- ν})$ . The total cross-sectional stretch is

f = f_{m} f_{θ} = λ_{m}^{- ν} λ_{θ} = λ^{- ν} λ_{θ}^{1 + ν},

so that

d f = - ν λ^{- 1 - ν} λ_{θ}^{1 + ν} d λ .

Using this in the above expression for the rate of change of virtual work, we find

\begin{matrix} d δ Π = \int_{ℓ} & (d F δ \bar{e} + d F_{α} δ {\bar{γ}}_{α} + d M^{α} λ^{- 1} δ b_{α} + d M_{t} λ^{- 1} δ b + d M_{w} f^{- 2} δ w_{p} + d W λ^{- 1} δ χ + \\ F (λ^{- 1} d δ λ - 2 ν λ^{- 2} d λ δ λ + f^{2} λ^{- 2} (\frac{I_{p}}{A} + S_{c}^{α} S_{c}^{α}) d b δ b) + \\ F_{α} d δ γ_{α} + {\bar{M}}^{α} λ^{- 1} d δ b_{α} + \bar{T} λ^{- 1} d δ b) d ℓ . \end{matrix}

This expression is symmetric if the material tensor $[D]$ is symmetric.

Section integration

The formulation presented in the previous pages is valid for all possible beam types. However, different classes of beams will result in different final formulations. We consider three different classes of beams:

Beams in which warping may be constrained. These beams generally have an open, thin walled section reinforced with some relatively solid parts or some relatively small closed cells and have a torsional constant that is considerably smaller than the polar moment of inertia. Hence, in the elastic range, the warping can be large, and warping prevention at the ends can contribute significantly to the torsional rigidity of the beam. In this case both the torsional shear stresses and the axial warping stresses can be of the same order of magnitude as the stresses due to axial forces and bending moments, and the complete theory must be used.
Beams in which warping is unconstrained. These beams generally have a solid section or a closed, thin walled section and have a torsional constant that is of the same order of magnitude as the polar moment of inertia of the section. Hence, in the elastic range the warping is rather small, and it is assumed that warping prevention at the ends can be neglected. The axial warping stresses are assumed to be negligible, but the torsional shear stresses are assumed to be of the same order of magnitude as the stresses due to axial forces and bending moments. In this case the warping is dependent on the twist and can be eliminated as an independent variable, which leads to a considerably simplified formulation.
Beams in which warping constraints dominate the torsional stiffness. These beams generally have an open, thin walled section and have a torsional constant that is much smaller than the polar moment of inertia. In the elastic range the warping is likely to be large, and warping constraints are essential to provide torsional stiffness for the beam. In this case the axial stresses may be of the same order of magnitude as the stresses due to axial forces and bending moments, but the torsional shear stresses are relatively small. Hence, the warping can be coupled to the twist with a relatively stiff elastic constraint but cannot be eliminated because it must be possible to prevent warping at the nodes.

Examples of cross-sections for the last two classes will be derived after the general discussion of sections with and without warping prevention.

In Abaqus we neglect the effect of shear stresses due to transverse shear forces at individual material points. We will consequently always assume elastic behavior of the section in transverse shear, leading to the relations

F_{α} = f_{p} G A {\bar{γ}}_{α}, d F_{α} = f_{p} G A d {\bar{γ}}_{α},

where $F_{α}$ are the transverse shear forces, working at the shear center, and $f_{p}$ is the “slenderness compensation factor” used to prevent the shear stiffness from becoming too big in slender beams. The slenderness compensation factor is defined as

f_{p} = {(1 + 0.25 \frac{ℓ^{2} A}{12 I})}^{- 1} for first-order Timoshenko beam elements and

f_{p} = {(1 + 0.25 \times 10^{- 4} \frac{ℓ^{2} A}{12 I})}^{- 1} for all other beam elements,

where $ℓ$ is the length of the element and I is the larger of the moments of inertia $I_{11}$ and $I_{22}$ . Hence, the transverse shear terms do not need to be considered in any further detail.

The fact that the transverse shear forces are considered separately allows us to write

\begin{matrix} γ_{(α + 1) 1} & ​ = γ_{(α + 1) 1}^{F} + γ_{(α + 1) 1}^{M}, \\ τ_{(α + 1) 1} & ​ = τ_{(α + 1) 1}^{F} + τ_{(α + 1) 1}^{M}, \end{matrix}

where $γ_{(α + 1) 1}^{F}$ and $τ_{(α + 1) 1}^{F}$ are the strains and stresses due to transverse shear forces and $γ_{(α + 1) 1}^{M}$ and $τ_{(α + 1) 1}^{M}$ are the strains and stresses due to a twisting around the shear center. Substitution in the expressions for the twisting and warping moments yields

\begin{matrix} M_{t} & = \int_{A} f γ_{α}^{t} (τ_{(α + 1) 1}^{F} + τ_{(α + 1) 1}^{M}) d A = \int_{A} f γ_{α}^{t} τ_{(α + 1) 1}^{M} d A, \\ M_{w} & ​ = \int_{A} f (ϵ_{β}^{α} (S^{β} - S_{c}^{β}) - γ_{α}^{t}) (τ_{(α + 1) 1}^{F} + τ_{(α + 1) 1}^{M}) d A \\ ​ = ϵ_{β}^{α} (S_{s}^{β} - S_{c}^{β}) F_{α} + \int_{A} f (ϵ_{β}^{α} (S^{β} - S_{c}^{β}) - γ_{α}^{t}) τ_{(α + 1) 1}^{M} d A, \end{matrix}

where we used the fact that $γ_{α}^{t}$ was calculated based on application of a twisting moment around the shear center and, hence, does not do any work on the shear stresses due to transverse shear forces.

The warping function $Ω$ is assumed to be determined based on isotropic, homogeneous elastic behavior of the section in shear. For this case the elastic energy due to twist is

U_{t} = \frac{1}{2} G \int_{A} γ_{(α + 1) 1}^{M} γ_{(α + 1) 1}^{M} d A .

For twist without warping constraints $w_{p}$ vanishes and the energy per unit length of the beam is

U_{t} = \frac{1}{2} f^{2} λ^{- 2} {(b - B)}^{2} G \int_{A} γ_{α}^{t} γ_{α}^{t} d A = \frac{1}{2} f^{2} λ^{- 2} {(b - B)}^{2} G J,

where we have introduced the torsion integral J given by

J = \int_{A} γ_{α}^{t} γ_{α}^{t} d A .

With complete warping prevention the $w_{p} = f^{2} λ^{- 1} (b - B)$ and the energy per unit length of the beam is

U_{t} = \frac{1}{2} f^{2} λ^{- 2} {(b - B)}^{2} G \int_{A} ϵ_{β}^{α} (S^{β} - S_{c}^{β}) ϵ_{γ}^{α} (S^{γ} - S_{c}^{γ}) d A = \frac{1}{2} f^{2} λ^{- 2} {(b - B)}^{2} G I_{p},

where we have introduced the polar moment of inertia

I_{p} = \int_{A} (S^{β} - S_{c}^{β}) (S^{β} - S_{c}^{β}) d A = \int_{A} ϵ_{β}^{α} (S^{β} - S_{c}^{β}) ϵ_{γ}^{α} (S^{γ} - S_{c}^{γ}) d A .

For unconstrained warping $τ_{(α + 1) 1} = G γ_{α}^{t}$ . Since the twisting moment must be equal to

M_{t} = \int_{A} f (S^{β} - S_{c}^{β}) ϵ_{β}^{α} τ_{(α + 1) 1} d A = G \int_{A} f (S^{β} - S_{c}^{β}) ϵ_{β}^{α} γ_{α}^{t} d A,

and since

M_{t} = \int_{A} f γ_{α}^{t} τ_{(α + 1) 1} d A = G \int_{A} f γ_{α}^{t} γ_{α}^{t} d A,

it follows that

\int_{A} (S^{β} - S_{c}^{β}) ϵ_{β}^{α} γ_{α}^{t} d A = \int_{A} γ_{α}^{t} γ_{α}^{t} d A = J .

Beams with unconstrained warping

For this beam type, warping prevention is not taken into consideration. Hence we assume $w_{p} = 0$ . In addition we assume that axial strains due to warping can be neglected: $χ Ω \approx 0$ .

For the strains at a material point this yields

\begin{matrix} e_{11} & ​ = \bar{e} - (S^{β} - S_{c}^{β}) ϵ_{β}^{α} (f λ^{- 1} b_{α} - B_{α}), \\ γ_{(α + 1) 1}^{M} & ​ = f λ^{- 1} γ_{α}^{t} (b - B), \end{matrix}

and for the variations in the strains

\begin{matrix} δ e_{11} & ​ = δ \bar{e} - f λ^{- 1} (S^{β} - S_{c}^{β}) ϵ_{β}^{α} δ b_{α}, \\ δ γ_{(α + 1) 1}^{M} & ​ = f λ^{- 1} γ_{α}^{t} δ b . \end{matrix}

Substitution in the virtual work statement yields

δ Π = \int_{ℓ} (F δ \bar{e} + F_{α} δ {\bar{γ}}_{α} + M^{α} λ^{- 1} δ b_{α} + M_{t} λ^{- 1} δ b) d ℓ,

where the expressions derived earlier for F, $F_{α}$ , $M^{α}$ , and $M_{t}$ apply.

Although there is no warping prevention in the section, the warping moment $M_{w}$ does not vanish. From the expression obtained earlier follows

M_{w} = ϵ_{β}^{α} (S_{s}^{β} - S_{c}^{β}) F_{α} .

This yields for the torque around the centroid

T = ϵ_{β}^{α} (S_{s}^{β} - S_{c}^{β}) F_{α} + M_{t}

and for the torque around the origin

\bar{T} = ϵ_{β}^{α} S_{s}^{β} F_{α} + M_{t} .

For the rate of change of virtual work we obtain

\begin{matrix} d δ Π = \int_{ℓ} & (d F δ \bar{e} + d F_{α} δ {\bar{γ}}_{α} + d M^{α} λ^{- 1} δ b_{α} + d M_{t} λ^{- 1} δ b + \\ F (λ^{- 1} d δ λ - 2 ν λ^{- 2} d λ δ λ + f^{2} λ^{- 2} (\frac{I_{p}}{A} + S_{c}^{α} S_{c}^{α}) d b δ b) + \\ F_{α} d δ γ_{α} + {\bar{M}}^{α} λ^{- 1} d δ b_{α} + \bar{T} λ^{- 1} d δ b) d ℓ . \end{matrix}

Beams with elastic torsion and constrained warping

Consider the case that the shear stresses are defined from the shear strains by linear elastic response, with a constant shear modulus G:

τ_{(α + 1) 1}^{M} = G γ_{(α + 1) 1}^{M} = f λ^{- 1} G γ_{α}^{t} (b - B) + f^{- 1} G (ϵ_{β}^{α} (S^{β} - S_{c}^{β}) - γ_{α}^{t}) .

This allows us to write for the twisting and warping moments:

\begin{matrix} M_{t} & ​ = f^{2} λ^{- 1} G (b - B) \int_{A} γ_{α}^{t} γ_{α}^{t} d A + G w_{p} \int_{A} γ_{α}^{t} (ϵ_{β}^{α} (S^{β} - S_{c}^{β}) - γ_{α}^{t}) d A \\ ​ = f^{2} λ^{- 1} G J (b - B) \\ M_{w} & ​ = f^{2} λ^{- 1} G (b - B) \int_{A} γ_{α}^{t} (ϵ_{β}^{α} (S^{β} - S_{c}^{β}) - γ_{α}^{t}) d A \\ + G w_{p} \int_{A} (ϵ_{β}^{α} (S^{β} - S_{c}^{β}) - γ_{α}^{t}) (ϵ_{γ}^{α} (S^{γ} - S_{c}^{γ}) - γ_{α}^{t}) d A \\ ​ = G (I_{p} - J) w_{p} . \end{matrix}

Substitution of these expressions in the part of the virtual work equation related to torsion yields

δ Π_{t} = \int_{ℓ} (λ^{- 1} M_{t} δ b + f^{- 2} M_{w} δ w_{p}) d ℓ = \int_{ℓ} (f^{2} λ^{- 2} G J (b - B) δ b + f^{- 2} G (I_{p} - J) w_{p} δ w_{p}) d ℓ .

Note that the torque T around the centroid is

T = M_{t} + M_{w} = G J f^{2} λ^{- 1} (b - B) + G (I_{p} - J) w_{p} .

Hence, we can write virtual work in terms of the primary variables b and w:

\begin{matrix} δ Π_{t} & ​ = \int_{ℓ} (λ^{- 1} T δ b - f^{- 2} M_{w} δ w) d A \\ ​ = \int_{ℓ} [(f^{2} λ^{- 2} G J (b - B) + λ^{- 1} G (I_{p} - J) w_{p}) δ b - f^{- 2} G (I_{p} - J) w_{p} δ w] d ℓ . \end{matrix}

The complete virtual work equation has the form

δ Π = \int_{ℓ} (F δ \bar{e} + F_{α} δ {\bar{γ}}_{α} + λ^{- 1} M^{α} δ b_{α} + f^{2} λ^{- 2} G J b δ b + f^{- 2} G (I_{p} - J) w_{p} δ w_{p} + λ^{- 1} W δ χ) d ℓ,

where F, $F_{α}$ , $M^{α}$ and W are defined as before.

For the rate of change of virtual work we obtain similarly

\begin{matrix} d δ Π = \int_{ℓ} & (d F δ \bar{e} + d F_{α} δ {\bar{γ}}_{α} + d M^{α} λ^{- 1} δ b_{α} + f^{2} λ^{- 2} G J d b δ b + f^{- 2} G (I_{p} - J) d w_{p} δ w_{p} + \\ d W λ^{- 1} δ χ + F (λ^{- 1} d δ λ - 2 ν λ^{- 2} d λ δ λ + f^{2} λ^{- 2} (\frac{I_{p}}{A} + S_{c}^{α} S_{c}^{α}) d b δ b) + \\ F_{α} d δ γ_{α} + {\bar{M}}^{α} λ^{- 1} d δ b_{α} + \bar{T} λ^{- 1} d δ b) d ℓ \end{matrix}

with

\bar{T} = T + f ϵ_{β}^{α} S_{c}^{β} F_{α} .

We now discuss some specific section types incorporated in Abaqus.

Circular section

For this type of section, warping is absent. Hence,

Ω = \frac{\partial Ω}{\partial S^{α}} = 0.

Solid noncircular sections

Solid sections such as rectangles or trapezoids are included in this category. The warping function $ψ$ is a harmonic function and is subject to the condition that no shear stress component can act normal to the boundary of the cross-section. Although it is possible to determine the warping function $ψ$ in this manner, we choose to work in terms of the Saint-Venant's stress function because of its simplicity. Following standard procedures we normalize this function so that the (elastic) shear strains can be derived directly from it. We introduce the function $φ (S^{β})$ , which is differentiable in the cross-section and has the property that

γ_{α}^{t} = ϵ_{α}^{β} \frac{\partial φ}{\partial S^{β}} .

The stress function is determined by solving the differential equation of the form

\nabla^{2} φ = - 2 in A, with φ = 0 on S,

where S represents the boundary of the section. This boundary condition ensures that no shear stress component can act normal to the boundary.

For the solid noncircular sections this differential equation is solved numerically using a second-order isoparametric finite element. The torsional constant of the bar is then equal to twice the volume under the normalized stress function surface.

Closed, thin walled cross-sections

In this case we assume that the shear strain perpendicular to the section must vanish so that

(ϵ_{β}^{α} S^{β} + \frac{\partial Ω}{\partial S^{α}}) n_{α} = 0.

Since $t_{β} = - ϵ_{β}^{α} n_{α}$ , this yields

\frac{\partial Ω}{\partial n} = \frac{\partial Ω}{\partial S^{α}} n_{α} = t_{β} S^{β},

which can be identically satisfied anywhere. The normalized shear strain along the section is

\bar{γ} = (ϵ_{β}^{α} S^{β} + \frac{\partial Ω}{\partial S^{α}}) t_{α} = n_{β} S^{β} + \frac{\partial Ω}{\partial s},

where s is the distance along the section. Integrating this around the circumference yields

\oint \bar{γ} d s = \oint n_{α} S^{α} d s + \oint \frac{\partial Ω}{\partial s} d s = 2 A_{C},

where $A_{C}$ is the area enclosed by the section. With the assumption that the shear modulus is constant along the section, the total torsional elastic strain energy is

U_{T} = \oint \frac{1}{2} G {\bar{γ}}^{2} t d s,

where t is the wall thickness. Minimization of the energy with the constraint enforced with a Lagrange multiplier $μ$ yields

δ P = \oint (G \bar{γ} h - μ) δ \bar{γ} d s - (\oint \bar{γ} d s - 2 A_{C}) δ μ = 0.

Hence, $\bar{γ} = \frac{μ}{G h} and \frac{μ}{G} \oint \frac{1}{h} d s - A_{C} = 0,$ which combine to define

\bar{γ} = \frac{2 A_{C}}{h \oint \frac{1}{h} d s} .

This allows calculation of $\bar{γ}$ at any point in the section based on the section geometry.

Thin walled open sections

The most important sections that exhibit substantial warping are the thin walled open sections. For a single branch section we can conveniently express $Ω$ as a function of the coordinate s along the section and the coordinate z perpendicular to the section. A suitable approximation for $Ω$ is

Ω (s, z) = Ω_{c} (s) + z η (s),

where $Ω_{c} (s)$ is the value at the centerline of the wall. Minimization of the torsional elastic energy yields

\int_{s_{1}}^{s_{2}} \int_{- h / 2}^{h / 2} [(t_{α} ϵ_{β}^{α} (S^{β} - S_{s}^{β}) + \frac{d Ω_{c}}{d s} + z \frac{d η}{d s}) (\frac{d δ Ω_{c}}{d s} + z \frac{d δ η}{d s}) + (n_{α} ϵ_{β}^{α} S^{β} + η) δ η] d z d s = 0.

We now introduce $S_{0}^{β}$ , which is the position of the middle of the wall. With $t_{α} ϵ_{β}^{α} = n_{β}$ and $S^{β} = S_{0}^{β} + z n_{β}$ , and after carrying out the integration over z, the minimization condition simplifies to

\int_{s_{1}}^{s_{2}} [h (n_{β} (S_{0}^{β} - S_{s}^{β}) + \frac{d Ω_{c}}{d s}) \frac{d δ Ω_{c}}{d s} + h (- t_{β} (S_{0}^{β} - S_{s}^{β}) + η) δ η + \frac{1}{12} h^{3} (1 + \frac{d η}{d s}) \frac{d δ η}{d s}] d s = 0.

Clearly, the dominant terms are of order h. This yields the equations

\frac{d Ω_{c}}{d s} = - n_{β} (S_{0}^{β} - S_{s}^{β}), η = t_{β} (S_{0}^{β} - S_{s}^{β}),

so that $Ω (s, z)$ is

Ω (s, z) = \int_{0}^{s} - n_{β} (S_{0}^{β} - S_{s}^{β}) d s + t_{β} (S_{0}^{β} - S_{s}^{β}) z + Ω_{s},

where $Ω_{s}$ is the value of $Ω$ at the start of the integration over the section. Observe that

Ω_{c} (s) = \int_{0}^{s} - n_{β} (S_{0}^{β} - S_{s}^{β}) d s + Ω_{s}

represents the sectorial area between two points on the midsection relative to the shear center. Therefore, it readily follows that

\int_{0}^{s_{e n d}} h Ω_{c} S^{α} d s = 0,

so there is no coupling between twist and bending in the section for linear elastic material behavior. As was discussed before, $Ω_{s}$ must be chosen such that

\bar{Ω} = \int_{0}^{s_{e n d}} h Ω_{c} d s / \int_{0}^{s_{e n d}} h d s = 0,

which eliminates the coupling between twist and axial extension in the section for linear elasticity.

Note that coupling terms still exist but that they are incorporated in the generalized strain displacement relations. The coupling between twist and extension is governed by $Ω_{0}$ , the value of the warping function at the origin of the cross-sectional coordinate system. If the origin is on the section, this value can be evaluated properly. If the origin is not on the section (which means that the node is not connected to the section), we assume that $Ω_{0} = 0$ .

The torsion integral J is readily obtained as

\begin{matrix} J & ​ = \int_{s_{1}}^{s_{2}} \int_{- h / 2}^{h / 2} [{(n_{β} (S^{β} - S_{s}^{β}) + \frac{d Ω}{d s})}^{2} + {(- t_{β} (S^{β} - S_{s}^{β}) + \frac{d Ω}{d z})}^{2}] d z d s \\ = \int_{s_{1}}^{s_{2}} \int_{- h / 2}^{h / 2} {(2 z)}^{2} d z d s = \int_{s_{1}}^{s_{2}} \frac{1}{3} h^{3} d s, \end{matrix}

and the polar moment of inertia is given by the expression

I_{p} = \int_{s_{1}}^{s_{2}} h (S^{β} - S_{c}^{β}) (S^{β} - S_{c}^{β}) d s .

The above derivations cover single branch open sections. Multibranch open sections can be transformed into single branch open sections by connecting the end of one branch with the beginning of the next branch with a section that has thickness $h = 0$ . Such dummy sections do not yield any contribution to the area, the moments of inertia, or the torsion integral and, hence, have no influence on the results.