Asymptotic Homogenization of Auxetic-Core Sandwich Shells for the Isogeometric Analysis of a Rear Fuselage Structure with Embedded Engines

1. Introduction and Motivation

The aviation industry is under significant pressure to reduce its environmental impact. Boundary layer ingestion (BLI) is a promising way to reduce fuel consumption. By embedding the engines into the fuselage, the boundary layer forming on the fuselage is ingested into the engine, which re-energizes the aircraft’s wake. As a result, the wake drag is reduced and the overall efficiency is improved, with an expected energy saving of approximately 5–8% [2,3,4,5]. However, the engine inlet no longer receives a clean free stream but the distorted boundary layer of the fuselage, whose state changes considerably over the flight envelope and can provoke flow separation in different regions of the intake [6]. This could be addressed by morphing the engine inlet skin to adapt the aerodynamic surface for different flight conditions. Structurally, the skin around the inlet must carry the loads introduced by the embedded engines into the rear fuselage while limiting the transmission of engine vibrations and the resulting cabin noise.

A promising solution for these requirements is the use of auxetic metamaterials, which are engineered cellular structures with a negative Poisson’s ratio, where the cell geometry rather than the base material governs the deformation. Used as a sandwich core, such cells interact with the non-auxetic face sheets such that energy is stored and dissipated more effectively than in a conventional sandwich [7,8]. This behavior has been exploited for enhanced energy dissipation [9], high energy absorption [10] and passive vibration damping [11]. Since the response follows from the cell geometry, it can be tuned locally and very precisely through the deformation mechanism of the cell [12], without changing the base material.

In BLI configurations, such sandwiches with an auxetic core can offer favorable characteristics for the inlet region, as the face sheets provide a smooth surface to guide the flow while the tailored mechanical behavior is preserved beneath. Beyond this passive behavior, auxetic materials also open up the possibility of a morphing surface [13,14], which can offer significant advantages for the flow characteristics [15], and they provide stabilizing effects on curved and double-curved surfaces [16].

Despite these advantages, significant obstacles still prevent the widespread use of auxetic metamaterials in aerospace applications, one reason being the extensive numerical modeling effort. For components of the size of the proposed embedded engine inlet region, the complex cell geometry increases the system size considerably. The unit cells are smaller than the macroscopic component by orders of magnitude, making it prohibitively expensive to resolve their fine details. Many auxetic structures additionally feature notches, which further increase the modeling effort and the complexity. The resulting system size leads to high computational costs and makes it difficult to investigate real-world engineering problems.

Various approaches have been proposed to reduce the computational effort associated with such structures without compromising the accuracy required for analyzing the fine-scale properties described above. Notable examples include methods for reduced-fidelity modeling [17], several approaches involving machine learning [18], asymptotic homogenization [19] and two-scale computational homogenization [20]. These approaches differ mainly in their field of application. Reduced-fidelity modeling methods are applied to lattices with large auxetic cells, which allow for scale separation only to a limited extent, if at all, and therefore must retain a simplified representation of the cell geometry. Machine learning approaches, on the other hand, often involve inverse design or the generation of response surfaces that map the cell parameters to the effective properties by nonlinear interpolation and extrapolation, so that objective functions can be evaluated efficiently in areas such as structural optimization. Lastly, a wide variety of homogenization approaches are used to predict the overall behavior of large-scale structures with a clear separation of scales. Since the auxetic cells considered here are orders of magnitude smaller than the fuselage skin, the present work fits to this last category.

The novel contribution of this article is the homogenization of an auxetic-core sandwich, which makes the advantages listed above accessible while keeping the macroscale analysis computationally manageable. The framework is shown in Figure 1 and the two scales are described in Section 2. At the mesoscale, the ABD matrix of the auxetic sandwich section is obtained by numerical homogenization (Section 3). Two approaches are compared with a resolved reference model under tension, shear and bending as well as in a modal analysis (Section 4). At the macroscale, the rear fuselage is modeled with isogeometric analysis (IGA), which allows the structural analysis to be performed directly on the CAD description of the fuselage without a separate meshing step. The homogenized stiffness is then assigned to the skin of the inlet region, and the influence of the core topology and its orientation on the global deformation under a $2.5g$ pull-up maneuver is evaluated (Section 5). This novel approach, combining the design of auxetic unit cells at the mesoscale with isogeometric analysis for the structural analysis at the macroscale, enables the core topology to be used as a design variable at the aircraft scale.

