Optimal Design of CNT-Nanocomposite Nonlinear Shells

Abstract

Carbon nanotube/polymer nanocomposite plate- and shell-like structures will be the next generation lightweight structures in advanced applications due to the superior multifunctional properties combined with lightness. Here material optimization of carbon nanotube/polymer nanocomposite beams and shells is tackled via ad hoc nonlinear finite element schemes so as to control the loss of stability and overall nonlinear response. Three types of optimizations are considered: variable through-the-thickness volume fraction of random carbon nanotubes (CNTs) distributions, variable volume fraction of randomly oriented CNTs within the mid-surface, aligned CNTs with variable orientation with respect to the mid-surface. The collapse load, which includes both limit points and deformation thresholds, is chosen as the objective/cost function. An efficient computation of the cost function is carried out using the Koiter reduced order model obtained starting from an isogeometric solid-shell model to accurately describe the point-wise material distribution. The sensitivity to geometrical imperfections is also investigated. The optimization is carried out making use of the Global Convergent Method of Moving Asymptotes. The extensive numerical analyses show that varying the volume fraction distribution as well as the CNTs orientation can lead to significantly enhanced performances towards the loss of elastic stability making these lightweight structures more stable. The most striking result is that for curved shells, the unstable postbuckling response of the baseline material can be turned into a globally stable response maintaining the same amount of nanostructural reinforcement but simply tailoring strategically its distribution.

1. Introduction

Thin-walled, lightweight composite structures are commonly used in a wide range of engineering applications, particularly in aerospace engineering, where they are often employed as primary structural components. Due to the high strength-to-weight ratio, the mechanical response is dominated by buckling and turns out to be mainly influenced by two factors: the geometry and the elastic properties. While the former is often imposed by the structural functionality and only little variations are possible, the spatial distribution of the material properties (e.g., fiber orientations in the layups) can be easily tailored in composite shells. Consequently, an efficient optimization process of the material distribution is required to obtain the desired structural response, usually defined in terms of deflections and load-carrying capacity. Many manufacturing options are also available to fine-tune the stiffness and the onset of buckling: grid stiffeners [1], multi-layered and variable thickness composites [2], variable angle tows (VATs) [3].

A promising direction in the context of material optimization is that offered by nanostructured materials which can exhibit multifunctional properties and are thus prone to more advanced multi-objective optimizations. In this field, nanocomposite materials made of thermosetting or thermoplastic polymers integrated with carbon nanotubes (CNTs) are currently subject to intense developments due to their superior mechanical/electrical/thermal performance, electromagnetic shielding or energy storage capacity [4,5,6,7]. Among other attractive properties exhibited by nanocomposites, a relatively high strength-to-weight ratio, unique damping capability [8,9] and high fatigue tolerance make them ideal candidates for whole new classes of multifunctional composite structures (e.g., high-performance vehicles, aerostructures and devices). Previous works addressed the optimal design of multilayer nanocomposites by fine-tuning the nonlinear interfacial properties regulating the CNT/polymer stick-slip in each of the dedicated layers to achieve maximum storage and damping capabilities [10,11,12]. The optimization problem was restricted to two- and three-layer nanocomposites embedding in selected polymers, CNTs with tunable properties such as the interfacial shear strength and the CNT volume fraction. A robust gradient-free optimization algorithm was developed employing the family of differential evolution optimizers. The objective function was chosen to be the product of the average damping ratio and stored energy of the multilayer nanocomposite within a predefined range of admissible strains.

Moreover, as known, optimal properties in composites are usually sought also by controlling the orientation of the fibers in each layer, since the fiber orientation significantly affects the stiffness distribution, hence, the load-carrying capability or the elastic limits states (see, for instance, [13,14]). Along the same lines, optimization of the CNT distribution within a composite can be achieved to maximize the buckling loads or to change the global stability of nanocomposite shells. Many optimization strategies proposed in the literature use the linearized buckling load as the objective function of the design. However, in this case, structures may suffer another elastic limit state (i.e., static bifurcation) known as buckling mode interaction, which leads to an unstable post-critical behavior [15] and a high sensitivity to imperfections, resulting in a deterioration of their load-bearing capacity due to geometrical, load and material deviations. For this reason, a more reliable design, which takes into account the geometrically nonlinear behavior, has also been investigated in previous works. In this framework, the collapse load of the structure can be defined as the first limit load, for the unstable cases, or as the load magnitudes giving rise to deformations which compromise the usability, taking into account the stiffness reduction that typically characterizes the post-buckling regime. Optimizing the post-buckling behavior in terms of collapse load is, however, a challenging task. In fact, a suitable mechanical model and its discrete approximation are required to describe with acceptable accuracy the geometry, the boundary conditions and the mechanical behavior. This means that the structural response is generally described by a high number of discrete nonlinear equations, whose solution describes the equilibrium path.

The Riks arc-length strategy [16,17] is a standard tool for path following the solutions of a set of nonlinear equations. Although this approach is effective for assigned data, it is not suitable for an optimization process, which requires the evaluation of the equilibrium path for any change in the design variables, and for an imperfection sensitivity analysis, because a single run is too time-consuming with current CPUs. Promising generalizations of the path following strategy have been presented [18,19,20] with the aim of performing parametric studies in a more efficient way.

In the optimization strategy presented in [21], the collapse load is evaluated by a nonlinear finite element (FE) buckling problem. The algorithm is extended in [22,23] in order to take into account the worst geometrical imperfection case. An interesting way of analyzing slender structures is offered by strategies based on Koiter’s theory of elastic stability [24]. They make use of an asymptotic expansion of the equilibrium equations which allows the description of the initial post-critical behavior in terms of some variables related to the slope and curvature of the bifurcated branches [25].

More recently, a solution algorithm based on Koiter’s theory implemented within a Finite Element environment was proposed in [26,27]. It allows structures to be optimized with general geometries, loading and boundary conditions. Moreover, the strategy exhibits good levels of accuracy in predicting the initial postbuckling response of several structures due to a multi-modal asymptotic expansion which accounts also for nonlinear buckling modal interactions [28]. The strategy is also capable of efficiently providing the worst equilibrium path through a statistical estimation of the worst-case imperfection, which is assumed to belong to the space of the buckling modes of the structure under consideration. A hybrid solution strategy, referred to as the Koiter-Newton approach, was further investigated in [29,30].

Despite the difficulties associated with the prediction of the nonlinear behavior, another challenge is the solution to the optimization problem. This is always expressed as a nonlinear, nonconvex mathematical programming problem whose solution is generally computationally expensive and extremely difficult because of the likely presence of multiple local minima. This is indeed the most penalizing aspect of the analysis. Among others, frequently employed solution strategies are the random search methods [31], genetic algorithms [32] and gradient-based techniques such as the method of moving asymptotes [33] or sequential linear programming [34].

This work addresses the optimization problem of nanocomposite shell-like structures with variable CNTs distributions. Three types of optimizations are tackled: (a) through-the-thickness CNT distribution in terms of volume fractions, (b) randomly oriented CNT distributions across the mid-surface in terms of volume fraction, (c) distribution of aligned CNTs across the mid-surface in terms of orientation. The collapse load, including limit points and deformation limits, is taken as the objective function. Its efficient estimate is carried out using a reduced order model (Koiter’s method) obtained starting from an isogeometric solid-shell model capable of accurately modeling the variable material distribution. The sensitivity to geometrical imperfections is also addressed as in [26]. The optimization is performed using the Global Convergent Method of Moving Asymptotes (GCMMA). A numerical investigation is carried out to assess how varying the volume fraction distribution, as well as the CNTs orientation, affects the stability of nanocomposite shells. It is shown that an unprecedented tuning of the shells stability can be achieved in different ways, either affecting the local bifurcation behavior (e.g., shifting the buckling loads to higher values) or by affecting the global behavior (e.g., suppressing the snap-through instability).

The paper is organized as follows: Section 2 describes the solid-shell model for elastic shell structures; Section 3 introduces the isogeometric NURBS-based discretization technique; the constitutive equations for a polymeric matrix reinforced with CNTs are reported in Section 4; Section 5 formulates the nonlinear optimization problem for variable distributions of CNTs in the shell domain; the numerical optimization strategy based on Koiter method and GCMMA is presented in Section 6; extensive numerical analyses are carried out in Section 7; conclusions are drawn in Section 8.

