Nonlinear finite element analysis of lattice core sandwich beams

TL;DR

This paper develops a nonlinear finite element model for analyzing the bending behavior of micropolar Timoshenko beams, specifically applied to various lattice core sandwich beams, capturing moderate rotations and complex core topologies.

Contribution

It introduces a nonlinear 1-D finite element model for micropolar beams with a novel application to different lattice core sandwich structures.

Findings

Good agreement with 2-D finite element results

Effective reduction of shear and membrane locking

Applicable to multiple lattice core topologies

Abstract

A geometrically nonlinear finite element model is developed for the bending analysis of micropolar Timoshenko beams using the principle of virtual displacements and linear Lagrange interpolation functions. The nonlinearity enters the model via a nonlinear von K\'arm\'an strain term that allows the micropolar beam to undergo moderate rotations. The nonlinear micropolar Timoshenko beam is used as an equivalent single layer model to study four different lattice core sandwich beams. A two-scale energy method is used to derive the micropolar constitutive equations for web, hexagonal, Y-frame and corrugated core topologies. Various bending cases are studied numerically using the developed 1-D finite element model. Reduced integration techniques are used to overcome the shear and membrane locking. The present 1-D results are in good agreement with the corresponding 2-D finite element beam frame results for global bending.

Equations

$$U_x(x,y)=u_x(x)+y\varphi (x),U_y(x,y)=u_y(x),\Psi (x,y)=\psi (x),$$

$${\epsilon }_x$$

$${\epsilon }_{xy}$$

$${\gamma }_s$$

$${\gamma }_a$$

$$U_x(x\pm l/2,\pm h/2)$$

$$U_y(x\pm l/2,\pm h/2)$$

$$\Psi (x\pm l/2,\pm h/2)$$

$${\mathbf{d}}^c={\mathbf{T}}_u^c\mathbf{u}+{\mathbf{T}}_{\epsilon }^c\mathbf{\epsilon },$$

$${\mathbf{d}}^c={\left\{U_{x,1}{U}_{y,1}{\Psi }_1{U}_{x,2}{U}_{y,2}{\Psi }_2{U}_{x,3}{U}_{y,3}{\Psi }_3{U}_{x,4}{U}_{y,4}{\Psi }_4\right\}}^{\text{T}}$$

$$\mathbf{u}$$

$$\mathbf{\epsilon }$$

$${\mathbf{d}}^m={\mathbf{T}}_u^m\mathbf{u}+{\mathbf{T}}_{\epsilon }^m\mathbf{\epsilon }.$$

$$\mathbf{d}=\left\{\begin{array}{c}{\mathbf{d}}^c\\ {\mathbf{d}}^m\end{array}\right\}=\left[\begin{array}{cc}{\mathbf{T}}_u^c& {\mathbf{T}}_{\epsilon }^c\\ {\mathbf{T}}_u^m& {\mathbf{T}}_{\epsilon }^m\end{array}\right]\left\{\begin{array}{c}\mathbf{u}\\ \mathbf{\epsilon }\end{array}\right\}.$$

$$W=\frac{1}{2}{\mathbf{d}}^{\text{T}}\mathbf{kd},$$

$$W=\frac{1}{2}{\mathbf{\epsilon }}^{\text{T}}{\mathbf{T}}_{\epsilon }^{\text{T}}\mathbf{k}{\mathbf{T}}_{\epsilon }^{}\mathbf{\epsilon }.$$

$$W_0^l\equiv \frac{W}{l}=\frac{1}{2}{\mathbf{\epsilon }}^{\text{T}}\mathbf{C}\mathbf{\epsilon }$$

$$\mathbf{C}=\frac{1}{l}{\mathbf{T}}_{\epsilon }^{\text{T}}\mathbf{k}{\mathbf{T}}_{\epsilon }^{}.$$

$$\hat{\mathbf{\epsilon }}=\{{\hat{\epsilon }}_x^0{\kappa }_x{\gamma }_s{\gamma }_a{\kappa }_{xz}{\}}^{\text{T}}.$$

$${\hat{\epsilon }}_x^0=u_x^{\prime }+\frac{1}{2}{\left(u_y^{\prime }\right)}^2.$$