2. Scale-Bridging Analysis Framework

2.1. Mesoscale: Auxetic Sandwich

In this article, a re-entrant auxetic cell is chosen as the representative unit cell, as shown in Figure 2. Its geometry is deliberately simple, essentially being a honeycomb with inverted diagonal struts. This allows it to be adapted to existing components with little geometric effort, while still producing an auxetic effect resulting from the re-entrant angle $\alpha$ of the diagonal struts. The corresponding dimensions are listed in Table 1, where w and h denote the total width and height of a cell, H the core height and t the strut thickness.

2.2. Macroscale: Structural Design Using Isogeometric Analysis

At the macroscale, the overall aim of the project is to identify unconventional stiffening topologies that minimize the structural weight of the rear fuselage using isogeometric analysis (IGA) [21]. IGA uses the same basis functions for the analysis that already describe the geometry, in this case non-uniform rational B-splines (NURBS). Geometry and simulation model are therefore no longer two separate representations. Instead, the CAD description of the outer fuselage shape is the analysis model itself, so that the meshing step and the geometric approximation it introduces are avoided [22]. Refinement reduces to knot insertion or degree elevation on that same description, which keeps the geometry exact at every refinement level and provides high solution accuracy with comparatively few degrees of freedom. Since the basis functions can be up to $C^{p-1}$-continuous across element boundaries, element formulations with continuity requirements beyond $C^0$ become accessible. In this study, isogeometric Kirchhoff–Love shell elements are used, which require $C^1$ continuity and therefore cannot be formulated with classical $C^0$-continuous Lagrange basis functions.

A further advantage of using IGA is the possibility of modeling the stiffeners as beams embedded into the fuselage shell, as proposed by Bauer et al. [23]. For the fuselage structural design, the stringer layouts are described as curves in the parameter space of the fuselage shell and mapped into physical space, so that the stiffeners introduce no degrees of freedom of their own. Instead, their structural response is expressed entirely through the shell. The stiffener layout is thereby decoupled from the shell discretization so that a stiffener can be moved, added or removed without remeshing. Together with the direct use of the CAD geometry, this makes the approach well suited for the early design stage, where the outer fuselage shape, the position of the embedded engines and the stiffening layout all change from one design variant to the next.

3. Homogenization Approaches

This section presents two different homogenization approaches that transfer the advantages of auxetic materials from the mesoscale (Section 2.1) to the macroscale (Section 2.2).

3.1. Computational Homogenization of a Periodic Unit Cell

The first approach consists of two steps, which are illustrated in Figure 3. First, the material properties of the auxetic sandwich core are determined on a unit cell acting as a representative volume element (RVE), using first-order computational homogenization. Second, the ABD matrix of the shell is calculated using classical laminate theory (CLT) by integrating the layer stiffnesses through the thickness.

3.1.1. Homogenization of the Sandwich Core Material

The goal of homogenization is to replace the periodic mesostructure, in this case the auxetic sandwich core structure, with an equivalent homogeneous medium that reflects the cell-averaged stiffness and deformation behavior. To this end, the displacement field of the unit cell is split into a linear macroscopic part and a periodic fluctuation, $\mathbf{u}\left(\mathbf{x}\right)=\overline{\epsilon }\mathbf{x}+\tilde{\mathbf{u}}\left(\mathbf{x}\right)$, where the fluctuation $\tilde{\mathbf{u}}$ takes the same value at any two opposite boundary points. Since the cell walls are prismatic in the z-direction, the deformation under periodic boundary conditions is uniform along z, and the homogenized properties do not depend on the height of the RVE. To reduce the model size, a slice of height $\overline{H}=H/4$ is therefore used instead of the full core height H. Opposite faces are separated by the period vectors $\mathbf{p}={(w,0,0)}^T$, ${(0,h,0)}^T$ or ${(0,0,\overline{H})}^T$ in the x-, y- and z-directions (see Figure 4). Subtracting the displacements of such a boundary pair therefore cancels the fluctuation and leaves

