Notes on beam models

TODO: explain variables

AeroFrame may be used to couple a so-called VLM model (vortex lattice method) and a beam FEM model as shown in [Dett19]. The VLM uses a surface-like mesh which is made up of quadrilaterals (panels), but the beam FEM mesh is line-like with discrete nodes. It is obvious that these discretisations do not coincide. A transfer and interpolation of loads and deformations becomes necessary. Since the load mapping from the VLM mesh onto the FEM beam is fairly illustrative, a few comments will be made here. Note that this page is based on or copied from [Dett19] where further information and context can be found.

Comments on beam load mapping schemes

When projecting aerodynamic loads onto the FEM beam, there are two potentially sensible mapping methods. The first method is to project loads onto the nearest node or beam element ([Bouc03], [Brau07], [RBBW10], [Seyw16]). The second method is to project loads chordwise (stripwise) onto the beam axis (e.g. [Seyw11], [BoEl16]). These two different mapping schemes are illustrated in Fig. 8. They are only geometrically equivalent in the case of unswept wings.

Fig. 8 Comparison of load mapping schemes for a swept wing (image from [Dett19])

In publications known to the author there is very little or no explanation of why either of these mapping schemes is chosen. Which of the projection methods is more sensible may depend on the real structure and the real load transfer. For instance, the orientation of stiff wing ribs can influence load paths inside the structure and therefore the deformation of the wing [RoWH14]. A simple beam model as indicated in Fig. 8 does not have any information about such structural details, and therefore no information of how loads are being transferred. As only the nodal loads matter, one has to make a choice of transferring loads onto the idealised elastic axis. In the following it will be shown that there is no equivalence between the two mentioned mapping methods. In general, work and displacements will differ.

Elastic work

Consider a straight cantilever beam in the global \(X\)-\(Y\) plane modelling the structure of a swept wing (Fig. 9). The beam axis \(x\) is oriented with an angle \(\phi\) with respect to the global \(Y\)-axis. Suppose there is a point force \(\mathbf{F}_\text{a} = (0, 0, F)^T\) acting at some off-axis point \(\mathbf{P}_\text{a}\). As discussed, the force may either be projected onto a point \(\mathbf{P}_1\) (chordwise projection) or onto a point \(\mathbf{P}_2\) (nearest element projection).

Fig. 9 Projection of off-axis loads onto a cantilever beam inclined with respect to the global coordinate system. The global coordinate system (blue axes) and the beam-local coordinate system (green axes) do not coincide. The off-axis loads can be projected onto the beam axis either using a parallel-to-X or a closest-element (closest-node) approach (image from [Dett19]).

For the following work consideration, it is convenient to directly express the projected loads in the beam-local coordinate system (axes denoted by lower-case \(x\), \(y\) and \(z\)). At point \(\mathbf{P}_1\) the equivalent, projected load is [1]

[1]	Notice that \(\phi < 0\).

(1)\[\begin{split}\mathbf{F}_{P_1, \text{loc}} = %% \mathbf{F}_\text{a} = %% \begin{pmatrix} 0 \\ 0 \\ F \\ \end{pmatrix} ~\text{and}\quad %% \mathbf{M}_{P_1, \text{loc}} = %% %%%%% \begin{pmatrix} - l_1 \cdot \sin \phi \\ - l_1 \cdot \cos \phi \\ 0 \\ \end{pmatrix} \times \begin{pmatrix} 0 \\ 0 \\ F \\ \end{pmatrix} = %% %%%%% F l_1 \begin{pmatrix} -\cos \phi \\ \sin \phi \\ 0 \\ \end{pmatrix}\end{split}\]

At point \(\mathbf{P}_2\) the equivalent load is

(2)\[\begin{split}\mathbf{F}_{P_2, \text{loc}} = %% \mathbf{F}_\text{a} = %% \begin{pmatrix} 0 \\ 0 \\ F \\ \end{pmatrix} ~\text{and} \quad \mathbf{M}_{P_2, \text{loc}} = %% %%%%% \begin{pmatrix} 0 \\ -l_2 \\ 0 \\ \end{pmatrix} \times \begin{pmatrix} 0 \\ 0 \\ F \\ \end{pmatrix} = %% %%%%% -F l_2 \begin{pmatrix} 1 \\ 0 \\ 0 \\ \end{pmatrix}\end{split}\]

Looking at the load components of the projected loads, it is apparent that bending is induced due to a transverse force component \(F_{z,\text{proj}} = F\) and due to a bending moment \(M_{y,\text{proj}}\) (\(M_{y,\text{proj}}=0\) when projecting onto \(\mathbf{P}_2\)), and twist is induced due to a torsional moment \(M_{x,\text{proj}}\). The internal elastic energy due to bending about \(y\) and torsion (equal to the work done by the external loads) is given as [Sund10]