$$\mathbf{S}\equiv \frac{\partial {\hat{W}}_0^l}{\partial \hat{\mathbf{\epsilon }}}=\frac{1}{2}\frac{\partial }{\partial \hat{\mathbf{\epsilon }}}\left({\hat{\mathbf{\epsilon }}}^{\text{T}}\mathbf{C}\hat{\mathbf{\epsilon }}\right)=\mathbf{C}\hat{\mathbf{\epsilon }},$$

$$\left\{\begin{array}{c}{N}_x\\ M_x\\ Q_s\\ Q_a\\ P_{xz}\end{array}\right\}=\left[\begin{array}{ccccc}{C}_{11}& C_{12}& 0& 0& C_{15}\\ C_{12}& C_{22}& 0& 0& C_{25}\\ 0& 0& C_{33}& C_{34}& 0\\ 0& 0& C_{34}& C_{44}& 0\\ C_{15}& C_{25}& 0& 0& C_{55}\end{array}\right]\left\{\begin{array}{c}{\hat{\epsilon }}_x^0\\ {\kappa }_x\\ {\gamma }_s\\ {\gamma }_a\\ {\kappa }_{xz}\end{array}\right\},$$

$$Q_s=\frac{Q_{xy}+Q_{yx}}{2}\text{and}{Q}_a=\frac{Q_{xy}-Q_{yx}}{2},$$

$$\left\{\begin{array}{c}{N}_x\\ M_x\\ Q_s\\ Q_a\\ P_{xz}\end{array}\right\}=\left[\begin{array}{ccccc}2E{A}_f& 0& 0& 0& 0\\ & \frac{E{A}_f{h}^2}{2}+\Theta & 0& 0& \Theta \\ & & \frac{6E{I}_f+\Theta }{l^2}& \frac{6E{I}_f-\Theta }{l^2}& 0\\ & \text{SYM}& & \frac{6E{I}_f+\Theta }{l^2}& 0\\ & & & & 2E{I}_f+\Theta \end{array}\right]\left\{\begin{array}{c}{\hat{\epsilon }}_x^0\\ {\kappa }_x\\ {\gamma }_s\\ {\gamma }_a\\ {\kappa }_{xz}\end{array}\right\}$$

$$\Theta =\frac{3E{I}_w{k}_{\theta }l}{6E{I}_w+k_{\theta }h}.$$

$${\mathbf{C}}_{\text{hex}}=\left[\begin{array}{ccccc}4.26438\cdot 10^7& 0& 0& 0& 0\\ 0& 340740& 0& 0& 27.1761\\ 0& 0& 94735.4& 8319.08& 0\\ 0& 0& 8319.08& 3315.25& 0\\ 0& 27.1761& 0& 0& 44.3467\end{array}\right]$$

$${\mathbf{C}}_{\text{Y}}=\left[\begin{array}{ccccc}1.26053\cdot 10^8& 11696.2& 0& 0& 5292.50\\ 11696.2& 6.10097\cdot 10^6& 0& 0& 1164.35\\ 0& 0& 42094.9& 9541.42& 0\\ 0& 0& 9541.42& 5302.30& 0\\ 5292.50& 1164.35& 0& 0& 1012.92\end{array}\right]$$

$${\mathbf{C}}_{\text{corr}}=\left[\begin{array}{ccccc}1.26018\cdot 10^8& -3902.45& 0& 0& -3902.45\\ -3902.45& 6.09926\cdot 10^6& 0& 0& 858.539\\ 0& 0& 2.09792\cdot 10^7& 3734.96& 0\\ 0& 0& 3734.96& 5078.42& 0\\ -3902.45& 858.539& 0& 0& 1334.05\end{array}\right],$$

$$\delta W=\delta W_I-\delta W_E=0,$$

$$\delta W_I$$

$$\displaystyle=\int_{x_{a}}^{x_{b}}\big{[}N_{x}(\delta u_{x}^{\prime}+u_{y}^{\prime}\delta u_{y}^{\prime})+M_{x}\delta\phi^{\prime}+P_{xz}\delta\psi^{\prime}$$

$$\displaystyle\qquad\quad+Q_{xy}(\delta u_{y}^{\prime}-\delta\psi)+Q_{yx}(\delta\phi+\delta\psi)\big{]}dx,$$

$$\delta W_E={\int }_{x_a}^{x_b}q\delta u_{y}dx+\sum _{i=1}^8{Q}_i\delta U_i,$$

$$U_1$$

Full text

††journal: Elsevier

1 Introduction