$$\mathbf{u}(\mathbf{x}+\mathbf{p})-\mathbf{u}\left(\mathbf{x}\right)=\overline{\epsilon }\mathbf{p}.$$

(1)

Consequently, the relative displacement of every boundary pair is prescribed by the macroscopic strain tensor $\overline{\epsilon }$ and the cell dimensions alone, independent of the position $\mathbf{x}$ on the face. The three resulting offsets are assigned to the displacements of one reference point per direction, ${\mathrm{RP}}_x$, ${\mathrm{RP}}_y$ and ${\mathrm{RP}}_z$, and enforced by linear constraint equations between each dependent node d and its independent periodic image n on the opposite face. This leads to the following equations

$$\begin{array}{cc}\hfill u_i^{d_x}-u_i^{n_x}-U_i^x& =0,\hfill \\ \hfill u_i^{d_y}-u_i^{n_y}-U_i^y& =0,\hfill \\ \hfill u_i^{d_z}-u_i^{n_z}-U_i^z& =0,\hfill \end{array}i=1,2,3,$$

(2)

with the reference-point displacements holding the columns of the macroscopic strain tensor scaled by the respective cell period,

$$\begin{array}{cc}\hfill {\mathbf{U}}^x& =\overline{\epsilon }\left(\begin{array}{c}w\\ 0\\ 0\end{array}\right),{\mathbf{U}}^y=\overline{\epsilon }\left(\begin{array}{c}0\\ h\\ 0\end{array}\right),{\mathbf{U}}^z=\overline{\epsilon }\left(\begin{array}{c}0\\ 0\\ \overline{H}\end{array}\right),\hfill \\ \hfill \overline{\epsilon }& =\left(\begin{array}{ccc}{\overline{\epsilon }}_{11}& \frac{1}{2}{\overline{\gamma }}_{12}& \frac{1}{2}{\overline{\gamma }}_{13}\\ \frac{1}{2}{\overline{\gamma }}_{12}& {\overline{\epsilon }}_{22}& \frac{1}{2}{\overline{\gamma }}_{23}\\ \frac{1}{2}{\overline{\gamma }}_{13}& \frac{1}{2}{\overline{\gamma }}_{23}& {\overline{\epsilon }}_{33}\end{array}\right).\hfill \end{array}$$

(3)

The RVE is discretized with shell elements having rotational degrees of freedom. Therefore, the boundary rotations must be periodic as well. Since there is no macroscopic curvature in first-order homogenization, the rotational fluctuation is periodic and the jump disappears entirely. For the pairs of opposite faces, the constraints read

$$\begin{array}{cc}\hfill {\phi }_i^{d_x}-{\phi }_i^{n_x}& =0,\hfill \\ \hfill {\phi }_i^{d_y}-{\phi }_i^{n_y}& =0,\hfill \\ \hfill {\phi }_i^{d_z}-{\phi }_i^{n_z}& =0,\hfill \end{array}i=4,5,6.$$

(4)

The green and blue edges in Figure 4 belong to two face pairs simultaneously, so nodes on these edges are paired in one direction only. Furthermore, one interior node $P_{RB}$ is fixed to remove the rigid-body translation. To compute the material properties, three tension and three shear load cases are applied, which are listed in Table 2. The reaction force ($RF$) at a reference point equals the resultant force transmitted across the corresponding cell face, and the areas of the RVE faces are $A_x=h\overline{H}$, $A_y=w\overline{H}$ and $A_z=wh$.

Uniaxial tension in x is imposed by prescribing $U_1^x=w$, which corresponds to a unit strain ${\overline{\epsilon }}_{11}=1$. The reference points ${\mathrm{RP}}_y$ and ${\mathrm{RP}}_z$ remain unloaded, and only their normal components $U_2^y$ and $U_3^z$ are left unconstrained (Table 2). As a result, the cell can deform freely in the transverse directions. The Young’s modulus and the Poisson’s ratios then follow from the effective stresses and the respective strains. For tension in the x-direction, the effective stresses and strains are

$$\begin{array}{cc}\hfill {\overline{\sigma }}_{11}& ={\displaystyle \frac{R{F}_1^x}{A_x}},{\overline{\epsilon }}_{11}={\displaystyle \frac{U_1^x}{w}},\hfill \\ \hfill {\overline{\epsilon }}_{22}& ={\displaystyle \frac{U_2^y}{h}},{\overline{\epsilon }}_{33}={\displaystyle \frac{U_3^z}{\overline{H}}}.\hfill \end{array}$$