2. Solid-Shell Model

This section describes the main equations of the solid-shell FE formulation (see [35]) used to construct the discrete model. The outset is the 3D Cauchy continuum whose geometry of deformation is described by the Green–Lagrange strains. By employing a solid-shell concept, a linear through-the-thickness interpolation is assumed for the kinematic unknowns. Euclidean vectors will be denoted by boldface, italic, lower case letters, algebraic vectors in ${\mathbb{R}}^3$ will be denoted by boldface, non-italic letters; tensors by boldface, italic, upper case letters, matrices by boldface, non-italic, upper case letters.

Stored and Complementary Energy

Convective curvilinear shell coordinates $\zeta =({\zeta }_1,{\zeta }_2,{\zeta }_3)$ are employed, with $({\zeta }_1,{\zeta }_2)$ representing mid-surface coordinates and ${\zeta }_3\in [-\frac{h}{2},\frac{h}{2}]$ being the thickness coordinate with h the shell thickness. The position vector $\mathit{p}\left(\mathit{\zeta }\right)$ of material points in the current configuration is given in terms of their position vector $\mathit{x}\left(\mathit{\zeta }\right)$ in the reference configuration and the displacement $\mathit{u}\left(\mathit{\zeta }\right)$,

$$\mathit{p}\left(\mathit{\zeta }\right)=\mathit{x}\left(\mathit{\zeta }\right)+\mathit{u}\left(\mathit{\zeta }\right).$$

(1)

Note that while $(\mathit{x},\mathit{p},\mathit{u})$ are treated as Euclidean vectors, $\mathit{\zeta }$ is a vector of ${\mathbb{R}}^3.$ The covariant basis vectors in the undeformed configuration are obtained from the corresponding partial derivatives $G_i=\mathit{x}{,}_i$ of the position vectors $\mathit{x}$, where $\left(\right){,}_i$ indicates partial differentiation with respect to the i-th component of $\mathit{\zeta }$. By letting ${\mathit{G}}^i$ denote the contravariant basis so that ${\mathit{G}}_i·{\mathit{G}}_j={\delta }_i^j$ with ${\delta }_i^j$ the Kronecker delta and $(·)$ the dot product, the Green–Lagrange strain tensor can be expressed as

$$\mathit{E}={\overline{E}}_{ij}{\mathit{G}}^i\otimes {\mathit{G}}^j,{\overline{E}}_{ij}=\frac{1}{2}\left(\mathit{x}{,}_i·\mathit{u}{,}_j+\mathit{u}{,}_i·\mathit{x}{,}_j+\mathit{u}{,}_i·\mathit{u}{,}_j\right),$$

(2)

where $(\otimes )$ indicates the tensor product.

Assuming a linear through-the-thickness interpolation, the position vector is expressed as

$$\mathit{x}\left(\mathit{\zeta }\right)={\mathit{x}}_0({\zeta }_1,{\zeta }_2)+\frac{2{\zeta }_3}{h}{\mathit{x}}_n({\zeta }_1,{\zeta }_2)$$

(3)

where ${\mathit{x}}_0:=\frac{1}{2}\left(\mathit{x}\left({\mathit{\zeta }}_{+}\right)+\mathit{x}\left({\mathit{\zeta }}_{-}\right)\right)$ and ${\mathit{x}}_n:=\frac{1}{2}\left(\mathit{x}\left({\mathit{\zeta }}_{+}\right)-\mathit{x}\left({\mathit{\zeta }}_{-}\right)\right)$, with ${\mathit{\zeta }}_{+}=({\zeta }_1,{\zeta }_2,\frac{h}{2})$ and ${\mathit{\zeta }}_{-}=({\zeta }_1,{\zeta }_2,-\frac{h}{2})$. Similarly, the displacement field is described as

$$\mathit{u}={\mathit{u}}_0({\zeta }_1,{\zeta }_2)+\frac{2{\zeta }_3}{h}{\mathit{u}}_n({\zeta }_1,{\zeta }_2)$$

(4)

with ${\mathit{u}}_0:=\frac{1}{2}\left(\mathit{u}\left({\mathit{\zeta }}_{+}\right)+\mathit{u}\left({\mathit{\zeta }}_{-}\right)\right)$ and ${\mathit{u}}_n:=\frac{1}{2}\left(\mathit{u}\left({\mathit{\zeta }}_{+}\right)-\mathit{u}\left({\mathit{\zeta }}_{-}\right)\right)$ being the coordinates of the upper and lower surfaces of the shell. The independent strain components in Equation (2) are collected in the six-dimensional strain vector $ϵ={[E_{11},E_{22},2E_{12},E_{33},2E_{23},2E_{13}]}^T$ and linearized with respect to ${\zeta }_3$ as

$$ϵ\approx \left[\begin{array}{c}e({\zeta }_1,{\zeta }_2)+{\zeta }_3\mathit{\chi }({\zeta }_1,{\zeta }_2)\\ E_{33}\left({\mathit{\zeta }}_0\right)\\ \mathit{\gamma }({\zeta }_1,{\zeta }_2)\end{array}\right]$$

(5)

with ${\mathit{\zeta }}_0=({\zeta }_1,{\zeta }_2,0)$ and the membrane strain vector $\mathit{e}$, the curvature vector $\mathit{\chi }$, and the transverse shear strains vector $\mathit{\gamma }$ given by

$$\mathbf{e}({\zeta }_1,{\zeta }_2)=\left[\begin{array}{c}{E}_{11}\left({\mathit{\zeta }}_0\right)\\ E_{22}\left({\mathit{\zeta }}_0\right)\\ 2E_{12}\left({\mathit{\zeta }}_0\right)\end{array}\right],\mathit{\chi }({\zeta }_1,{\zeta }_2)=\left[\begin{array}{c}{E}_{11}{,}_3\left({\mathit{\zeta }}_0\right)\\ E_{22}{,}_3\left({\mathit{\zeta }}_0\right)\\ 2E_{12}{,}_3\left({\mathit{\zeta }}_0\right)\end{array}\right],\mathit{\gamma }({\zeta }_1,{\zeta }_2)=\left[\begin{array}{c}2E_{23}\left({\mathit{\zeta }}_0\right)\\ 2E_{13}\left({\mathit{\zeta }}_0\right)\end{array}\right].$$

The constitutive equations, in Voigt notation, of the 3D continuum are assumed to be those of a transversally isotropic material with the axis of transverse isotropy collinear with the axis of the CNTs. Hence, they can be rewritten as the following block matrix:

$$\mathbf{L}=\left[\begin{array}{ccc}{\mathbf{L}}_m& 0& \mathbf{0}\\ \mathbf{0}& L_{33}& \mathbf{0}\\ \mathbf{0}& 0& {\mathbf{L}}_s\end{array}\right],$$

by neglecting the coupling between the membrane strains $\mathit{e}$ and the thickness strain $E_{33}$ and blocks defined by Equation (5). The block matrix ${\mathbf{L}}_s$ describes the shear elastic constants. The membrane-related block matrix is evaluated assuming plane stress conditions in order to avoid thickness locking, while the transverse elastic constant $L_{33}$ is maintained to avoid zero-energy modes (thickness stretch).

By denoting with V the region occupied by the shell in the reference configuration and $\mathsf{\Omega }$ the area of its mid-plane and performing the integration over the thickness, we obtain

$${\int }_V{ϵ}^{\top }\mathbf{L}ϵdV={\int }_{\mathsf{\Omega }}{\int }_{-h/2}^{h/2}{ϵ}^{\top }\mathbf{L}ϵd{\zeta }_{3}d\mathsf{\Omega }={\int }_{\mathsf{\Omega }}{\mathit{\epsilon }}^{\top }\mathbf{C}\mathit{\epsilon }d\mathsf{\Omega },$$

(6)

where vector $\mathit{\epsilon }({\zeta }_1,{\zeta }_2):={\left[e^{\top },E_{33},{\mathit{\chi }}^{\top },{\mathit{\gamma }}^{\top }\right]}^{\top }$ collects the generalized strains and $\mathbf{C}$ is the generalized constitutive matrix expressed as