A typical sandwich panel consists of a thick, low-stiffness core between two relatively thin but stiff face sheets. The face sheets take bending and in-plane loads while the core carries transverse shear loads (2; 44). The face sheets and core can be made of the same or different materials and some possible core structures include, for example, foam, solid, honeycomb, and truss cores (45). A number of manufacturing techniques are available for sandwich panels (22; 46). A sandwich panel generally has a high bending stiffness compared to a single solid plate of the same dimensions made of either the face sheet or core material and the panel weighs considerably less than the solid plate making it a weight-efficient structure (44). Most early sandwich panels had isotropic face sheets but with the development of fiber reinforced composites, construction of sandwich panels with composite face sheets has become possible. Because of their high stiffness-to-weight ratios, structural efficiency and durability, sandwich panels are suitable candidates for high-speed aircraft and space applications (41). They have also found applications, for example, in shipbuilding (37; 26) and other marine applications (9; 29; 25) and in the construction of bridges and buildings (15; 10; 11; 30; 12).

The number of applications for sandwich panels is increasing rapidly. The required accuracy in the structural analysis of the panels depends on the type of the application considered. For example, in air-crafts a very detailed response of the sandwich structure may be required, whereas an overall global response may suffice in residential buildings when the natural vibration frequencies are of interest, for example. In any case, there is a need for appropriate modeling tools for different applications. Reviews on the modeling of sandwich structures have been given by several authors (31; 20; 13; 40; 8). Modeling methods for sandwich panels can be broadly classified as: (a) Complete 3-D analysis (computational or analytical), with complete details of the face sheets and the core structure considered; (b) layer-wise modeling with the faces and core considered as separate continuum layers (33); (c) statically equivalent single layer (ESL) models. Although computational 3-D and layer-wise analyses give very detailed stress distributions for the panels, they come with the inherent disadvantage of including a large number of variables and, thus, the computational analysis of them can be very burdensome. Therefore, equivalent single layer theories such as the ESL first-order shear deformation (FSDT) beam and plate models are attractive especially when the global response of the structure is of main interest without accounting for every small detail. Extensive literature exists on the modeling of sandwich beam, plates and shells by ESL theories, see, for example (42; 7; 6; 20; 1).

It was shown recently that all-steel web-core sandwich beams deform so that when an ESL-FSDT model is used, the model has to take into account anti-symmetric shear deformations in order for the response of the sandwich structure to be captured accurately in some applications (23). In more detail, the constituents of a 2-D or 3-D web-core sandwich beam model do not exhibit any anti-symmetric shear strains, but when the problem is reduced to a 1-D ESL beam problem essentially by thickness integration, the anti-symmetric behavior needs to be considered via a 1-D micropolar Timoshenko beam model. In this paper, we develop a geometrically nonlinear finite element model based on the micropolar Timoshenko beam theory.

With the revived interest in micropolar elasticity (17), considerable work has been put into developing appropriate finite element models for micropolar continua in general; see, for example, (32; 27; 38; 47). To list a few recent finite element models for micropolar plates we mention the works of (4; 3) and (18). Various finite element models have been proposed for the bending analysis of micropolar beams as well. (21) used 3-D non-compatible finite elements to analyze the bending of beams, and (27) proposed three different elements for plane micropolar elasticity and used them to analyze thin in-plane beams. (19) developed a 1-D micropolar beam finite element model using Lagrange interpolation functions. (36) derived a finite element model for a micropolar Timoshenko beam with the microrotation assumed to be equal to the cross-sectional rotation. More recently, (23) proposed nodally-exact 1-D finite element to analyze micropolar Timoshenko beams and (3) proposed a 27-node 3-D finite element for the analysis of beams. Only linear strains were considered in developing the finite element models in all the above mentioned papers.

In this study, in order to develop the geometrically nonlinear micropolar 1-D beam finite element, we start from the principle of virtual displacements and use the linear Lagrange interpolation functions. The originally linear micropolar Timoshenko beam model of (23) is enriched by nonlinear von Kármán strains in order to account for the moderate rotations of the beam (16). Lattice cores, namely, hexagonal, corrugated and Y-frame cores in addition to the web-core topology studied earlier are considered. To this end, the two-scale constitutive modeling method presented by (24) is first generalized to cores other than the web-core.