(5)

From these relations, the Young’s modulus $E_1$ and the two Poisson’s ratio components ${\nu }_{12}$ and ${\nu }_{13}$ can be calculated as follows:

$$\begin{array}{cc}\hfill E_1& ={\displaystyle \frac{{\overline{\sigma }}_{11}}{{\overline{\epsilon }}_{11}}},{\nu }_{12}=-{\displaystyle \frac{{\overline{\epsilon }}_{22}}{{\overline{\epsilon }}_{11}}},{\nu }_{13}=-{\displaystyle \frac{{\overline{\epsilon }}_{33}}{{\overline{\epsilon }}_{11}}}.\hfill \end{array}$$

(6)

Tension in y and z is applied in the same manner and yields $E_2$, ${\nu }_{21}$, ${\nu }_{23}$, as well as $E_3$, ${\nu }_{31}$, ${\nu }_{32}$. For the shear cases, the strain jump is introduced at two reference points simultaneously. For shear in the x–y plane, this leads to $U_2^x=\frac{1}{2}w$ and $U_1^y=\frac{1}{2}h$. The shear stress is then available from either reference point,

$$\begin{array}{cc}\hfill {\overline{\tau }}_{12}& ={\displaystyle \frac{R{F}_1^y}{A_y}}={\displaystyle \frac{R{F}_2^x}{A_x}},\hfill \\ \hfill {\overline{\gamma }}_{12}& ={\displaystyle \frac{U_2^x}{w}}+{\displaystyle \frac{U_1^y}{h}},G_{12}={\displaystyle \frac{{\overline{\tau }}_{12}}{{\overline{\gamma }}_{12}}}.\hfill \end{array}$$

(7)

The remaining shear moduli follow analogously as

$$G_{13}={\displaystyle \frac{R{F}_1^z/A_z}{U_3^x/w+U_1^z/\overline{H}}},G_{23}={\displaystyle \frac{R{F}_2^z/A_z}{U_3^y/h+U_2^z/\overline{H}}}.$$

(8)

3.1.2. Computation of the ABD Matrix

The homogenized core is now treated as one thick ply in the composite stack (see Figure 3). Its orthotropic constants are the RVE results computed in the previous section. The ABD matrices used for the isogeometric shell model are computed according to CLT [22], which was implemented in the IGA solver of Kratos Multiphysics [24]. For each layer k, the reduced plane-stress stiffness is computed as

$${\mathbf{Q}}^{\left(k\right)}=\left(\begin{array}{ccc}{\displaystyle \frac{E_1}{1-{\nu }_{12}{\nu }_{21}}}& {\displaystyle \frac{{\nu }_{12}{E}_2}{1-{\nu }_{12}{\nu }_{21}}}& 0\\ {\displaystyle \frac{{\nu }_{12}{E}_2}{1-{\nu }_{12}{\nu }_{21}}}& {\displaystyle \frac{E_2}{1-{\nu }_{12}{\nu }_{21}}}& 0\\ 0& 0& G_{12}\end{array}\right).$$

(9)

The reduced plane-stress stiffness matrix is then rotated by the ply angle with the strain transformation matrix ${\mathbf{T}}_{\epsilon }$:

$${\overline{\mathbf{Q}}}^{\left(k\right)}={\mathbf{T}}_{\epsilon }^T\left({\theta }_k\right){\mathbf{Q}}^{\left(k\right)}{\mathbf{T}}_{\epsilon }\left({\theta }_k\right).$$

(10)

To compute the ABD matrices, the rotated matrices are integrated over the through-thickness coordinates $z_{k-1},z_k$ of the layer boundaries, measured from the laminate midplane:

$$\begin{array}{cc}\hfill \mathbf{A}& =\sum _k{\overline{\mathbf{Q}}}^{\left(k\right)}(z_k-z_{k-1}),\hfill \\ \hfill \mathbf{B}& =\frac{1}{2}\sum _k{\overline{\mathbf{Q}}}^{\left(k\right)}(z_k^2-z_{k-1}^2),\hfill \\ \hfill \mathbf{D}& =\frac{1}{3}\sum _k{\overline{\mathbf{Q}}}^{\left(k\right)}(z_k^3-z_{k-1}^3).\hfill \end{array}$$

