C0 Timoshenko Beam Element

Rotation increment

Let $^0R$ be the rotation matrix from the global coordinate system to the original beam local coordinate system. This is calculated using the unit vector along the axis of the beam (assumed to be the beam local x axis), the user-provided y_orientation that is orthonormal to the beam axis, and the cross product of the beam-axis and y_orientation. If the beam undergoes small deformation/rotation, the rotation matrix is not incrementally updated and the rotation matrix at $t+\Delta t$ ($^{t+\Delta t}R$) is assumed to be same as $^0R$.

If the beam undergoes finite deformation, a rotation increment ($\Delta R$) that transforms the beam from beam local configuration at $t$ to the beam local configuration at $t + \Delta t$ is calculated at each time step. This calculation is performed using Euler angles as described below.

Let $^{t+\Delta t} {u_1}^i$, $^{t+\Delta t} {u_2}^i$ and $^{t+\Delta t} {u_3}^i$ denote the translational displacements at node $i$ at $t+\Delta t$. These displacements are the difference between the nodal positions at time $t+\Delta t$ and $t$ and are in the global coordinate system. Here 1, 2 and 3 refers to the x, y and z coordinate, respectively. Similarly, let $^{t+\Delta t} {\theta_1}^i$, $^{t+\Delta t} {\theta_2}^i$ and $^{t+\Delta t} {\theta_3}^i$ denote the rotational displacements at node $i$ at $t+\Delta t$ in the global coordinate system. Let $^{t+\Delta t} \{ {\Delta u} \} = ^{t+\Delta t} {\{\Delta u\}}^1 - ^{t+\Delta t} {\{\Delta u\}}^0$ and $^{t+\Delta t} \{ {\Delta \theta} \} = ^{t+\Delta t} {\{\Delta \theta\}}^1 - ^{t+\Delta t} {\{\Delta \theta\}}^0$ be the differential displacements along the length of the beam.

These global displacements are converted to the beam local configuration at $t$ using the rotation matrix $^tR$ as follows: $$^{t+\Delta t} {\{ {\Delta u} \}}^{local} = {}^tR \;\; ^{t+\Delta t} \{ {\Delta u} \}$$

$$^{t+\Delta t} {\{ {\Delta \theta} \}}^{local} = {}^tR \;\; ^{t+\Delta t} \{ {\Delta \theta} \}$$

These beam local displacements are then used in the calculation of the Euler angles $\alpha$, $\beta$ and $\gamma$ as in Bathe and Bolourchi (1979). These Euler angles are used to calculate the rotation increment $\Delta R$. Using $\Delta R$, the rotation matrix at $t+\Delta t$ can be calculated as: $${}^{t+\Delta t}R = \Delta R \; {}^tR$$

The rotation matrix is incrementally updated if ComputeFiniteBeamStrain is used. If ComputeIncrementalBeamStrain is used, the rotation matrix is same as the initially computed rotation matrix ($^0R$).

Strain increment

To calculate the beam strain increment, the incremental beam displacements at time $t+\Delta t$ in the local configuration at $t+\Delta t$ are obtained using ${}^{t+\Delta t} R$ to transform the displacements from global to local configuration at $t+\Delta t$. For simplifying notation, let ${u_{j}}^i$ and ${\theta_{j}}^i$ be the translational and rotational displacements at node $i$ at $t+\Delta t$ in the beam local configuration at $t+\Delta t$. Interpolating the nodal DOFs along the axis of the beam using first order Lagrange shape functions gives $u_{j}(x) = u_{j}^0 \frac{x}{{}^0L} + u_{j}^1 \frac{{}^0L-x}{{}^0L}$ and $\theta_{j}(x) = \theta_{j}^0 \frac{x}{{}^0L} + \theta_{j}^1 \frac{{}^0L-x}{{}^0L}$. However, $u_j(x)$ only gives the translational displacement of the beam axis. The translational displacement at any point on the beam is obtained as follows: $${u_1}(x,y,z) = u_1(x) - \theta_3(x) y + \theta_2(x) z$$ $${u_2}(x,y,z) = u_2(x) - \theta_1(x) z$$ $${u_3}(x,y,z) = u_3(x) + \theta_1(x) y$$

The axial strain ($\epsilon_{11}$) and the shear strains ($\epsilon_{12}$ and $\epsilon_{13}$) are then obtained as follows: $$\epsilon_{11}(x,y,z) = \frac{\partial u_1(x,y,z)}{\partial x} = \frac{\partial u_1(x)}{\partial x} - \frac{\partial \theta_3(x)}{\partial x} y + \frac{\partial \theta_2(x)}{\partial x} z$$