(3)\[W = \displaystyle\int_0^{L_2} \left[ \frac{1}{2} \cdot E \cdot I_y \left(\frac{\text{d}{}^2 u_z}{\text{d}{x}^2}\right)^2 + \frac{1}{2} \cdot G \cdot J \left(\frac{\text{d}{\Theta_x}}{\text{d}{x}}\right)^2 \right] \text{d}{x}\]

where the notation from TODO is used, and \(L_2\) is the beam length. For the sake of keeping the following relations somewhat simpler, the bending stiffness \(E \cdot I_y\) and the torsional stiffness \(G \cdot J\) are assumed to be constant. It can be shown that the projected point loads acting at \(x = x_\text{i}\) will introduce work given by the following expression (c.f. [Sund10], [Megs16]).

(4)\[\begin{split}W &= %%% \displaystyle\int_0^{L_2} \left( \frac{\widetilde{M}^2_{y}}{2 \cdot E \cdot I_y} + \frac{\widetilde{M}^2_x}{2 \cdot G \cdot J} \right) \text{d}{x} = %%% \dots \nonumber \\ %%% &= \frac{{M}_{y,\text{proj}}^2 \cdot x_\text{i}}{2 \cdot E \cdot I_y} - \frac{{M}_{y,\text{proj}} \cdot {F}_{z,\text{proj}} \cdot x_\text{i}^2}{2 \cdot E \cdot I_y} + \frac{{F}_{z,\text{proj}}^2 \cdot x_\text{i}^3}{6 \cdot E \cdot I_y} + \frac{{M}_{x,\text{proj}}^2 \cdot x_\text{i}}{2 \cdot G \cdot J}\end{split}\]

\(\widetilde{M}_y\) is the internal bending moment about the \(y\)-axis and \(\widetilde{M}_x\) is the internal torsional moment. \(F_{z,\text{proj}}\), \(M_{x,\text{proj}}\) and \(M_{y,\text{proj}}\) are the projected loads (external loads) introduced at a beam position \(x = x_\text{i}\). After inserting the loads from eqs. (1) and (2) into (4) and some algebraic manipulations, the total work done by the external loads for the two projection cases can be expressed. In the first case, when projecting onto \(\mathbf{P}_1\), the total work is

(5)\[W_{P_1} = %% %% \underbrace{\frac{F^2 \cdot L_1 \cdot l_1^2 \cdot \sin^2 \phi}{2 \cdot E \cdot I_y}}_{\text{a}} %% + \underbrace{\frac{-F^2 \cdot L_1^2 \cdot l_1 \cdot \sin \phi}{2 \cdot E \cdot I_y}}_{\text{b}} %% + \underbrace{\frac{F^2 \cdot L_1^3}{6 \cdot E \cdot I_y}}_{\text{c}} %% + \underbrace{\frac{F^2 \cdot L_1 \cdot l_1^2 \cdot \cos^2 \phi}{2 \cdot G \cdot J}}_{\text{d}} %%\]

In the second case, when projecting onto \(\mathbf{P}_2\), the total work is

(6)\[\begin{split}\begin{align} W_{P_2} &= %% \frac{F^2 \cdot L_2^3}{6 \cdot E \cdot I_y} + \frac{F^2 \cdot L_2 \cdot l_2^2}{2 \cdot G \cdot J} \\ &= \left\{ \text{using} \quad L_2 = L_1 - l_1 \cdot \sin \phi \quad \text{and} \quad l_2 = l_1 \cdot \cos \phi \right\} = \dots \\[3mm] &= \underbrace{\frac{F^2 \cdot L_1 \cdot l_1^2 \cdot \sin^2 \phi}{2 \cdot E \cdot I_y}}_{\text{a}} %% + \underbrace{\frac{-F^2 \cdot L_1^2 \cdot l_1 \cdot \sin \phi}{2 \cdot E \cdot I_y}}_{\text{b}} %% + \underbrace{\frac{F^2 \cdot L_1^3}{6 \cdot E \cdot I_y}}_{\text{c}} %% + \underbrace{\frac{F^2 \cdot L_1 \cdot l_1^2 \cdot \cos^2 \phi}{2 \cdot G \cdot J}}_{\text{d}} %% \\ % Additional terms &+ \underbrace{\frac{-F^2 \cdot l_1^3 \cdot \sin^3 \phi}{6 \cdot E \cdot I_y}}_{\text{e}} %% + \underbrace{\frac{-F^2 \cdot l_1^3 \cdot \sin \phi \cdot \cos^2 \phi}{2 \cdot G \cdot J}}_{\text{f}} %% \end{align}\end{split}\]