(11)

Although the approach is simple to implement, it is based on assumptions that the real sandwich does not satisfy. First, the continuum assumption requires the cell to be small compared to the region it fills. In plane, this can be achieved with sufficiently small cells. Through the thickness, however, the core height H is of the same order as the cell dimensions and the core is bounded by the face sheets, whereas the periodic constraints in z treat it as an infinitely extended medium. Second, classical laminate theory reduces every layer to a plane stress state, so that only $E_1$, $E_2$, ${\nu }_{12}$ and $G_{12}$ enter the reduced stiffness ${\mathbf{Q}}^{\left(k\right)}$, while $E_3$, ${\nu }_{13}$, ${\nu }_{23}$, $G_{13}$ and $G_{23}$ are discarded. An alternative approach to calculating the ABD matrix that avoids these problems is asymptotic homogenization, which is described in the next section.

3.2. Asymptotic Homogenization

The second approach is based on the NIAH (Novel Implementation of Asymptotic Homogenization) method [25,26], which originates from the homogenization theory of plates with in-plane periodic microstructure, developed by Kalamkarov et al. [19]. Like the first approach, it requires only standard FE solutions. However, it homogenizes the complete periodic sandwich section (auxetic core and both face laminates), which directly results in the ABD matrix without an intermediate core material (see Figure 5).

The idea of asymptotic homogenization of periodic plate structures is based on defining a periodic unit cell whose thickness $\delta$ and in-plane dimensions are small compared to the dimensions of the component. The stresses and displacements are then expressed through an asymptotic expansion, which leads to equilibrium equations for the unit cell under in-plane periodic boundary conditions. These equations can be solved analytically or by transforming the problem into a finite element formulation, which, however, requires access to the internal structure of the FE solver.

The idea of NIAH is to circumvent this by expressing every step as a standard static solve of a commercial code [25,26]. First, a finite element model of the unit cell is constructed. Then, for each of the six unit macroscopic states (three membrane strains and three curvatures) the corresponding nodal displacement field $\chi$ is applied to every node of the unconstrained unit cell. The reaction force vectors for this nodal field can then be computed as $\mathbf{f}=\mathbf{K}\chi$. These forces are then applied to the unit cell with in-plane periodic constraints and free top and bottom surfaces, yielding the characteristic displacements $\tilde{\mathbf{a}}$ of every node. In a third solve, the characteristic displacements $\tilde{\mathbf{a}}$ are applied to the unconstrained cell, which returns the nodal reaction forces $\mathbf{P}=\mathbf{K}\tilde{\mathbf{a}}$ and allows for the computation of the effective material properties, which are the entries of the ABD matrix [26].

The original method is restricted to elements with translational degrees of freedom. Recently, Liu et al. developed an Abaqus plug-in that also allows the use of elements with rotational degrees of freedom [27]. This makes it directly applicable to the present cell meshed with S4 shell elements. The plug-in requires a unit cell that is mirror-symmetric about the coordinate axes, with the nodes on opposite boundary faces matching one to one. The shared cell walls in the x–z plane were therefore split into two halves along the cell end faces, one belonging to each neighboring cell. Each half is a shell of thickness $t/2$ whose nodes lie on the cell face and whose section is offset by half of its own thickness toward the cell interior (see Figure 6). Since the plug-in enforces periodicity for all six degrees of freedom of the boundary nodes, the two halves deform as one continuous wall.

3.3. Homogenization of the Mass

The modal analysis in Section 4 additionally requires the inertia of the section. Unlike for the stiffness, no cell problem has to be solved, since the kinetic energy is a plain volume average. The effective core density is the density of the wall material scaled by the relative density of the cell, ${\rho }^{*}={\rho }_{\mathrm{rel}}{\rho }_{\mathrm{wall}}$, and the core is treated as an additional homogeneous layer of the section. The relevant quantity for a shell is the mass per unit area,

$$I_0=\int \rho \mathrm{d}z=\sum _k{t}_k{\rho }_k=2t_{\mathrm{f}}{\rho }_{\mathrm{f}}+H{\rho }^{*},$$

(12)