$$\mathbf{C}=\left[\begin{array}{cccc}{\mathbf{C}}_{ee}& \mathbf{0}& {\mathbf{C}}_{e\mathit{\chi }}& \mathbf{0}\\ \mathbf{0}& C_{33}& \mathbf{0}& \mathbf{0}\\ {\mathbf{C}}_{e\mathit{\chi }}^T& \mathbf{0}& {\mathbf{C}}_{\mathit{\chi }\mathit{\chi }}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& \mathbf{0}& {\mathbf{C}}_{\mathit{\gamma }\mathit{\gamma }}\end{array}\right]$$

(7)

with

$${\mathbf{C}}_{ee}={\int }_{-\frac{h}{2}}^{\frac{h}{2}}{\mathbf{L}}_{m}d{\zeta }_3,{\mathbf{C}}_{e\mathit{\chi }}={\int }_{-\frac{h}{2}}^{\frac{h}{2}}{\zeta }_3{\mathbf{L}}_{m}d{\zeta }_3,{\mathbf{C}}_{\mathit{\chi }\mathit{\chi }}={\int }_{-\frac{h}{2}}^{\frac{h}{2}}{\zeta }_3^2{\mathbf{L}}_{m}d{\zeta }_3,$$

$$C_{33}={\int }_{-\frac{h}{2}}^{\frac{h}{2}}{L}_{33}d{\zeta }_3,{\mathbf{C}}_{\mathit{\gamma }\mathit{\gamma }}={\int }_{-\frac{h}{2}}^{\frac{h}{2}}{\mathbf{L}}_{s}d{\zeta }_3.$$

3. Isogeometric Solid-Shell Model

The continuum solid-shell model is discretized by using NURBS functions. In particular, according to IGA, the same interpolation is used for the geometry and displacements [35,36].

NURBS basics. A B-Spline curve is represented as

$$\mathit{g}\left(\zeta \right)=\sum _{i=1}^n{N}_i^p\left(\zeta \right){\mathit{y}}_i=\mathit{N}\left(\zeta \right)·\mathit{y},$$

(8)

where ${\mathit{y}}_i$ ($i=1\cdots n$) are control points and $N_i^p\left(\zeta \right)$ are the set of B-Spline basis functions taken as piecewise polynomial functions of order p. The latter are defined by a set of nondecreasing real numbers $\mathsf{\Xi }=[{\zeta }^1,{\zeta }^2,\cdots ,{\zeta }^{n+p+1}]$ known as knot vectors. B-Spline basis functions are calculated recursively by using

$$N_i^p\left(\zeta \right)=\frac{\zeta -{\zeta }^i}{{\zeta }^{i+p}-{\zeta }^i}{N}_i^{p-1}\left(\zeta \right)+\frac{{\zeta }^{i+p+1}-\zeta }{{\zeta }^{i+p+1}-{\zeta }^{i+1}}{N}_{i+1}^{p-1}\left(\zeta \right),$$

(9)

for $p\ge 1$ and starting from piecewise constant functions ($p=0$) defined as