$$\epsilon_{12}(x,y,z) = \frac{\partial u_1(x,y,z)}{\partial y} + \frac{\partial u_2(x,y,z)}{\partial x} = -\theta_3 + \frac{\partial u_2(x)}{\partial x} - \frac{\partial \theta_1(x)}{\partial x} z$$

$$\epsilon_{13}(x,y,z) = \frac{\partial u_1(x,y,z)}{\partial z} + \frac{\partial u_3(x,y,z)}{\partial x} = \theta_2 + \frac{\partial u_3(x)}{\partial x} + \frac{\partial \theta_1(x)}{\partial x} y$$

It should be noted here that using a linear interpolation of the rotational variables in the calculation of shear strain leads to shear locking of the beam. Using $\theta_1(x) = ({\theta_1}^0 + {\theta_1}^1)/2$ in the calculation of $\epsilon_{12}$ and $\epsilon_{13}$ ensures that the beam doesn't lock under shear deformation (Prathap and Bhashyam, 1982).

The above translational strain increments are functions of x, y and z. However, since the beam cross-section does not deform, these strains can be integrated over the cross-section to obtain strain increments as a function of only x. $$\epsilon_{1}(x) = \int_A \epsilon_{11}(x,y,z) dA= \frac{\partial u_1(x)}{\partial x} A - \frac{\partial \theta_3(x)}{\partial x} A_y + \frac{\partial \theta_2(x)}{\partial x} A_z$$

$$\epsilon_{2}(x) = \int_A \epsilon_{12}(x,y,z) dA= -\theta_3 A + \frac{\partial u_2(x)}{\partial x} A - \frac{\partial \theta_1(x)}{\partial x} A_z$$

$$\epsilon_{3}(x) = \int_A \epsilon_{13}(x,y,z) dA= \theta_2 A + \frac{\partial u_3(x)}{\partial x} A + \frac{\partial \theta_1(x)}{\partial x} A_y$$

where $A$ is the cross-sectional area, $A_y = \int_A y dA$ and $A_z = \int_A z dA$. $\epsilon_{1}(x)$, $\epsilon_{2}(x)$ and $\epsilon_{3}(x)$ are the axial and shear increments.

Apart from the translational strain increments, rotational strain increments also need to be calculated. $$\kappa_1(x) = \int_A (\epsilon_{13}(x,y,z) y - \epsilon_{12}(x,y,z) z) dA = \theta_2 A_y + \frac{\partial u_3(x)}{\partial x} A_y + \frac{\partial \theta_1(x)}{\partial x} I_x + \theta_3 A_z - \frac{\partial u_2(x)}{\partial x} A_z$$

$$\kappa_2(x) = \int_A \epsilon_{11}(x,y,z) z dA = \frac{\partial u_1(x)}{\partial x} A_z - \frac{\partial \theta_3(x)}{\partial x} Iyz + \frac{\partial \theta_2(x)}{\partial x} I_z$$

$$\kappa_3(x) = -\int_A \epsilon_{11}(x,y,z) y dA = -\frac{\partial u_1(x)}{\partial x} A_y + \frac{\partial \theta_3(x)}{\partial x} I_y - \frac{\partial \theta_2(x)}{\partial x} Iyz$$

where $I_y = \int_A y^2 dA$, $I_z = \int_A z^2 dA$, $I_x = \int_A (y^2 + z^2) dA$ and $Iyz = \int_A yz dA$. Note that $Iyz$ is assumed to be zero for simplicity and this assumption is valid for symmetric cross-sections such as square, rectangular and circular. $I_x$ is automatically calculated from $I_y$ and $I_z$ unless $I_x$ is provided as input by the user.

The above strain increments are calculated using the small strain assumption. If the large_strain option in ComputeIncrementalBeamStrain or ComputeFiniteBeamStrain is set to true, then the strains are calculated using the equations below:

$$\epsilon_{11}(x,y,z) = \frac{\partial u_1(x,y,z)}{\partial x} + \frac{1}{2} {\left[ {\left( \frac{\partial u_1(x,y,z)}{\partial x}\right)}^2 + {\left( \frac{\partial u_2(x,y,z)}{\partial x} \right)}^2 + {\left( \frac{\partial u_3(x,y,z)}{\partial x} \right)}^2 \right]}$$