where the density $\rho$ is integrated along the thickness coordinate z, which reduces to the sum of thickness $t_k$ times density ${\rho }_k$ over the layers k of the section, i.e., the two face laminates with thickness $t_{\mathrm{f}}$ and density ${\rho }_{\mathrm{f}}$ and the homogenized core with height H and density ${\rho }^{*}$. For the symmetric section the first moment $\int \rho z\mathrm{d}z$ vanishes, and rotary inertia is neglected. Hence, $I_0$ is the only inertia term of the thin-shell model.

4. Comparison at the Panel Level

4.1. Comparison of the ABD Matrices

First, the ABD matrices resulting from the two homogenization approaches of Section 3.1 and Section 3.2 are compared. For all comparisons below, the auxetic unit cell measures $10\times 10$ mm in plane, with a height of $H=8$ mm, a wall thickness of $t=0.7$ mm and a re-entrant angle of ${60}^{\circ }$ (see Figure 2). A stainless spring steel (EN 1.4310) with $E=185,000$ MPa, $\nu =0.30$ and $\rho =7.9$ g/cm³ is used as the auxetic core material.

The face laminates are quasi-isotropic ${[45/-45/0/90]}_{\mathrm{s}}$ layups with a ply thickness of $0.125$ mm, which gives a thickness of $1.0$ mm per face. The $0^{\circ }$ ply is oriented along the x-axis of the auxetic cell, which coincides with the x-axis of the panel in the examples below. An S2-glass/epoxy with $E_1=52,000$ MPa, $E_2=19,000$ MPa, ${\nu }_{12}=0.30$, $G_{12}=G_{13}=G_{23}=6,700$ MPa and $\rho =1.96$ g/cm³ is used for the plies ([28] [Table 1]).

For the homogenized auxetic core, the method of Section 3.1 yields $E_1=436.4$ MPa, $E_2=1324.4$ MPa, $G_{12}=50.4$ MPa, ${\nu }_{12}=-0.54$ and ${\nu }_{21}=-1.637$. The reduced plane-stress stiffness computed from these constants with Equation (9) is then inserted into the CLT routine of the IGA code as an additional core layer between the two face laminates, which results in a stack of 17 layers. The ABD matrix that follows from Equation (11) is shown on the left in Figure 7, next to the matrix obtained with the asymptotic homogenization of Section 3.2 on the right. The entries of $\mathbf{A}$ are given in MN/m and those of $\mathbf{D}$ in N m. The color of each block is scaled by its maximum absolute value with the same scale used for both approaches.

To isolate the core contribution, the face-laminate contribution (Equation (11)) is subtracted from both ABD matrices (Figure 8).

A comparison of the two core contributions shows that the first-order homogenization yields the more compliant core, with $A_{11}$ being approximately $27\%$ lower than the NIAH result. At the same time, it indicates a more pronounced auxetic response under in-plane loading, since $A_{12}=A_{21}$ is about $24\%$ larger in magnitude. Next, both homogenized models are compared with the resolved model under different load cases.

4.2. Square Plate Under Different Load Cases

4.2.1. Axial Tension

As a first example, a $500\times 500$ mm plate under a uniaxial tension of $N_{xx}=200$ N/mm is considered. The resolved model is analyzed with Abaqus using S4 shell elements and $7,083,614$ degrees of freedom. Figure A1(a) in Appendix A shows the boundary conditions of the resolved model. The axial resultant is introduced as a force of $F=100$ kN on a reference point that is tied to all nodes of the loaded end cross-section by linear constraint equations, so that this end remains plane. The opposite end is tied to a second reference point held at $u_x=0$. Additionally, three pins on the outer face remove the remaining rigid-body motion.

The homogenized models are analyzed with Kratos Multiphysics using isogeometric Kirchhoff–Love shell elements with three translational degrees of freedom per control point. The mesh of $5\times 5$ control points results in 75 degrees of freedom. As shown in Figure A1(b), the IGA model introduces the same resultant as a line load on the control-point row at one end and keeps the ends plane by multi-point constraints. One point on the edge opposite the loaded edge is constrained in all three translations, and two further points are constrained in the z-direction to remove the remaining rigid-body motion.