In (6) simple geometric relations between lengths \(L_1\), \(L_2\), \(l_1\) and \(l_2\) have been utilised. Comparing eqs. (5) and (6) shows that the first four terms (a, b, c and d) occur in both equations (with the assumption of constant \(E \cdot I_y\) and \(G \cdot J\)). However, in the mapping case \(\mathbf{P}_2\), there are two additional terms (e and f). Term e results from the bending deformation and term f from the torsional deformation. In the case illustrated in Fig. 9 (\(\phi < 0\)), these terms will be positive, hence the beam will store more elastic energy when transferring loads onto point \(\mathbf{P}_2\) instead of \(\mathbf{P}_1\).

Suppose that the off-axis force was located on the opposite side of the beam axis (force shifted in negative \(X\)-direction, like \(\mathbf{F}_\text{b}\)) at the same negative beam inclination \(\phi\). With the closest-node mapping, the force would be projected closer to the wing root than with the projection parallel to \(X\). With analogous reasoning, it can be shown that the beam takes up less energy when the load is projected closer to the root.

Example

The FEM tool FramAT allows to assess work and displacements in a more convenient and general way. The API for off-axis loads allows to project onto the nearest node or parallel to \(X\) (node with the closest \(Y\)-coordinate). In the FEM formulation, the elastic energy is computed as (c.f. [CMPW02])

(7)\[W = \frac{1}{2} \cdot \mathbf{U}^T \cdot \mathbf{K} \cdot \mathbf{U}\]

where \(\mathbf{U}\) is the vector of nodal deformations and \(\mathbf{K}\) the global stiffness matrix.

A straight cantilever beam with length \(L_2 = 1.5 \text{m}\) and variable inclination \(\phi\) (as illustrated in Fig. 9) was analysed. The bending stiffness \(E \cdot I_y\) and the torsional stiffness \(G \cdot J\) were both set to \(1 \text{N/m}^2\). The beam was loaded with off-axis forces \(\mathbf{F}_\text{a} = (0 , 0 , F_{z,\text{a}})^T\) at point \(\mathbf{P}_\text{a}\) and \(\mathbf{F}_\text{b} = (0 , 0 , F_{z,\text{b}})^T\) at \(\mathbf{P}_\text{b}\). The points of attack were computed as \(\mathbf{P}_\text{a} = \mathbf{P}_1 + l_1 \cdot \mathbf{e}_x\) and \(\mathbf{P}_\text{b} = \mathbf{P}_1 - l_1 \cdot \mathbf{e}_x\) where \(\mathbf{P}_1 = L_1 \cdot (-\sin \phi , \cos \phi , 0)^T\). The length \(L_1\) was set to \(1 \text{m}\) and \(l_1\) to \(0.5 \text{m}\). Two different load cases were analysed. In the first load case (A), \(F_{z,\text{a}}\) was set to \(1 \text{N}\) and \(F_{z,\text{b}} = 0\). In the second load case (B), \(F_{z,\text{a}}\) and \(F_{z,\text{b}}\) were both set to \(1 \text{N}\). All analyses used 200 beam elements (1206 d.o.f.). Fig. TODO1 and TODO2 show the elastic energy (or work) and the beam tip deflection \(u_z\) for different beam inclinations \(\phi\).

TODO add figure 1

TODO add figure 2

In case B, a load pair is applied. With the parallel-to-\(X\) projection the two point forces only introduce a force \(F_{z,\text{proj}}\) into the beam since torsion and bending moments (\(M_{x,\text{proj}}\) and \(M_{y,\text{proj}}\)) are cancelled out by contributions acting in opposite directions. In other words, the beam is only loaded in bending due to a point force which is of equal magnitude for all beam inclinations \(\phi\). Hence, deflection and work are constant. When projecting with the closest-node approach, there is always one projection point which is closer and one projection point which is further away from point \(\mathbf{P}_1\) (with the exception \(\phi=0\)). Both work and tip deflection appear to always be larger than for the parallel-to-\(X\) mapping scheme. The differences grow larger for larger angles \(|\phi|\). With analytical methods it can be easily shown that two individual forces \(F_z\) applied at \(x = x_\text{i} + \Delta x\) and at \(x = x_\text{i} - \Delta x\) will cause a cantilever beam to deflect more than if \(2 \cdot F_z\) were applied at \(x_i\) (given constant bending stiffness).

The presented example illustrates that work and deflections generally differ depending on the load mapping choice. Additional differences may arise when the beam has a variable stiffness. Also the overall wing span and the exact locations of the points of attack of the aerodynamic loads will affect results. In the shown example, differences become larger for larger beam inclinations, as the distance between the projection points \(\mathbf{P}_1\) and \(\mathbf{P}_2\) grows larger. A general recommendation cannot be made. The preferred (more accurate) mapping scheme may in the end depend on the real internal load transfer which is affected by components such as stiff ribs which are not modelled with a single beam. Nevertheless, one should be aware of differences that can arise in these modelling decisions and their general implications.

Note

This summary is based on/copied from [Dett19] with the authors permission.