$$\epsilon_{12}(x,y,z) = \frac{\partial u_1(x,y,z)}{\partial y} + \frac{\partial u_2(x,y,z)}{\partial x} + \frac{\partial u_1(x,y,z)}{\partial x} \frac{\partial u_1(x,y,z)}{\partial y} + \frac{\partial u_3(x,y,z)}{\partial x} \frac{\partial u_3(x,y,z)}{\partial y}$$

$$\epsilon_{13}(x,y,z) = \frac{\partial u_1(x,y,z)}{\partial z} + \frac{\partial u_3(x,y,z)}{\partial x} + \frac{\partial u_1(x,y,z)}{\partial x} \frac{\partial u_1(x,y,z)}{\partial z} + \frac{\partial u_2(x,y,z)}{\partial x} \frac{\partial u_2(x,y,z)}{\partial z}$$

Note that for the large_strain calculation, $A_y$, $A_z$ and $J$ are assumed to be zero for simplicity. This is again valid for symmetric cross-sections such as square, rectangular and circular cross-sections.

If displacement or rotational eigenstrains are provided as input, those eigenstrain increments would be subtracted from the total displacement and rotational strain increments to get the mechanical displacement and rotational strain increments.

Consistent mass matrix

If $A_y$, $A_z$ and $Iyz$ are zero in InertialForceBeam, then there is no coupling between the rotational and translational DOFs. The residual for the $j^{th}$ translational degree of freedom at node $i$ can be obtained as follows:

$$R_j^0 = \int_{0}^{{}^0L} \rho A \left(\frac{x}{{}^0L} {\ddot{u}_j}^0 + \frac{{}^0L - x}{{}^0L} {\ddot{u}_j}^1 \right) \frac{x}{{}^0L} dx$$ $$R_j^1 = \int_{0}^{{}^0L} \rho A \left(\frac{x}{{}^0L} {\ddot{u}_j}^0 + \frac{{}^0L - x}{{}^0L} {\ddot{u}_j}^1 \right) \frac{{}^0L - x}{{}^0L} dx$$ where $\rho$ is the density of the beam, which is assumed to be constant through the length of the beam, and ${\ddot{u}_j}^i$ is the translational acceleration at node $i$ in $j^{th}$ direction.

The residual for the $j^{th}$ rotational degree of freedom can be obtained as follows: $$R_j^0 = \int_{0}^{{}^0L} \rho I \left(\frac{x}{{}^0L} {\ddot{\theta}_j}^0 + \frac{{}^0L - x}{{}^0L} {\ddot{\theta}_j}^1 \right) \frac{x}{{}^0L} dx$$ $$R_j^1 = \int_{0}^{{}^0L} \rho I \left(\frac{x}{{}^0L} {\ddot{\theta}_j}^0 + \frac{{}^0L - x}{{}^0L} {\ddot{\theta}_j}^1 \right) \frac{{}^0L - x}{{}^0L} dx$$ where $I = I_y + I_z$ or user provided $I_x$ for $j = 1$, $I = I_z$ for $j = 2$ and $I = I_y$ for $j = 3$.

If $A_y$ and $A_z$ are non-zero, then there is coupling between $u_1$, $\theta_2$ and $\theta_3$, $u_2$ and $\theta_1$, and $u_3$ and $\theta_1$.

Point mass/inertia

NodalTranslationalInertia and NodalRotationalInertia are used to apply mass/inertia at a node. The mass and inertia terms contribute only to the residual of the node at which it is applied. Mass (m) behaves isotropically in all three coordinate directions resulting in the following $j^{th}$ translational nodal residual: $$R_j = m \ddot{u}_j$$

NodalRotationalInertia takes 6 rotational inertia components as input ($I_{xx}$, $I_{yy}$, $I_{zz}$, $I_{xy}$, $I_{xz}$, $I_{yz}$). If the coordinate system in which the moment of inertia components are provided is different from the global coordinate system, then x_orientation and y_orientation need to be provided as input so that a transformation matrix from the local coordinate system to the global coordinate system can be calculated.

Once the moment of inertia tensor in the global coordinate system is obtained, the residual for the rotational DOF in the $j^{th}$ direction is: $$R_j = \sum_{i=1}^3 I_{ji} \; \ddot{\theta}_i$$

Time-integration and damping

Rayleigh damping and Newmark and HHT time integration are calculated as in Damping but with StressDivergenceBeam and InertialForceBeam while using consistent mass/inertia matrix, and StressDivergenceBeam, NodalTranslationalInertia and NodalRotationalInertia while using point mass/inertia matrix.