Figure 9 shows the displacement field on the top surface of the resolved plate (a), compared with the IGA results obtained with the ABD matrix from the first-order computational homogenization (b) and from the asymptotic homogenization (c). The field deviations are shown in Figure A5 in Appendix B. As anticipated from the ABD matrix comparison, the asymptotic homogenization leads to smaller deviations, while the first-order homogenization is too compliant and overestimates the auxetic behavior. The same is observed in the comparison of the axial and lateral strains in Table 3, where the effective strains are evaluated from the edge-averaged displacements $\langle ·\rangle$ as

$$\begin{array}{cc}\hfill {\overline{\epsilon }}_{xx}& ={\displaystyle \frac{{\langle u_x\rangle }_{x=L_x}-{\langle u_x\rangle }_{x=0}}{L_x}},\hfill \\ \hfill {\overline{\epsilon }}_{yy}& ={\displaystyle \frac{{\langle u_y\rangle }_{y=L_y}-{\langle u_y\rangle }_{y=0}}{L_y}},\hfill \\ \hfill {\nu }_{\mathrm{eff}}& =-{\displaystyle \frac{{\overline{\epsilon }}_{yy}}{{\overline{\epsilon }}_{xx}}},\hfill \end{array}$$

(13)

with the plate edge lengths $L_x=L_y=500$ mm.

4.2.2. In-Plane Shear

For the comparison under shear load, both models suppress the out-of-plane displacement $u_z$ along all four edges, and two corner points are pinned in plane to remove the remaining rigid-body motion (see Figure A2). The shear resultant of $N_{xy}=105$ N/mm is applied as a tangential edge load q that forms a closed shear flow around the four edges, so that the load is self-equilibrated and each edge carries a total force of $52.5$ kN.

Figure 10 shows the displacement field of the resolved plate (a), compared with the IGA results obtained with the ABD matrix from the first-order computational homogenization (b) and from the asymptotic homogenization (c). The field deviations are shown in Figure A6 in Appendix B. As already seen for the axial case, the first-order homogenization is too compliant, while the asymptotic homogenization matches the resolved response closely. The effective shear strain in Table 4 is evaluated as a pair of central differences between the four mid-edge points of the plate, denoted by the superscripts $\mathrm{t}$, $\mathrm{b}$, $\mathrm{l}$ and $\mathrm{r}$ for the top, bottom, left and right mid-edge points,

$$\begin{array}{cc}\hfill {\gamma }_{xy}& ={\displaystyle \frac{u_x^{\mathrm{t}}-u_x^{\mathrm{b}}}{L_y}}+{\displaystyle \frac{u_y^{\mathrm{r}}-u_y^{\mathrm{l}}}{L_x}},\hfill \\ \hfill A_{66,\mathrm{eff}}& ={\displaystyle \frac{N_{xy}}{{\gamma }_{xy}}},\hfill \end{array}$$

(14)

so that the two terms of ${\gamma }_{xy}$ are the discrete counterparts of $\partial u_x/\partial y$ and $\partial u_y/\partial x$. For both IGA models, the effective shear stiffness $A_{66,\mathrm{eff}}$ matches the $A_{66}$ entry of the corresponding ABD matrix in Figure 7.

4.2.3. Pure Bending

Pure bending is applied to the free plate with equal and opposite moments of $M={10}^5$ N mm at the two ends, and three corner pins on the outer face remove rigid-body motion in both models (see Figure A3). In the resolved model, all nodes of each end face are coupled to a reference point that carries the moment. Since the Kirchhoff–Love shells of the IGA model have no rotational degrees of freedom, the same moment is introduced there as work-conjugate nodal loads on the control points of each end row.

Figure 11 shows the displacement field of the resolved plate (a), compared with the IGA results obtained with the ABD matrix from the first-order computational homogenization (b) and from the asymptotic homogenization (c). The field deviations are shown in Figure A7 in Appendix B. The curvatures in Table 5 are obtained by fitting a paraboloid to the out-of-plane displacement field by least squares, and the bending Poisson’s ratio follows as ${\nu }_{\mathrm{bending}}=-{\kappa }_{yy}/{\kappa }_{xx}$. Since the lateral edges are moment-free, ${\nu }_{\mathrm{bending}}$ equals the ratio $D_{12}/D_{22}$ of the corresponding ABD matrix in Figure 7. The first-order homogenization is again too compliant, overestimating ${\kappa }_{xx}$ by $12.6\%$ and underestimating ${\nu }_{\mathrm{bending}}$ by $48.8\%$, while the asymptotic homogenization stays within $3\%$ of the resolved values.