$$N_i^0\left(\zeta \right)=\left\{\begin{array}{cc}\hfill 1& ,\mathrm{if}{\zeta }^i\le \zeta \le {\zeta }^{i+1}\hfill \\ \hfill 0& ,\mathrm{otherwise}.\hfill \end{array}\right.$$

The regularity of B-Spline basis functions is given by $r=p-s$, where p and s are the order used for the basis functions and the multiplicity of the knot ${\zeta }_i$, respectively. NURBS functions are obtained by a projective transformation of the B-splines by extending Equation (8) with the following shape functions:

$$R_i^p\left(\zeta \right)=\frac{N_i^p\left(\zeta \right)w_i}{{\sum }_i^n{N}_i^p\left(\zeta \right)w_i}.$$

(10)

It is worth noting that all properties of B-Splines are retained and, in particular, a B-Spline is retrieved when all weights are equal.

By applying the tensor product, the NURBS surface is constructed in a way similar to Equation (8) as

$$\mathbf{q}({\zeta }_1,{\zeta }_2)=\sum _{i=1}^n\sum _{j=1}^m{R}_i^p\left({\zeta }_1\right)M_j^q\left({\zeta }_2\right){\mathbf{Y}}_{ij}=\mathbf{N}({\zeta }_1,{\zeta }_2)\mathbf{Y},$$

(11)

where $\mathcal{R}=[{\zeta }_1^1,{\zeta }_1^2\cdots {\zeta }_1^{n+p+1}]$ and $\mathcal{M}=[{\zeta }_2^1,{\zeta }_2^2\cdots {\zeta }_2^{m+q+1}]$ are two knot vectors, $R_i^p$ and $M_j^q$ are the one-dimensional basis functions over them and ${\mathbf{Y}}_{ij}$ defines a set of $n\times m$ control points. The tensor product between the knot vectors $\mathcal{R}$ and $\mathcal{M}$ defines a mesh of quadrilateral isogeometric elements.

Weights, as well as control points of the initial geometry, are provided by the CAD model while suitable algorithms exist for the refinement required to approximate the unknown solution [37]. The geometry is represented exactly regardless of the adopted mesh.

Isogeometric interpolation. In this subsection, the discrete isogeometric model used within the optimization strategy is summarized. The geometry is described by NURBS interpolation functions as

$$\mathbf{x}\left(\zeta \right)={\mathbf{N}}_u\left(\mathit{\zeta }\right){\mathbf{x}}_e$$

(12)

where ${\mathbf{x}}_e=[{\mathbf{x}}_{0e},{\mathbf{x}}_{ne}]$ collects the control points of the geometry corresponding to ${\mathit{x}}_0$ and ${\mathit{x}}_n$, respectively. The matrix ${\mathbf{N}}_u\left(\mathit{\zeta }\right)$ collects the interpolation functions

$${\mathbf{N}}_u\left(\mathit{\zeta }\right):=\left[\begin{array}{c}\mathbf{N}({\zeta }_1,{\zeta }_2),\frac{2{\zeta }_3}{h}\mathbf{N}({\zeta }_1,{\zeta }_2)\end{array}\right]$$

(13)

where $\mathbf{N}({\zeta }_1,{\zeta }_2)$ are bi-dimensional NURBS (11) functions of the mid-surface coordinates only.

By following the isogeometric concept, the displacement field is interpolated using the same shape functions of the geometry

$$\mathbf{u}\left(\mathit{\zeta }\right)={\mathbf{N}}_u\left(\mathit{\zeta }\right){\mathbf{u}}_e$$

(14)

where ${\mathbf{u}}_e=[{\mathbf{u}}_{0e},{\mathbf{u}}_{ne}]$ collects the control points for the displacement fields ${\mathbf{u}}_0$ and ${\mathbf{u}}_n$.

The Green–Lagrange strains in Equation (5), upon considering Equations (12) and (14), become

$$\mathit{\epsilon }({\zeta }_1,{\zeta }_2,{\mathbf{u}}_e)=\left(\mathcal{L}({\zeta }_1,{\zeta }_2)+\frac{1}{2}\mathcal{Q}({\zeta }_1,{\zeta }_2,{\mathbf{u}}_e)\right){\mathbf{u}}_e,$$

(15)

where $\mathcal{L}({\zeta }_1,{\zeta }_2):=\mathcal{Q}({\zeta }_1,{\zeta }_2,{\mathbf{x}}_e)$ and $\mathcal{Q}({\zeta }_1,{\zeta }_2,{\mathit{u}}_e)$ has a linear dependence from ${\mathbf{u}}_e$ and its expression can be found in [35].

Stored energy and equilibrium path. The stored energy of the shell can be evaluated using a numerical integration as

$$\Phi =\sum _e{\Phi }_e,{\Phi }_e=\frac{1}{2}\sum _g\left({\mathit{\epsilon }}_g{\left({\mathbf{u}}_e\right)}^{\top }{\mathbf{C}}_g{\mathit{\epsilon }}_g\left({\mathbf{u}}_e\right)\right)w_g,$$

(16)

where e denotes the FE, g indicates the integration point and $w_g$ is the corresponding weight. By taking advantage of the high continuity of the NURBS functions, patch-wise integration schemes can be adopted, thereby reducing the number of integration points. Moreover, well-tuned patch-wise reduced scheme can avoid locking (for more details on these aspects, see [35,36]).

The system of discrete equilibrium equations is then obtained through enforcement of the stationarity of the total potential energy according to

$$\mathbf{r}(\lambda ;\mathbf{u})=\frac{\partial \Phi }{\partial \mathbf{u}}-\lambda \overline{\mathbf{f}}=\mathbf{n}\left(\mathbf{u}\right)-\lambda \overline{\mathbf{f}}=\mathbf{0},$$

(17)

where $\mathbf{r}$ is the residual vector, $\mathbf{n}\left(\mathbf{u}\right)$ is the vector of generalized stress resultants (i.e., restoring forces), $\overline{\mathbf{f}}$ is the load vector per unit multiplier, $\mathbf{u}$ collects the discrete displacement variables of the isogeometric model and $\lambda$ is the load multiplier. The solutions of Equation (17) define the equilibrium paths of the structure in the $\mathbf{u}-\lambda$ space.

4. Constitutive Formulation for CNT Nancomposite Shells

Our aim is to tailor the nanostructured material distribution within the shell that leads to the optimization of a certain property of the nonlinear equilibrium path obtained (e.g., maximizing the collapse load) by solving the equilibrium problem later expressed by Equation (27). Before delving into further details of the material tailoring problem, we pause to discuss the constitutive model for nanocomposites with randomly oriented or perfectly aligned CNTs (Figure 1).

A compact expression of the effective elastic tensor for the constitutive law of the 3D continuum is given by

$$\overline{\mathit{L}}={\mathit{L}}_M+{\varphi }_C\langle \mathit{B}\rangle {\left[(1-{\varphi }_C)\mathit{I}+{\varphi }_C\langle \mathit{A}\rangle \right]}^{-1}$$

(18)

with

$$\mathit{B}=〚\mathit{L}〛\mathit{A},〚\mathit{L}〛={\mathit{L}}_C-{\mathit{L}}_M,\mathit{A}={\left(\mathit{I}+\mathit{S}{\mathit{L}}_M^{-1}〚\mathit{L}〛\right)}^{-1},$$

where $\mathit{I}$ is the identity tensor, ${\mathit{L}}_C$ and ${\mathit{L}}_M$ are the elastic tensors of the CNT inclusions and the matrix, respectively, $〚\mathit{L}〛$ is the elastic mismatch, and $\mathit{S}$ is the Eshelby tensor.

Tensor $\mathit{B}$ must be transformed in order to account for a generic orientation of the material frame $\{{\overline{e}}_1,{\overline{e}}_2,{\overline{e}}_3\}$ according to

$${\overline{B}}_{ijkl}=c_{i\mathit{p}}{c}_{jq}{c}_{kr}{c}_{ls}{B}_{pqrs}$$

(19)

where the explicit expression of the transformation coefficients $c_{ip}$ is reported in [13]. Subsequently, the notation $\langle ·\rangle$ indicates an averaging over the range of the CNTs orientations and their expression is given in [13]. Two cases are considered here:

(i)

randomly orientated CNTs

$$\langle \mathit{B}\rangle =\frac{1}{8{\pi }^2}{\int }_0^{2\pi }{\int }_0^{\pi }{\int }_0^{2\pi }\overline{\mathit{B}}sin\vartheta d\phi d\vartheta d\beta$$

(ii)

CNTs alignment along ${\overline{e}}_1$

$$\langle \mathit{B}\rangle =\frac{{\int }_0^{2\pi }{\int }_0^{\pi }{\int }_0^{2\pi }\overline{\mathit{B}}f(\phi ,\beta )sin\vartheta d\phi d\vartheta d\beta }{{\int }_0^{2\pi }{\int }_0^{\pi }{\int }_0^{2\pi }f(\phi ,\beta )sin\vartheta d\phi d\vartheta d\beta },$$

where $(\phi ,\vartheta ,\beta )$ are the Euler angles providing the CNTs orientation with respect to a fixed frame, $f(\phi ,\beta )$ is the orientation distribution function which, for the case of fibers aligned along the ${\overline{E}}_1$ axis, becomes $f(\phi ,\beta )=\delta (\phi -0)\delta (\beta -0)$ (see [13,14] for details).

The constitutive relation of the shell is obtained as $\mathbf{L}=\mathbf{R}\overline{\mathbf{L}}{\mathbf{R}}^{\top }$ where $\mathbf{R}=\mathbf{R}\left(\theta \right)$ is the rotation matrix with $\theta$ indicating the angle between ${\overline{e}}_1$ and the $e_1$ axis of the shell, while $e_3:={\overline{e}}_3$.

It is also possible to enhance the accuracy of the Eshelby–Mori–Tanaka model by taking into account the actual CNTs aspect ratio and the CNTs macrodispersion in the considered nanocomposites [14].

5. Postbuckling Optimization of CNT Nanocomposite Shells

Several types of material optimizations will be sought, namely, (1) through-the-thickness optimization of the aligned CNTs volume fraction; (2) optimization of randomly orientated CNTs volume fraction; (3) optimization of the in-plane CNTs orientation. We will discuss each of these problems next.

5.1. Through-the-Thickness Optimization of the Aligned CNTs Volume Fraction

In the first case, the optimal distribution of CNTs parameters is obtained by changing the through-the-thickness CNTs distribution across the shell mid-surface.

The volume fraction of CNTs is assumed to be distributed as

$${\varphi }_C\left(\mathbf{\zeta }\right)=\left[1+\frac{2}{h}a({\zeta }_1,{\zeta }_2){\zeta }_3\right]{\varphi }_C^{*},$$

which means that the through-the-thickness average ${\varphi }_C^{*}$ of ${\varphi }_C\left[\mathbf{\zeta }\right]$ is constant over the shell domain, while $a({\zeta }_1,{\zeta }_3)$ is a function of the through-the-thickness distribution such that $-1\le a({\zeta }_1,{\zeta }_2)\le 1$. In this way, we limit the volume fraction variability to lie in the range $0\le {\varphi }_C\le 2{\varphi }_C^{*}$. The constitutive matrix $\mathit{L}=\overline{\mathit{L}}$ is that obtained for the case of uniformly aligned CNTs. The thickness-wise variability function $a({\zeta }_1,{\zeta }_2)$ can be described using Bernstein polynomials over the whole shell with discrete parameters $\mathbf{a}$:

$$a({\zeta }_1,{\zeta }_2)={\mathbf{N}}_a({\zeta }_1,{\zeta }_2)\mathbf{a}.$$

(20)

By simply increasing the polynomial order, more complex material distributions can be obtained. Considering a certain feature $\mathcal{P}\left[\mathbf{a}\right]$ of the shell response, its minimization is expressed as

$$\begin{array}{cccc}\hfill & \underset{\mathbf{a}}{\mathrm{minimise}}\hfill & \hfill & \mathcal{P}\left[\mathbf{a}\right]\hfill \\ \hfill & \mathrm{subject}\mathrm{to}\hfill & \hfill & -1\le {\mathbf{N}}_a({\zeta }_{1g},{\zeta }_{2g})\mathbf{a}\le 1\forall g.\hfill \end{array}$$

(21)

5.2. Optimization of Randomly Orientated CNTs Volume Fraction

The optimal distribution of CNTs in this instance is obtained by changing the volume fraction ${\varphi }_C({\zeta }_1,{\zeta }_2)$ over the mid-surface of the shell, while maintaining it constant through the thickness direction ${\zeta }_3$. For this optimization problem, all parameters describing the domain variation of the CNTs volume fraction ${\varphi }_C({\zeta }_1,{\zeta }_2)$ are collected in vector $\mathit{\varphi }$ and interpolated as

$${\varphi }_C({\zeta }_1,{\zeta }_2)={\mathbf{N}}_{\varphi }({\zeta }_1,{\zeta }_2)\mathit{\varphi }.$$

As in the previous case, the CNTs volume faction distribution at a point of the domain can be obtained using Bernstein polynomials throughout the whole shell. The constitutive matrix $\mathit{L}=\overline{\mathit{L}}$ is that obtained for the randomly oriented CNTs.

The minimization problem is cast in the form

$$\begin{array}{cccc}\hfill & \underset{\mathit{\varphi }}{\mathrm{minimize}}\hfill & \hfill & \mathcal{P}\left(\mathit{\varphi }\right)\hfill \\ \hfill & \mathrm{subject}\mathrm{to}\hfill & \hfill & {\varphi }_C^{min}\le {\mathbf{N}}_{\varphi }({\zeta }_{1g},{\zeta }_{2g})\mathit{\varphi }\le {\varphi }_C^{max}\forall g\hfill \\ \hfill & \hfill & & V_C=h{\int }_{\mathsf{\Omega }}{\varphi }_{C}d\mathsf{\Omega }={\overline{V}}_C\hfill \end{array}$$

(22)

where the point-wise value ${\varphi }_C$ is constrained to vary between the upper bound ${\varphi }_C^{max}$ and the lower bound ${\varphi }_C^{min}$ and the total CNT volume is set to the fixed value ${\overline{V}}_C$.

5.3. Optimization of the in-Plane CNTs Orientation

The CNTs are assumed to be aligned along a direction parallel to the shell mid-surface at an angle $\theta$ with respect to $e_1$. The constitutive matrix at each point of the shell is $\mathit{L}=\mathit{R}\overline{\mathit{L}}{\mathit{R}}^{\top }$ with the rotation matrix $\mathit{R}=\mathit{R}\left(\theta \right)$, with $\theta$ being constant through the thickness, and the constitutive matrix $\overline{\mathit{L}}$ obtained for the aligned case. The spatial distribution of orientation angles $\theta ({\zeta }_1,{\zeta }_2)$ is described using Bernstein polynomials over the whole shell. The angle distribution can be expressed as

$$\theta ({\zeta }_1,{\zeta }_2)={\mathbf{N}}_{\theta }({\zeta }_1,{\zeta }_2)\mathit{\vartheta }.$$

upon collecting the discrete angle parameters in the vector $\mathit{\vartheta }$.

The optimization problem reads

$$\begin{array}{cccc}\hfill & \underset{\mathit{\vartheta }}{\mathrm{minimize}}\hfill & \hfill & \mathcal{P}\left(\mathit{\vartheta }\right)\hfill \\ \hfill & \mathrm{subject}\mathrm{to}\hfill & \hfill & \frac{-\pi }{2}\le {\mathbf{N}}_{\theta }({\zeta }_{1g},{\zeta }_{2g})\mathit{\vartheta }\le \frac{\pi }{2}\forall g\hfill \end{array}$$

(23)

6. Objective Function and Optimization Algorithm

The objective function. The optimization process is aimed at maximizing the collapse load of the nanocomposite shells. In buckling problems, the collapse load can be defined as the lower bound between the critical limit load ${\lambda }_{lim}$ and the load associated with a deformation limit ${\lambda }_{def}$. We denote with $\mathit{\alpha }$ the vector collecting the generic design optimization parameters. In this work, $\mathit{\alpha }$ coincides with the material variables: either $\mathbf{a}$, $\mathit{\varphi }$ or $\mathit{\vartheta }$ introduced in the previous section for the three stated optimization problems, respectively. The objective function can thus be written as

$$\mathcal{P}\left(\mathit{\alpha }\right)=-{\lambda }_c=-\mathrm{m}\mathrm{i}\mathrm{n}\left({\lambda }_{lim},{\lambda }_{def}\right).$$

(24)

The evaluation ${\lambda }_{lim}$ and ${\lambda }_{def}$, and thus the objective function computation, requires the construction of the equilibrium path of the structure for assigned design variables $\mathit{\alpha }$. A common approach for path following the equilibrium curve is the Riks arc-length method [16,35]. In this case, the nonlinear equations in the kinematic unknowns are solved step-by-step using the Newton-Raphson method. However, this kind of analysis bears a significant computational cost due to the large size of the matrices associated with a the high number of DOFs. Furthermore, a reliable evaluation of the equilibrium path should take into account the sensitivity of the structure to imperfections, in particular, to geometrical imperfections. In this case, the nonlinear analysis has to be repeated for a large number of imperfection shapes [26,27], in order to detect the worst imperfection scenario [22]. For this reason, in this work we use an alternative approach called Koiter’s method [25]. Hence, a reduced order model based on Koiter’s theory of elastic stability is assembled for the assigned material configuration. Then, the corresponding reduced nonlinear equations, usually in a lower number of unknowns, are solved to obtain a good estimate of the equilibrium path with the benefit of a much lower computational cost. Moreover, the most interesting feature of the method is the possibility of obtaining the equilibrium path of the imperfect structure by including imperfections a-posteriori in the reduced system [28] of the perfect structure, thus enabling inexpensive sensitivity analyses.

Koiter’s method. The structure is first discretized using the isogeometric environment described in Section 3. Then, the stored energy of each element is rewritten in a mixed form using the stresses at each integration point ${\mathit{\sigma }}_g$ as independent variables [17]

$${\Phi }_e=\sum _g\left({\mathit{\sigma }}_g^T{\mathit{\epsilon }}_g\left({\mathbf{u}}_e\right)-\frac{1}{2}{\mathit{\sigma }}_g^T{\mathbf{C}}_g^{-1}\left(\mathit{\alpha }\right){\mathit{\sigma }}_g\right)w_g$$

(25)

As shown in [36], this is necessary to improve Koiter’s method accuracy, because of the direct prediction of the stress, and efficiency, due to the vanishing of the forth order strain energy variations.

By collecting in vector $\mathbf{z}$ the global discrete displacements $\mathbf{u}$ and stresses ${\mathit{\sigma }}_g$ at each integration point g, Koiter’s method is based on the following reduced order model:

$$\mathbf{z}(\lambda ;{\psi }_i)=\lambda \overline{\mathbf{z}}+\sum _{i=1}^m{\psi }_i{\dot{\mathbf{v}}}_i+\frac{1}{2}\sum _{i,j=1}^m{\psi }_i{\psi }_j{\mathbf{w}}_{ij}+\frac{1}{2}{\lambda }^2\widehat{\widehat{\mathbf{w}}}$$

(26)

where ${\psi }_i$ are the scalar values defined as modal amplitudes, $\widehat{\mathbf{z}}$ is the linear elastic solution (path tangent to the stress-free configuration), ${\dot{\mathbf{v}}}_i$ denotes the ith of the n buckling modes, ${\mathbf{w}}_{ij}$ and $\widehat{\widehat{\mathbf{w}}}$ are quadratic corrections. The evaluation of these vectors requires the solution of linear systems for $\widehat{\mathbf{z}}$, ${\mathbf{w}}_{ij}$ and $\widehat{\widehat{\mathbf{w}}}$ and a linearized buckling analysis for ${\dot{\mathbf{v}}}_i$. Details can be found in [25,36].

According to this choice and after a few asymptotic expansions, it is possible to obtain the following nonlinear algebraic reduced system of equilibrium equations where the unknowns are ${\psi }_i$ and $\lambda$

$$\begin{array}{cc}\hfill r_k[\lambda ,{\psi }_1,\cdots ,{\psi }_m]& ={\mu }_k\left[\lambda \right]+({\lambda }_k-\lambda ){\psi }_k-\frac{1}{2}{\lambda }^2\sum _{i=1}^m{\psi }_i{\mathcal{C}}_{ik}+\frac{1}{2}\sum _{i,j=1}^m{\psi }_i{\psi }_j{\mathcal{A}}_{ijk}\hfill \\ \hfill & +\frac{1}{6}\sum _{i,j,h=1}^m{\psi }_i{\psi }_j{\psi }_h{\mathcal{B}}_{ijhk}=0,k=1\cdots m\hfill \end{array}$$

(27)

The coefficients ${\mathcal{A}}_{ijk}$, ${\mathcal{C}}_{ik}$, ${\mathcal{B}}_{ijhk}$ and ${\mu }_k\left[\lambda \right]$ are scalar quantities evaluated as the sum of elemental contributions of the stored energy variations. Their explicit expressions can be found in [25,36].

A remarkable advantage afforded by Koiter’s method is the possibility of performing computationally efficient and robust imperfection sensitivity analysis of elastic structures. Once Equation (27) has been solved for the perfect structure, the imperfect structure can be studied by only perturbing a posteriori the same system, or rather by adding to it the imperfection term ${\tilde{\mu }}_k$. Therefore, the system in Equation (27) becomes

$$r_k+{\tilde{\mu }}_k=0,k=1\cdots m.$$

(28)

This means that all coefficients of (28) coincide with those evaluated for the perfect structure and thus the analysis of the effect of a new geometrical imperfection simply requires to update ${\tilde{\mu }}_k$ and solve again Equation (28). In this fashion, it is possible to test thousands of imperfections in a few seconds. It is worth noting that, for a given structure, it is impossible to obtain the same amount of results in such a short time if the nonlinear analysis is performed with the full discrete model.

For the evaluation of ${\tilde{\mu }}_k$, two strategies were previously proposed [28]. A first solution is very fast but its range of validity is limited to small imperfection amplitudes and almost linear pre-buckling path. A second strategy was recently proposed by [28] and, by overcoming many drawbacks of the first formulation, it is recommended when the effects of the imperfection are more important. In this work, the geometrical imperfection $\tilde{\mathbf{u}}$ is assumed to be a linear combination of modal displacement shapes ${\dot{\mathbf{u}}}_i$ with combination factors $\widetilde{{\psi }_i}$,

$$\tilde{\mathbf{u}}=\sum _{i=1}^m\widetilde{{\psi }_i}{\dot{\mathbf{u}}}_i,\left|\right|\tilde{\mathbf{u}}\left|\right|\le {\tilde{u}}_{max},$$

(29)

scaled in order to have an assigned maximum amplitude ${\tilde{u}}_{max}$ chosen, for example, from experimental measurements.

The Global Convergent Method of Moving Asymptotes. The optimal design problem is solved using a gradient-based optimizer, i.e., the Global Convergent Method of Moving Asymptotes (GCMMA) [33,38,39,40]. This algorithm devised for the optimization of objective functions requires a relatively high computational cost to evaluate the gradient and is characterized by many optimization variables. It is based on convex subsequent approximations of the objective function.

Constructing the gradient of the objective function with respect to the design parameters is not a simple task and is thus computed numerically. The i-th component of the gradient is evaluated according to a forward finite difference scheme as

$$\nabla {\mathcal{P}}_i\left(\mathit{\alpha }\right)\approx \frac{P(\mathit{\alpha }+d{\mathbf{e}}_i)-P\left(\mathit{\alpha }\right)}{d},$$

(30)

where d is a conveniently small finite difference and ${\mathbf{e}}_i$ is a basis vector whose ith component is 1 while the others are zero.

Although the task of evaluating the objective function is extremely time-consuming when using the full discrete model, the efficiency of Koiter’s method allows its computational cost to be reduced to within an acceptable level. Moreover, since GCMMA generally converges with relatively few iterations, the number of gradient evaluations is small and the overall computational cost of the optimization is sustainable.

7. Numerical Results

In this section, the CNTs distribution is optimized in order to improve the nonlinear response under buckling. The three types of optimizations discussed in Section 3 are here summarized:

-

OPT1: through-the-thickness optimization, i.e., variable through-the-thickness distribution of aligned CNTs across the shell mid-surface with assigned mid-surface value;

-

OPT2: optimization across the mid-surface, i.e., the volume fraction of randomly oriented CNTs is kept constant through the thickness but can vary within the shell mid-surface with a constraint on the overall CNTs volume;

-

OPT3: optimization of the orientation, i.e., the CNTs volume is assigned, the orientation of the CNTs can vary across the shell mid-surface but is constant through the thickness.

Three applications are considered. The first samples are an Euler beam and a simply supported plate in compression exhibiting a stable post-buckling behavior. In this case, the optimizations are expected to provide a higher buckling load and a less compliant post-buckling response. The last test regards a curved panel, whose post-buckling behavior can be either stable or unstable depending on the material distribution. This is a more significant test to show the formidable potential of a variable CNTs distribution in the design of thin nanocomposite structures.

We consider shells made of an isotropic polymeric matrix with Young modulus $E=2.47$ GPa and Poisson’s ratio $\nu =0.36$. CNTs are considered as an isotropic elliptical inclusion with aspect ratio equal to 732, Young modulus $E=970$ GPa and Poisson’s ratio $\nu =0.1$.

Two material descriptions over the mid-surface are considered: polynomials of order 4 and 9. Due to the problem symmetry, only the control points of a quarter of the structure are considered as independent variables, leading to $3\times 3$ and $5\times 5$ design control points, respectively, for the two descriptions. The number of optimization variables is 9 and 25 for the two cases, respectively.

For the OPT2 problem, the volume fraction bounds are set to ${\varphi }_C^{min}=0.1\%$ and ${\varphi }_C^{max}=15\%$ while the CNTs volume is given by $V_C={\varphi }_C^{*}V$.

It is worth noting that the effective CNTs volume fraction decreases by increasing ${\varphi }_C$ according to [14] as shown in

Table 1. Effective volume fraction (in %).

nominal	0.50	0.75	1.00	2.00	5.00	10.0	15.0
effective	0.48	0.72	0.96	1.89	4.45	7.570	8.50

due to CNTs agglomeration phenomena which are detrimental for the load transfer.

7.1. Nanocomposite Beam Under Compression

This test consists of a simply supported Euler beam under compression on two opposite ends as shown in Figure 2.

This simple test is considered to show the correctness of the proposed numerical approach. A single linearized buckling mode is used to construct the ROM of the Koiter method and the geometrical imperfection. The imperfection shape is scaled such that the maximum deviation ${\tilde{u}}_{max}$ is $0.1$ of the shell thickness h. The deformation limit is set to $v=2$ mm with v the mean value of the flexural displacement.

Optimization of the CNT volume fraction for random orientations. The optimization problem OPT2 only is considered because OPT1 and OPT3 do not yield any improvement for this specific problem. The results of OPT2 are compared with the performance of a uniform distribution of ${\varphi }_C$ (hereafter referred to as UD). Two polynomial orders are used to describe the distribution of ${\varphi }_C$ across the surface: order 4 and order 9. First, we compare the linearized buckling loads reported in

Table 2. Euler beam: linearized buckling loads for the optimal volume fraction distribution ${\varphi }_C({\zeta }_1,{\zeta }_2)$ with assigned average value ${\varphi }_C^{*}$ described by Bernstein polynomials of order 4 normalized with respect to ${\lambda }_r=0.0005313$ KN/mm.

Mode	$0.50\%$		$0.75\%$		$1.00\%$		$2.00\%$		$5.00\%$		$10.0\%$
	OPT	UD	OPT	UD	OPT	UD	OPT	UD	OPT	UD	OPT	UD
1	0.5899	0.5386	0.6951	0.6163	0.7977	0.6936	1.1867	1.0000	2.2341	1.8662	3.1037	2.9634
2	2.1437	2.1621	2.4200	2.4737	2.6829	2.7843	3.7129	4.0141	7.3660	7.4911	12.027	11.894
3	4.8159	4.8906	5.4150	5.5955	5.9816	6.2982	8.1653	9.0798	15.715	16.944	26.985	26.904
4	8.6074	8.7516	9.6712	10.013	10.675	11.270	14.512	16.248	27.450	30.322	48.143	48.142

for order 4. We note an increase of the first load for all average volume fractions ${\varphi }_C^{*}$ and in particular for the case ${\varphi }_C^{*}=5\%$ where the improvement reaches $20\%$. Similar considerations hold for the description with order 9 as shown in

Table 3. Euler beam: linearized buckling loads for the optimal volume fraction distribution ${\varphi }_C({\zeta }_1,{\zeta }_2)$ with assigned average value ${\varphi }_C^{*}$ described by Bernstein polynomials of order 9 normalized with respect to ${\lambda }_r=0.0005313$ KN/mm.

Mode	$0.50\%$		$0.75\%$		$1.00\%$		$2.00\%$		$5.00\%$		$10.0\%$
	OPT	UD	OPT	UD	OPT	UD	OPT	UD	OPT	UD	OPT	UD
1	0.5988	0.5386	0.7020	0.6163	0.8026	0.6936	1.1897	1.0000	2.2694	1.8662	3.2106	2.9634
2	2.0477	2.1621	2.3072	2.4737	2.5889	2.7843	3.6998	4.0141	6.8262	7.4911	12.350	11.894
3	4.7410	4.8906	5.2998	5.5955	5.8482	6.2982	8.0605	9.0798	14.120	16.944	27.260	26.904
4	8.4905	8.7516	9.4941	10.013	10.486	11.270	14.363	16.248	24.777	30.322	48.199	48.142

. A similar improvement is highlighted by a full nonlinear analysis in Figure 3. Finally, the optimal distribution of CNTs volume fraction is depicted in Figure 4 and Figure 5 for the two orders used to describe ${\varphi }_C$ over the surface. As expected, it is possible to observe that, by comparison with the initial uniform distribution, the CNTs volume fraction is greater at the midspan of the beam and lower at the end sections in order to maximize the flexural stiffness. Moreover, it is worth noting that the effectiveness ratio of the CNTs decreases by increasing ${\varphi }_C$ as shown in Table 1. This is the reason why the optimal distribution tends to be more uniform and equal to the maximum admissible fraction near the midspan as shown in Figure 4 and Figure 5 for high values of ${\varphi }_C^{*}$.

7.2. Nanocomposite Plate under Compression

This test deals with a simply supported square plate under compression on two opposite sides as shown in Figure 6.

The thickness was set to 5 mm while the span is 508 mm, and the uniform compression load was $q=0.01$ KN/mm.

A single linearized buckling mode is used to construct the ROM of the Koiter method and the geometrical imperfection. The imperfection shape is scaled such that the maximum deviation ${\tilde{u}}_{max}$ is $0.1h$, with h being the shell thickness. The deformation limit is set to $v=2$ mm with v being the mean value of the out-of-plane displacement. Due to the symmetry of the test, OPT1 does not provide any improvement compared to a uniform volume fraction distribution and is omitted. Instead, the focus is on OPT2 and OPT3, respectively.

Optimization of the CNTs volume fraction for random orientations. The optimization problem OPT2 is discussed here, making use of a comparison with the solution obtained for the uniform distribution of ${\varphi }_C$ (UD). The linearized buckling loads are reported in

Table 4. Square plate: linearized buckling loads for the optimal volume fraction distribution ${\varphi }_C({\zeta }_1,{\zeta }_2)$ with assigned average value ${\varphi }_C^{*}$ described by Bernstein polynomials of order 4.

Mode	$0.50\%$		$0.75\%$		$1.00\%$		$2.00\%$		$5.00\%$		$10.0\%$
	OPT	UD	OPT	UD	OPT	UD	OPT	UD	OPT	UD	OPT	UD
1	0.6736	0.6358	0.7835	0.7276	0.8892	0.8191	1.2947	1.1810	2.4128	2.2030	3.5968	3.4943

for the polynomial description of ${\varphi }_C({\zeta }_1,{\zeta }_2)$ of order 4. We note a general but very slight increase of the first load for all average volume fractions ${\varphi }_C^{*}$. Similar considerations hold for the parametrization of order 9 where the optimized solution is slightly better. The same considerations hold if we consider the full nonlinear paths shown in Figure 7. Finally, the optimal distribution of CNTs volume fraction is depicted in Figure 8 for 9th-order polynomials employed to describe ${\varphi }_C$ across the surface. For the simply supported square plate, the use of a nonuniform distribution of randomly CNTs within the mid-surface does not yield significant improvements.

Optimization of the CNTs orientation within the mid-surface. The optimization problem OPT3 is investigated next. It consists in optimizing the orientation of the CNTs in the plate plane, while keeping the volume fraction constant. The results of OPT3 are compared with those corresponding to all CNTs uniformly aligned with the load direction (UD). The linearized buckling loads are reported in

Table 5. Square nanocomposite plate: linearized buckling loads for the optimal CNTs orientation $\theta ({\zeta }_1,{\zeta }_2)$ described by Bernstein polynomials of order 4 for assigned volume fraction ${\varphi }_C$.

Mode	$0.50\%$		$0.75\%$		$1.00\%$		$2.00\%$		$5.00\%$		$10.0\%$
	OPT	UD	OPT	UD	OPT	UD	OPT	UD	OPT	UD	OPT	UD
1	0.8611	0.6907	1.0369	0.8038	1.2020	0.9131	1.7696	1.3199	3.1007	2.3369	4.5583	3.4909

for the polynomial parametrization of the orientation of order 4. It is possible to observe that the first buckling load is increased notably by the optimization and the improvement gets better with the volume fractions. Figure 9 with the full nonlinear path shows that a higher post-buckling stiffness can be obtained for the optimized solutions compared to the uniform orientation and an even better behavior is obtained with the parametrization of order 9. The results are completed with Figure 10 depicting the optimal orientation distribution of order 9.

7.3. Cylindrical Nanocomposite Panel under Compression

The last test in Figure 11 features a cylindrical panel under compression. The structure presents a simple geometry but complex post-buckling behavior, involving multi-modal interactions with possible unstable path and significant imperfection sensitivity. This kind of panel was already subject to optimization for composite laminates by means of a stochastic optimizer [26,27]. Additionally, it was taken as a benchmark example to test the accuracy of the isogeometric solid-shell model [36].

The optimizations seek the material distributions that maximize the collapse load considering as deformation limit the axial displacement of the loaded section $v_{ls}=0.2$ mm. The maximum amplitude of the geometrical imperfection is ${\tilde{u}}_x=0.5h$. The first eight linearized buckling modes are used to construct the ROM of the Koiter method and the geometrical imperfection shape.

Through-the-thickness optimization of aligned CNTs volume fraction. We start with the through-the-thickness optimization problem referred to as OPT1. The CNTs are aligned along the load direction. The results are compared with those obtained for a uniform through-the-thickness distribution (UD). First,

Table 6. Cylindrical nanocomposite panel: linearized buckling loads for the optimal through-the-thickness CNTs distribution with function $a({\zeta }_1,{\zeta }_2)$ in Equation (20) described by Bernstein polynomials of order 4 and assigned volume fraction ${\varphi }_C^{*}$ normalized with respect to ${\lambda }_r=0.006466$ KN/mm.

Mode	$0.50\%$		$0.75\%$		$1.00\%$		$2.00\%$		$5.00\%$		$10.0\%$
	OPT	UD	OPT	UD	OPT	UD	OPT	UD	OPT	UD	OPT	UD
1	0.6325	0.6065	0.7230	0.6878	0.7920	0.7605	1.0336	1.0000	1.4586	1.4978	1.9950	1.9873
2	0.7947	0.7861	0.9166	0.9108	1.0286	1.0283	1.4351	1.4507	2.1035	2.4574	3.3456	3.5521
3	1.1102	1.0811	1.2726	1.2461	1.4225	1.4029	1.9248	1.7974	2.8299	2.6611	3.7066	3.5924
4	1.3092	1.2268	1.4444	1.3339	1.5427	1.4354	1.9806	1.9787	2.9550	3.3660	4.6859	4.8864
5	1.3720	1.3669	1.6309	1.6260	1.8700	1.8737	2.7485	2.7820	4.2409	5.0083	6.9759	7.4824
6	1.5385	1.5380	1.8357	1.8323	2.1058	2.1160	3.1257	3.1646	4.8939	5.5607	7.6127	7.7939
7	2.2954	2.2719	2.5844	2.5079	2.8048	2.7305	3.6036	3.5502	5.1469	5.7560	8.2031	8.6383
8	2.3122	2.3043	2.6621	2.6583	2.9917	3.0021	4.2227	4.2762	6.3765	7.3119	10.0125	9.1880

shows that the first linearized buckling load turns out to be almost unaffected by the optimization process. On the contrary, the full nonlinear analysis reported in Figure 12 shows how the optimization globally turns the behavior from unstable into stable using the same CNTs volume. The snap-through behavior is completely suppressed at the cost of a very slight stiffness reduction in the pre-critical range. Similar results in terms of equilibrium paths are obtained with a polynomial description of the function $a({\zeta }_1,{\zeta }_2)$ of order 9. The shape of $a({\zeta }_1,{\zeta }_2)$ is very similar for the two orders and is reported in Figure 13 for order 4.

Optimization of CNTs volume fraction for random orientations. The optimization problem OPT2 is considered next. The results are compared with the performance of a uniform distribution of randomly oriented CNTs (UD). First, we compare the linearized buckling loads reported in

Table 7. Cylindrical nanocomposite panel: linearized buckling loads for the optimal volume fraction distribution ${\varphi }_C({\zeta }_1,{\zeta }_2)$ with assigned average value ${\varphi }_C^{*}$ described by Bernstein polynomials of order 4 normalized with respect to ${\lambda }_r=0.006466$ KN/mm.

Mode	$0.50\%$		$0.75\%$		$1.00\%$		$2.00\%$		$5.00\%$		$10.0\%$
	OPT	UD	OPT	UD	OPT	UD	OPT	UD	OPT	UD	OPT	UD
1	0.6075	0.5552	0.6942	0.6353	0.7832	0.7151	1.1173	1.0310	1.9557	1.9237	3.1329	3.0530
2	0.7684	0.7010	0.8776	0.8022	0.9892	0.9030	1.4251	1.3019	2.4889	2.4290	3.9793	3.8548
3	1.0120	1.0087	1.1678	1.1543	1.3151	1.2994	1.8771	1.8736	3.3529	3.4952	5.3862	5.5450
4	1.1833	1.1107	1.3590	1.2708	1.5277	1.4305	2.2466	2.0625	3.9364	3.8483	6.2432	6.1075
5	1.2680	1.2401	1.4611	1.4190	1.6476	1.5974	2.3504	2.3032	4.1880	4.2968	6.7559	6.8168
6	1.3703	1.4200	1.5833	1.6249	1.7772	1.8292	2.5711	2.6374	4.6157	4.9201	7.3073	7.8052
7	1.9120	1.7370	2.1885	1.9875	2.4664	2.2372	3.5623	3.2256	6.2360	6.0183	9.9263	9.5510
8	2.0283	1.9699	2.3364	2.2541	2.6371	2.5374	3.7589	3.6584	6.6961	6.8254	10.807	10.83039

for order 4. We note a general increase of the first load for all average volume fractions ${\varphi }_C^{*}$ and, in particular, for the intermediate ${\varphi }_C^{*}$ whose improvement is about $10\%$. Similar considerations hold for the parametrization of order 9 as shown in

Table 8. Cylindrical nanocomposite panel: linearized buckling loads for the optimal volume fraction distribution ${\varphi }_C({\zeta }_1,{\zeta }_2)$ with assigned average value ${\varphi }_C^{*}$ described by Bernstein polynomials of order 9 normalized with respect to ${\lambda }_r=0.006466$ KN/mm.

Mode	$0.50\%$		$0.75\%$		$1.00\%$		$2.00\%$		$5.00\%$		$10.0\%$
	OPT	UD	OPT	UD	OPT	UD	OPT	UD	OPT	UD	OPT	UD
1	0.6237	0.5552	0.7167	0.6353	0.8029	0.7151	1.1353	1.0310	2.0376	1.9237	3.1647	3.0530
2	0.7799	0.7010	0.8936	0.8022	1.0059	0.9030	1.4267	1.3019	2.5885	2.4290	4.0217	3.8548
3	1.0178	1.0087	1.1699	1.1543	1.3147	1.2994	1.8990	1.8736	3.4311	3.4952	5.4352	5.5450
4	1.1623	1.1107	1.3317	1.2708	1.5204	1.4305	2.1848	2.0625	4.0413	3.8483	6.3539	6.1075
5	1.2778	1.2401	1.4679	1.4190	1.6554	1.5974	2.3892	2.3032	4.3172	4.2968	6.8085	6.8168
6	1.3526	1.4200	1.5470	1.6249	1.7422	1.8292	2.5344	2.6374	4.6328	4.9201	7.4096	7.8052
7	1.9405	1.7370	2.2259	1.9875	2.5083	2.2372	3.5554	3.2256	6.4522	6.0183	10.0560	9.5510
8	2.0844	1.9699	2.3961	2.2541	2.6806	2.5374	3.8509	3.6584	6.9260	6.8254	10.8735	10.8303

where the optimized solution is also slightly improved. However, the great benefit of the variable volume fraction distribution is highlighted by a full nonlinear analysis in Figure 14. The unstable behavior of the uniform distribution is made stable by the optimal CNTs distribution using the same overall amount of CNTs. The slight stiffness reduction in the pre-critical range is compensated by the complete elimination of snap-through, at least in the range of interest. Similar results in terms of equilibrium paths are obtained with a polynomial description of ${\varphi }_C({\zeta }_1,{\zeta }_2)$ of order 9. Indeed, the analysis leads to a very similar optimal distribution for the two orders, reported in Figure 15 for order 9.

Optimization of CNTs orientation. The optimization problem OPT3 is finally discussed here. It consists of optimizing the CNTs orientation within the shell mid-surface, while keeping constant the volume fraction at each point of the structure. The results of OPT3 are compared with those obtained for uniformly aligned CNTs collinear with the load direction (UD). The linearized buckling loads are reported in

Table 9. Cylindrical nanocomposite panel: linearized buckling loads for the optimal CNTs orientation $\theta ({\zeta }_1,{\zeta }_2)$ described by Bernstein polynomials of order 4 for assigned volume fraction ${\varphi }_C$ normalized with respect to ${\lambda }_r=0.006466$ KN/mm.

Mode	$0.50\%$		$0.75\%$		$1.00\%$		$2.00\%$		$5.00\%$		$10.0\%$
	OPT	UD	OPT	UD	OPT	UD	OPT	UD	OPT	UD	OPT	UD
1	0.5956	0.6065	0.6710	0.6878	0.7368	0.7605	0.9519	1.0000	1.6582	1.4978	2.2547	1.9873
2	0.7753	0.7861	0.8947	0.9108	1.0060	1.0283	1.4116	1.4507	2.5809	2.4574	3.4228	3.5521
3	1.0966	1.0811	1.2574	1.2461	1.4100	1.4029	1.7751	1.7974	2.8041	2.6611	3.6373	3.5924
4	1.2540	1.2268	1.3588	1.3339	1.4587	1.4354	1.9597	1.9787	3.5398	3.3660	4.7805	4.8864
5	1.3938	1.3669	1.6453	1.6260	1.8821	1.8737	2.7433	2.7820	4.8703	5.0083	6.7610	7.4824
6	1.5363	1.5380	1.8221	1.8323	2.0941	2.1160	3.0839	3.1646	5.5986	5.5607	7.5653	7.7939
7	2.2102	2.2719	2.5331	2.5079	2.7802	2.7305	3.6117	3.5502	5.8414	5.7560	7.8166	8.6383
8	2.3373	2.3043	2.6906	2.6583	3.0311	3.0021	4.2693	4.2762	7.3583	7.3119	9.7666	9.1880

for the parametrized orientation distribution of order 4. It is possible to observe that the first buckling load becomes notably larger for high volume fractions, while it remains almost the same for low CNT contents. However, looking at Figure 16, the radical change of mechanical behavior is appreciable also for low volume fractions. The optimal CNTs distribution leads to the suppression of the snap-through instability, at the cost of initial stiffness reduction for low ${\varphi }_C$, reduction that becomes negligible for higher volume fractions. Similar results are obtained with the parametrization of order 9. The results are completed with Figure 17 depicting the optimal orientation distribution of order 9.

8. Conclusions

This work tackled the new challenging problem of finding optimal distributions of the nanostructural reinforcement phase (here, carbon nanotubes) within polymeric composite panels subject to buckling and snap-through.

• A numerical strategy for the optimization of the buckling and postbuckling response of nanocomposite shells with variable CNTs distribution was proposed and investigated. The method is based on an integrated isogeometric framework that employs NURBS functions to describe the geometry and displacements while the optimization variables deal with the CNT distributions within the polymer hosting matrix.
• Various CNTs distributions were investigated either through the thickness or within the mid-surface for both aligned CNTs, aligned but varying within the surface or randomly oriented. The obtained through-the thickness distributions can be practically realized in multilayer nanocomposite structures since a continuous law can be reasonably approximated by piece-wise functions when the multilayers are considered sufficiently thin.
• The outcomes of extensive numerical tests have proved that the limit load can be largely improved for optimal CNTs distributions in the sense of strategically deploying the nanofibers where the maximum elastic stiffness fighting against the negative stiffness can be attained.
• Most importantly, it has been shown that shallow shells, which are dangerously prone to snap-through, can become globally stable if the CNTs are optimally distributed. This is a remarkable result on the global stability of nanocomposite shallow shells which, when properly designed, do not show any snap-through and thus can be safely employed in engineering applications. Mention must be made of the fact that these nanocomposite panels also show the additional advantage of exhibiting enhanced damping capability thanks to the CNT/polymer interfacial dissipation which makes these structures generally more stable against dynamic loads.
• This work has shown the potential of optimizing nonlinear structural behaviors using the unprecedented flexibility afforded by the CNTs nanoreinforcement which not only acts to shift the elastic loss of stability towards higher stresses, but can also either suppress snap-through or make the response less compliant in the postbuckling range. The next step of the research will be the optimal design of high-performance and lightweight vehicles, aerostructures and devices.