In more detail, the rest of the paper is organized as follows. In Section 2, a brief review on the linear 1-D micropolar Timoshenko beam model is given followed by the two-scale constitutive modeling of the four different lattice cores after which the geometrically nonlinear micropolar Timoshenko beam equations are derived. In Section 3, the displacement based finite element model is formulated using the principle of virtual displacements. The element stiffness matrices are derived and the used nonlinear iterative procedures and techniques to avoid numerical locking are discussed. In Section 4, numerical examples are solved using the developed 1-D micropolar beam finite element model and are compared with the corresponding results from 2-D finite element analyses. Finally, concluding remarks are given in Section 5.

2 Geometrically nonlinear micropolar Timoshenko beam

Two-scale constitutive modeling of lattice materials in the context of linear micropolar Timoshenko beam theory is first carried out. The strain energy density for a geometrically nonlinear micropolar Timoshenko beam is obtained by retaining the linear constitutive matrix while the strain vector is augmented with von Kármán nonlinearity. The hyperelastic constitutive relations, i.e., the stress resultant equations for the nonlinear beam are derived from the strain energy for four lattice materials. The equilibrium equations for the nonlinear beam in terms of the stress resultants, and the corresponding weak form for finite element developments are attained by employing the principle of virtual displacements.

2.1 Displacements and linear strains

The two-dimensional displacements $U_{x}$ and $U_{y}$ and the independent microrotation $\Psi$ of a micropolar Timoshenko beam can be written in terms of central axis kinematic variables as (24)

[TABLE]

where $u_{x}$ is the axial displacement, $\phi$ is the rotation of the cross-section, $u_{y}$ is the transverse deflection, and $\psi$ is the microrotation. The nonzero infinitesimal strains of the beam are

[TABLE]

where the prime “′” on the variables denotes differentiation with respect to $x$. The symmetric and anti-symmetric shear strains of the beam are defined as

[TABLE]

respectively. The anti-symmetric part is defined by the difference between the macrorotation $\omega_{z}$ and the microrotation $\psi$. For $\omega_{z}=\psi$, the relative strains reduce to their classical definitions (5), for example, $\varepsilon_{xy}=u_{y}^{\prime}-\omega_{z}=(u_{y}^{\prime}+\phi)/2$.

2.2 Two-scale constitutive modeling

Figure 1 shows a rectangular unit cell attached to an arbitrary cross section of a micropolar Timoshenko beam. The unit cell of length $l$ and height $h$ represent the periodic microstructure of the macrostructural beam ($l\leq\textrm{beam length}$). The unit cell corner displacements in Fig. 1 are expressed in terms of the micropolar cross-sectional displacements $U_{x}$ and $U_{y}$ and rotation $\Psi$. With distance from an arbitrary beam cross section, Taylor series expansions of Eqs. (1) lead to

[TABLE]

where the micropolar strains (2)–(4) have been imposed on the cross-sectional rotation $\phi$ and the displacement gradients. Using the node numbering of Fig. 1, we can write the discrete-to-continuum transformation given by Eqs. (5)–(7) for the corner nodes in matrix form

[TABLE]

where the generalized discrete displacement vector is

[TABLE]

and the vectors for the continuous variables read

[TABLE]

The transformation matrices $\mathbf{T}^{c}_{u}$ and $\mathbf{T}^{c}_{\varepsilon}$ are given in Appendix A. For additional mid-nodes located at $\pm l/2$ and connected to the neighboring unit cells, a similar transformation can be written as

[TABLE]

The complete discrete-to-continuum transformation reads

[TABLE]

The transformation by Eq. (13) for nodes located at the unit cell edges $x=\pm l/2$ may be applied to different finite element based lattice unit cells once static condensation has been applied at all nodes located between $l/2<x<l/2$. The unit cell can be modeled, for example, by using Euler-Bernoulli or Timoshenko beam finite elements as both choices result in a system that is consistent with the generalized displacement vector (13). As for their material composition, the elements can be heterogeneous and anisotropic. Several lattice cores will be considered in Section 2.3.

The strain energy of a unit cell modeled by beam elements can be written as

[TABLE]

where $\mathbf{k}$ is the global stiffness matrix of the unit cell that corresponds to the master degrees of freedom after the static condensation. It is straightforward to verify that the displacement terms (10) do not contribute to the strain energy of any of the lattice core unit cells considered in this paper and we obtain

[TABLE]

We define the 1-D density of the unit cell strain energy as

[TABLE]

where the constitutive matrix is given by

[TABLE]