4.2.4. Modal Analysis

For the modal analysis, the plate is simply supported along all edges, with two points pinned in plane to remove the remaining rigid-body motion (see Figure A4). Since the coarse $5\times 5$ control-point mesh of the static cases cannot represent the higher mode shapes, the discretization of the IGA models is refined for the modal analysis, resulting in the 1200 degrees of freedom listed in Table 6. Figure 12 compares the first four bending mode shapes and Table 6 the corresponding frequencies. Mode 2 of the resolved model is an in-plane mode and is omitted. All mode pairs reach a modal assurance criterion of $\mathrm{MAC}=1.00$, which rates the similarity of two mode shapes from 0 (no correspondence) to 1 (identical shapes), so that both homogenized models reproduce the mode shapes of the resolved plate. The first-order homogenization is again too compliant and therefore underestimates the frequencies by up to $6.6\%$, while the asymptotic homogenization stays within $0.4\%$. Since the asymptotic homogenization reproduces the resolved plate in all four comparisons, it is used for the fuselage analysis in Section 5.

5. Application to the Rear Fuselage Structure

The previous section has shown that shell stiffnesses homogenized using NIAH yield reliable results. This approach is now applied to the structural analysis of a rear fuselage with embedded engines. To this end, part of the stiffened outer skin near the embedded engine inlet is replaced by the homogenized auxetic sandwich material (marked in red in Figure 13).

Figure 14 shows the effect of the core topology and of its orientation within the patch on the deformation of the rear fuselage for a $2.5g$ pull-up maneuver at ultimate load. Figure 14 illustrates the displacement field for the auxetic sandwich core material introduced earlier (see Figure 2). With the approach proposed above, it is now possible to change the core topology, repeat the homogenization and apply the updated sandwich properties to the skin area. Three variants are compared with the reference: the reference cell rotated by ${45}^{\circ }$ within the patch, and re-entrant angles of $\alpha ={75}^{\circ }$ and $\alpha ={120}^{\circ }$, the latter corresponding to a conventional, non-auxetic honeycomb. Figure 14Figure 14 show the relative deviation $\mathrm{\Delta }$ of the displacement magnitude in percent with respect to the reference core topology in Figure 14. The core topology thus directly influences the global behavior of the rear fuselage. This opens up the possibility of optimizing the rear fuselage at both scales: at the macroscale by adapting the stiffener layout and at the mesoscale by adapting the sandwich core topology.

6. Summary and Outlook

The aim of this work was to lay the foundation for using auxetic metamaterials as the core of sandwich structures in order to exploit their advantages, such as passive vibration damping, high energy dissipation and tailorable deformation behavior. This could be particularly beneficial for the integration of engines into the fuselage for BLI configurations. However, a resolved simulation of the sandwich material is not feasible in this context.

Two homogenization approaches were investigated to determine the homogenized shell stiffnesses of auxetic sandwich materials. The first combines a first-order homogenization of the core, which is then treated as a laminate layer, with classical laminate theory to obtain the ABD matrix. The second approach is based on asymptotic homogenization. The resolved plate model and the two homogenized models were compared under tension, shear and bending as well as in a modal analysis. The asymptotic homogenization deviates from the resolved model by less than $6\%$, while the first-order homogenization predicts an overly compliant response and is therefore not suitable for the macroscale analysis. The homogenized shell stiffnesses can now be used to perform the global structural simulations directly on the CAD geometry with few degrees of freedom using isogeometric analysis. This was demonstrated on a rear fuselage with embedded engines, where the stiffened skin of the inlet region was replaced by the homogenized sandwich material. Varying the core topology and its orientation showed a direct influence on the global deformation of the rear fuselage under a $2.5g$ pull-up maneuver.

In the future, the homogenized transverse shear stiffnesses will also be determined, so that a shear-deformable shell theory can be used, which is particularly important for thick sandwich structures. Furthermore, these methods will be used for a simultaneous optimization of the stiffener topology at the macroscale and the auxetic core topology at the mesoscale.