For geometrically nonlinear unit cell analysis in which strains are small but rotations may be moderate, we assume that the constitutive matrix (17) remains the same but the axial normal strain is modified so that the strain vector becomes

[TABLE]

in which the nonlinear von Kármán term is included in

[TABLE]

The unit cell represents a lattice material of which the micropolar beam is made of. We write the hyperelastic constitutive relations for the geometrically nonlinear micropolar beam continuum as

[TABLE]

where $\mathbf{S}$ is now the stress resultant vector of the 1-D micropolar beam and for which the general explicit form considered in this paper is

[TABLE]

where $N_{x}$ is the axial force, $M_{x}$ and $P_{xz}$ are the bending and couple-stress moments, respectively, and $Q_{s}$ and $Q_{a}$ are the symmetric and anti-symmetric shear forces defined as (24)

[TABLE]

respectively. Note that Eq. (20) implies that the bridging of the two scales, i.e., the macroscale (beam) and the microscale (unit cell), is founded on an assumption of strain energy equivalence between them.

2.3 Modeling of lattice materials

In this section, we consider four different lattice material sandwich beam cores, namely, the web-core studied earlier in the linear context (24) and the hexagonal, Y-frame and corrugated cores presented in Fig. 2. All cores are made of steel and are modeled using linearly elastic isotropic nodally-exact Euler-Bernoulli beam finite elements. All unit cell beam constituents have a rectangular cross-section and in this paper the width of all unit cells is 0.05 m.

The web-core is modeled using four Euler-Bernoulli beam elements that can represented by the dashed lines in Fig. 1. In this case, only the corner node transformation (8) needs to be considered and the application of this in Eq. (14) results in the constitutive equations

[TABLE]

where

[TABLE]

In the constitutive matrix, $EA_{f}$ and $EI_{f}$ are the axial and bending stiffnesses of the horizontal faces, respectively. For the vertical webs we have $EA_{w}$, $EI_{w}$ and $k_{\theta}$ for the axial, bending and rotational joint stiffnesses, respectively. The webs are modeled using special-purpose Euler-Bernoulli elements with rotational springs at both ends to account for the flexibility of the laser-welded joints (28; 14; 39). The numerical values of the core parameters are $E_{f}=212$ GPa, $E_{w}=200$ GPa, $k_{\theta}=2675$ Nm and $\nu=0.3$ for the face and web Young’s moduli, rotational joint stiffness and Poisson ratio, respectively. The face and web thicknesses are $t_{f}=3$ mm and $t_{w}=4$ mm, respectively. The height, i.e, the distance between the face central axes is $h=43$ mm.

For the other cores displayed in Fig. 2, Young’s modulus and Poisson ratio are $E=210$ GPa and $\nu=0.3$, respectively. Other relevant parameters are given in Fig. 2. The hexagonal core includes two mid-nodes in addition to the four corner nodes that need to be taken into account essentially to ensure connectivity between neighboring unit cells on the micropolar continuum level. Static condensation is applied at the inner nodes of the unit cell. It is difficult to obtain a meaningful symbolic form for the hexagonal constitutive matrix, in numerical form we have

[TABLE]

The constitutive matrix (25) is of the same form as that of the web-core with the exception that for the hexagonal core $C_{33}\neq C_{44}$. For the Y-frame and corrugated cores we obtain

[TABLE]

and

[TABLE]

respectively. The axial and classical sandwich bending stiffnesses in Eqs. (25)–(27) are practically given by $C_{11}\approx 2EA_{f}$ and $C_{22}\approx EA_{f}h^{2}/2$, respectively. We see that due to the lack of symmetry about the $x$-axis, the coupling terms $C_{12}$ and $C_{15}$ appear in the constitutive matrices of the Y-frame and corrugated cores. In addition, the symmetric shear stiffness $C_{33}$ of the corrugated core is very high in comparison to that of the other cores because the corrugated lattice core has a stretch-dominated shear-carrying mechanism while the others cores are bending-dominated. In other words, when bent, the constituents of the corrugated lattice core act as axial rods without significant bending so that the core is very stiff. The differences between bending- and stretch-dominated cores in geometrically nonlinear bending problems will be studied further by numerical examples in Section 4. Finally, it is easy verify that all eigenvalues of each constitutive matrix above are positive which means that the matrices are positive definite. It follows that each lattice core material is stable in the conventional sense (i.e., strain energy is positive for nonzero strains).

2.4 Geometrically nonlinear beam equations

The principle of virtual displacements can be stated for a geometrically nonlinear micropolar Timoshenko beam as

[TABLE]

The virtual strain energy stored in a typical beam element, $\Omega=(x_{a},x_{b})$, and consistent with the strains used in the previous section, is

[TABLE]

where $\hat{\varepsilon}_{x}=\hat{\varepsilon}_{x}^{0}+y\kappa_{x}$ [cf. Eqs. (2) and (19)]. The external virtual work is

[TABLE]

where $q(x)$ is the distributed transverse load, $Q_{i}$ are the generalized external forces at the beam ends and $U_{i}$ are the associated generalized displacements defined as

[TABLE]

Equation (28) conforms to the micropolar beam finite element presented in Fig. 3 and provides the weak form for the finite element formulation in Section 3. On the other hand, Eq. (28) yields the Euler-Lagrange equations

[TABLE]

where the effective transverse shear force is

[TABLE]

Finally, when the corresponding generalized displacements are not defined, the natural (or force) boundary conditions become

[TABLE]

The equilibrium equations (32) may be presented in terms of displacements by using the constitutive relations (21). The governing equations could also have been derived by employing the principle of minimum total potential energy and the strain energy density $\hat{W}_{0}^{l}$ given in Eq. (20).

3 Geometrically nonlinear micropolar Timoshenko beam element

3.1 Finite element formulation

The finite element model is developed using the statement of the principle of virtual displacements (28) over a typical element $\Omega^{e}=(x_{a},x_{b})$, which takes the explicit form

[TABLE]

After using the constitutive equations in (21), the following four statements for a typical finite element $\Omega^{e}=(x_{a},x_{b})$ which are equivalent to Eq. (35) can be written

[TABLE]

[TABLE]

The primary variables are $u_{x},u_{y},\phi$ and $\psi$. These variables are approximated using Lagrange interpolation functions $L_{j}^{(J)}$, $(J=1,2,3,4)$ (see, for example, (35)) so that

[TABLE]

By substituting Eq. (40) for $u_{x},u_{y},\phi$ and $\psi$ and putting $\delta u_{x}=L_{j}^{(1)},\delta u_{y}=L_{j}^{(2)},\delta\phi=L_{j}^{(3)}$ and $\delta\psi=L_{j}^{(4)}$ into the weak-form statements in Eqs. (36)–(39), the finite element equations for a typical beam element can be expressed as

[TABLE]

The stiffness coefficients $K_{ij}^{\alpha\beta}$ and force coefficients $F_{i}^{\alpha}$ ($\alpha,\beta=1,2,3,4$ and $i,j=1,2$) are defined as

[TABLE]

[TABLE]

[TABLE]

3.2 Solution of nonlinear equations

The nonlinear finite element equations (41) are solved iteratively using the Newton’s iteration procedure (see (34)). The linearized element equation at the beginning of the $r^{th}$ iteration will take the form:

[TABLE]

Where $\mathbf{U}^{(e)(r-1)}$ is the generalized nodal displacement vector of element $e$ at the end of $(r-1)^{th}$ iteration and $\Delta\mathbf{U}^{(e)(r)}$ is the incremental displacement vector of element $e$ at the $r^{th}$ iteration defined as

[TABLE]

The residual vector $\mathbf{R}^{(e)}(\mathbf{U}^{(e)(r-1)})$, computed after the end of $(r-1)^{th}$ iteration, is defined as

[TABLE]

Once the residual vector is computed using Eq. (46), the tangent stiffness matrix $\mathbf{T}^{(e)}$ can be calculated using the following definition

[TABLE]

For a typical element, writing Eq. (44) in a fashion similar to Eq. (41), we get

[TABLE]

where the notation $\mathbf{U}^{1}=\mathbf{u_{x}}$, $\mathbf{U}^{2}=\mathbf{u_{y}}$, $\mathbf{U}^{3}=\boldsymbol{\phi}$ and $\mathbf{U}^{4}=\boldsymbol{\psi}$ is used. Then the coefficients of the tangent stiffness matrix, $T_{ij}^{\alpha\beta}$ ($\alpha,\beta=1,2,3,4$ and $i,j=1,2$), in the above equation can be computed as

[TABLE]

In explicit terms, the coefficients read

[TABLE]

[TABLE]

The element equations computed from Eq. (48) are assembled and solved (after the imposition of the boundary conditions) to obtain the global incremental displacement vector $\Delta\mathbf{U}$ at the $r^{th}$ iteration. The normalized difference between solution vectors from two consecutive iterations, measured with Euclidean norm, is computed at the end of each iteration. If the value computed is less than a preselected tolerance $\varepsilon$ further iterations are terminated

[TABLE]

Further, acceleration of convergence for some type of nonlinearities may be achieved by using weighted-average of solutions from the last two iterations rather than the solution from the last iteration to evaluate the coefficient matrix:

[TABLE]

where $\beta$ is known as the acceleration parameter. The value of $\beta$ depends on the nature of nonlinearity and the type of problem considered. Once the convergence is obtained the final solution is computed using

[TABLE]

3.3 Shear and Membrane locking

The finite element model used in this paper uses linear interpolation on both $u_{y}$ and $\phi$. However, in the thin beam limit, when linear interpolation is used for $u_{y}$, the cross-sectional rotation $\phi$ should approach $-\left(du_{y}/dx\right)$, which is necessarily constant. But since $\phi$ is also interpolated as linear, it can never be constant. This inconsistency causes what is known as shear locking (see (34)). To avoid this inconsistency, we may use equal interpolation on both $u_{y}$ and $\phi$ but treat $\phi$ as constant while evaluating the symmetric $\gamma_{s}$ and anti-symmetric $\gamma_{a}$ shear strains. This amounts to using reduced Gauss quadrature rule in evaluating the integrals containing constants $C_{33}$, $C_{34}$ and $C_{44}$ while computing the element coefficient matrices of Eq. (41) and Eq. (48).

When von Kármán nonlinearity is included, there is coupling between $u_{x}$ and $u_{y}$ which causes the beam to undergo axial displacement even when there are no axial forces. But in the case of hinged-hinged beam, there are no constraints on $u_{x}$ at the boundaries, thus causing the beam to roll over freely without axial strain, i.e,

[TABLE]

In order to satisfy this we need

[TABLE]

In essence, we need to have the same degree of polynomial variation on both $\left(du_{x}/dx\right)$ and $\left(du_{y}/dx\right)^{2}$. But when equal interpolation of degree greater than one is used for both $u_{x}$ and $u_{y}$ this criteria cannot be satisfied and leads to what is known as membrane locking (see (34)). To overcome this we have to treat $\left(du_{y}/dx\right)^{2}$ as same order as $\left(du_{x}/dx\right)$. This is achieved using reduced integration while evaluating all the nonlinear terms of the element coefficients matrices of Eqs. (41) and (48).

4 Numerical bending examples

4.1 General Setup

The developed micropolar beam finite element model is used for bending analysis of lattice core sandwich beams. The four structural cores considered in Section 2.3 are used in the calculations. Both geometrically linear and nonlinear cases are analyzed using the 1-D beam model. 2-D reference solutions are computed using Euler–Bernoulli FE beam frames modeled by Abaqus; the pins in simply-supported cases are at the central axis of the 2-D frame so that the model corresponds to 1-D cases.

4.2 Bending of a web-core beam

A beam consisting of 24 web-core unit cells is considered first. The length of each web-core unit cell is $l=0.12$ m resulting in a total beam length of $L=2.88$ m. The beam is analyzed for two different boundary conditions, namely, a fixed-fixed case and a three-point-bending setup. For the fixed-fixed case the boundaries are subjected to the following conditions:

[TABLE]

A uniformly distributed load $q_{0}$ is exerted on the beam. The load is applied in increments of $\Delta q_{0}=50$ N/m until a maximum load of $1000$ N/m is reached. The maximum deflection, which occurs at the center of the beam, is recorded against the corresponding applied load. The results for both linear and nonlinear cases are plotted in Fig. 4a. The nonlinear deflections are smaller than the linear deflections at large loads because, as the load increases, the internal forces resisting the deformation increase in a nonlinear fashion.

For the three-point-bending case the boundaries are subjected to the following conditions:

[TABLE]

Here, instead of a uniformly distributed load, a point load $F_{0}$ is applied at the center of the beam. The point load is applied in increments of $\Delta F_{0}=50$ N until a maximum load of $1000$ N is reached. The maximum deflection, which occurs at the center of the beam, is recorded against the corresponding applied load. The results from the finite element model developed here for the 1-D equivalent single layer beam are compared with the 2-D FE results (see Fig. 4b). Note that ABAQUS uses a more complete Green-Lagrange strain tensor for the geometrically nonlinear beam element, whereas in the present finite element model developed in this paper the nonlinearity is included in the form of von Kármán strains.

4.3 Fixed-fixed hexagonal and Y-frame core beams

Here we consider two beams, one made of 48 hexagonal core unit cells and the other made of 30 Y-frame unit cells (see Fig. 2). Thus, the total length of the hexagonal core beam is $L=7.2$ m and the length of Y-frame core beam is $L=15.9$ m. Both the beams are subjected to a uniformly distributed load $q_{0}$. Fixed-fixed boundary conditions (55) are applied at the beam ends. For the hexagonal core beam the load is applied in increments of $\Delta q_{0}=75$ N/m until a maximum load of $1500$ N/m is reached, while for the Y-frame core beam the load is applied in increments of $\Delta q_{0}=5$ N/m until a maximum load of $100$ N/m is reached. The maximum transverse deflections, which occur at the beam centers, are recorded and plotted against the corresponding applied load in Figs. 5a and 5b.

4.4 Fixed-fixed corrugated core beam

A beam consisting of 30 corrugated core unit cells is considered. Since the length of each corrugated unit cell is $l=0.53$ m, the total length of the beam is $L=15.9$ m (see Fig. 2). The beam is subjected to fixed-fixed boundary conditions (55). A uniformly distributed load $q_{0}$ is applied on the beam. The load is applied in increments of $\Delta q_{0}=50$ N/m until a maximum load of $1000$ N/m is reached. The maximum vertical deflection of the beam is plotted against the corresponding applied load in Fig. 6a. The error in the maximum vertical deflection is calculated using,

[TABLE]

and is plotted against the applied load for both the linear and nonlinear cases in Fig. 6b.

Unlike the other structural cores, we see that the nonlinear deflections of the corrugated core beam, calculated using the finite element model developed for the 1-D equivalent single layer beam, are not in good agreement with the 2-D beam frame results all the way. This is due to the local buckling of the stretch-dominated corrugated core that occurs in the 2-D model as displayed in Fig. 7. The presented 1-D equivalent single layer model cannot account for this local buckling.

We also note that even though the lengths and heights of both corrugated core and Y-frame core beams are equal, the corrugated core beam is much stiffer than the Y-frame core beam. The maximum nonlinear deflection for the Y-frame core beam subjected to fixed-fixed boundary conditions is $88$ mm at a uniformly distributed load of $100$ N/m (see Fig. 5a), while for the corrugated core beam the maximum deflection for a uniformly distributed load of $100$ N/m is only $2.9$ mm (see Fig. 6a). The high stiffness of the corrugated core beam is attributed to its stretch-dominated behavior unlike Y-frame core which is bending-dominated. The corrugated core has a very high shear stiffness because of the fact that the elements (the Euler-Bernoulli beam elements within the core structure) of the corrugated core act essentially like rods and do not bend much, where as this is not the case in the Y-frame core. Although the elements of the part which resemble the corrugated core (the upper ‘V’ part of ‘Y’) in the Y-frame core do not exhibit lot of bending, the remaining part, consisting of lower element, undergoes significant bending.

5 Concluding remarks

In this paper, a displacement-based geometrically nonlinear finite element model for a micropolar Timoshenko beam which can undergo moderate rotations was developed. Appropriate reduced integration techniques were used to prevent the shear and membrane locking of the elements. The beam was used as an equivalent single layer model for lattice core sandwich beams. Constitutive equations were derived for four different lattice cores (hexagonal, corrugated, Y-frame and web-core) using a two-scale energy approach. The global bending results obtained from the developed 1-D finite element model were in good agreement with 2-D beam frame finite element results calculated using the finite element software ABAQUS.

Although the developed model predicts the global deflections with good accuracy, it does not capture the local buckling of stretch-dominated lattice cores. Nevertheless, the model may be extended to local buckling by allowing stiffness reduction inside the microscale unit cells within a computational multiscale finite element framework. Such considerations, and micropolar plates, are left for future studies.

Acknowledgements

The second author acknowledges that this work has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska–Curie grant agreement No 745770. The financial support is greatly appreciated. The authors also wish to acknowledge CSC – IT Center for Science, Finland, for computational resources (Abaqus usage).

Appendix A Transformation matrices

The displacement and strain transformation matrices in Eq. (8) are

[TABLE]

and

[TABLE]

respectively.