A cost-effective isogeometric approach for composite plates based on a stress recovery procedure

TL;DR

This paper presents a cost-effective isogeometric method for simulating laminated composite plates, utilizing stress recovery from coarse solutions to accurately predict out-of-plane stresses, especially effective for slender multilayer stacks.

Contribution

It introduces a novel stress recovery procedure within an isogeometric framework that reduces computational costs while maintaining accuracy for composite plate analysis.

Findings

Effective stress recovery from coarse solutions

High accuracy in out-of-plane stress prediction

Significant computational savings for multilayer stacks

Abstract

This paper introduces a cost-effective strategy to simulate the behavior of laminated plates by means of isogeometric 3D solid elements. Exploiting the high continuity of spline functions and their properties, a proper out-of-plane stress state is recovered from a coarse displacement solution using a post-processing step based on the enforcement of equilibrium in strong form. Appealing results are obtained and the method is shown to be particularly Peffective on slender composite stacks with a large number of layers.

Equations

$$\mathbf{\sigma }=\mathbb{C}\mathbf{\epsilon }$$

$$\mathbb{C}={\left[\begin{array}{cccccc}\frac{1}{E_1}& -\frac{{\nu }_{12}}{E_2}& -\frac{{\nu }_{13}}{E_1}& 0& 0& 0\\ -\frac{{\nu }_{21}}{E_2}& \frac{1}{E_2}& -\frac{{\nu }_{23}}{E_2}& 0& 0& 0\\ -\frac{{\nu }_{31}}{E_3}& -\frac{{\nu }_{32}}{E_3}& \frac{1}{E_3}& 0& 0& 0\\ 0& 0& 0& \frac{1}{G_{23}}& 0& 0\\ 0& 0& 0& 0& \frac{1}{G_{31}}& 0\\ 0& 0& 0& 0& 0& \frac{1}{G_{12}}\end{array}\right]}^{-1},$$

$$E_2$$

$$G_{12}$$

$${\nu }_{12}$$

$$p(x,y)={\sigma }_0\sin \left(\frac{\pi x}{St}\right)\sin \left(\frac{\pi y}{St}\right)$$

$${\bar{\sigma }}_{ij}$$

$${\bar{\sigma }}_{i3}$$

$${\bar{\sigma }}_{33}$$

$$\nabla \cdot \mathbf{\sigma }=\mathbf{b}$$

$${\sigma }_{11,1}+{\sigma }_{12,2}+{\sigma }_{13,3}$$

$${\sigma }_{12,1}+{\sigma }_{22,2}+{\sigma }_{23,3}$$

$${\sigma }_{13,1}+{\sigma }_{23,2}+{\sigma }_{33,3}$$

$${\sigma }_{ij,k}=\frac{\partial {\sigma }_{ij}}{\partial x_k}.$$

$${\sigma }_{13}(z)$$

$${\sigma }_{23}(z)$$

$${\sigma }_{33}(z)$$

$${\sigma }_{33}(z)={\int }_{z_0}^z({\sigma }_{11,11}+{\sigma }_{22,22}+2{\sigma }_{12,12}+b_3-b_{1,1}-b_{2,2})\mathrm{d}z+{\sigma }_{33}(z_0).$$

$${\sigma }_{ij,k}$$

$${\sigma }_{ij,kl}$$

$${\epsilon }_{mn,k}$$

$${\epsilon }_{mn,kl}$$

$$\text{error}=\frac{\max ({\sigma }_{\text{analytic}}-{\sigma }_{\text{recov}})}{\max ({\sigma }_{\text{analytic}})}.$$

Full text

1 Introduction

Composite materials are used in a wide variety of fields such as aerospace or automotive. The study of composite laminates has become more and more important along the past years especially because of their light weight and very resistant mechanical properties. A composite laminate is usually made of several layers of highly resistant fibers embedded in soft matrix. Laminate structures tend to be prone to damage at the interfaces between layers, this mode of failure being referred to as delamination. The prediction and evaluation of damage in composite laminates demands an accurate evaluation of the three-dimensional stress state through the thickness, although most of the studies available in the literature consider the laminate as a two-dimensional object. Classical two-dimensional theories such as shell approaches are not accurate enough to reliably predict interlaminar damage and delamination. Accordingly, a number of layerwise theories have been developed to compute more accurately the mechanical state inside the laminate. Such methods often rely on heavy computations and hybrid approaches in order to be able to capture the complex behaviour of the interlaminar interfaces. The numerical counterpart of such layerwise or hybrid theories was approached mostly using standard finite element discretization (see, e.g.,  [1, 2] and references therein).

Over the last decade, many novel methods have been proposed to enhance the standard finite element framework. Among them, Isogeometric Analysis (IGA) is a concept proposed by Hughes et al. [3] where the shape functions used in computer aided design are also approximating physical fields and state variables. It thus relies on spline functions, like, e.g., NURBS (Non Uniform Rational B-Splines), to approximate both the geometry and the solution. In addition to improve the pipeline between design and analysis, spline shape functions possess properties which can be used to drastically improve the performance of numerical analysis compared to standard finite element method. For example, spline shape functions can provide easily high order approximation and highly simplify the refinement process. Moreover, their smoothness guarantees higher accuracy and opens the door to the discretization of high-order PDEs in primal form. IGA has been successfully used to tackle a large variety of problems, including solid and structures ( see, e.g.,  [4, 5, 6, 7, 8]), fluids (see, e.g.,  [9, 10]), fluid-structure interaction (see, e.g.,  [11, 12]), and other fields (see e.g.,  [13, 14]).

IGA methods have already been used to solve composite laminate problems. A wide variety of laminated models have been implemented and studied within the IGA framework, especially relying on high-order theories [15, 16, 17]. These techniques mostly utilize enhanced shell and plate theories. In addition to these 2D approaches, some methods can be used to compute the full 3D stress state using 3D isogeometric analysis such as in [17, 18, 19]. In such approaches, each ply of the laminate is delimited by a C0-continuous interface. Such an approach (referred in the following as “layerwise” ) involves a large number of degrees of freedom when a lot of layers are composing the laminate. Thus it may be not completely satisfactory, despite its accuracy, due to its high cost.

This paper presents an approach consisting in using 3D computations with a reduced number (namely, one) of elements through the thickness and a layerwise integration rule, which is then post-processed in order to obtain an accurate 3D stress state. The post-processing relies on the integration of the equilibrium equations (such as in the recovery techniques proposed in [16, 20, 21, 22, 23, 24, 25]) to compute the stresses through the thickness from the in-plane ones. This allows to drastically reduce the computational time compared to the 3D layerwise approach as the number of degrees of freedom is significantly decreased. The post-processing operation is in fact very fast and its cost does not increase significantly with the number of degrees of freedom. Solutions obtained using the proposed technique are actually close to those provided by full “layerwise” 3D isogeometric analysis, at a fraction of the cost, in particular when many layers are present.

The effectiveness of the proposed approach relies on two ingredients, both granted by the peculiar properties of IGA, namely:

1. (i)

The capability of obtaining accurate in-plane results with a coarse mesh with only one element through the thickness; 2. (ii)

The higher continuity allowing an accurate computation of stresses and of their derivatives from displacements.

The structure of the paper is as follows. First, the considered isogeometric strategies are briefly introduced. Then, the 3D geometry and the mechanical problem used as a test case is detailed along with its analytical solution. Results using the two considered standard isogeometric approaches are then presented. Finally, the post-processing approach is proposed along with numerical results and a comparison with those provided by full “layerwise” 3D isogeometric analysis. Several cases using different numbers of quadrature points per material layer, numbers of layers and thicknesses are considered to show the effectiveness of the method.

2 Standard IGA strategies for 3D analysis of laminates

In this section we present a preliminary discussion on standard IGA strategies for 3D analysis of laminates. For the sake of conciseness we avoid reporting an introduction on the basic concepts of IGA for 3D elasticity, including basics of B-Splines, NURBS, etc. Readers are instead refered to [26] and references herein. Here, we just recall that we will consider a standard IGA consisting of an isoparametric Galerkin formulation where splines are used for approximating both geometry and displacements. As opposed to the framework typically used to compute laminate shell and plate solutions [1, 16, 27], we choose here to perform a full “layerwise” 3D computation. Such an approach has already been used in combination with IGA in order to solve multilayered plate problems [17, 18, 19] and along those lines we start considering two different sets of shape functions, as follows:

Layerwise approach: The first considered approach is similar to high order finite element methods. Each layer is in fact modeled by one patch through the thickness and $C^{0}$ continuity is kept between elements (see Figure 1(a)). Such an approach is completely layerwise and, as a consequence, the number of degrees of freedom is proportional to the number of layers. A standard integration rule is adopted, namely $p$ + 1 points in each direction ($p$ being the degree of the shape functions). This approach is equivalent to the one proposed in [17, 19].

Single-element approach: The second approach actually uses a single element through the thickness (see Figure 1(b)), strongly reducing the number of degrees of freedom with respect to the previously mentioned method. To account for the presence of the layers, a special integration rule is adopted, consisting of a $q$-point Gauss rule ($q\geq 1$) over each layer. This can be considered as a special homogenized approach with a layerwise integration. Such a method would a priori not give sufficiently good results (in terms of through-the-thickness stress description) but could be easily coupled with a post-processor to improve the solution, and this will be the object of Section 4.

It should be noted that the in-plane continuity of the shape functions in both approaches is the same.

3 Test case: The Pagano layered plate

The test case used in this study is the classical one proposed by Pagano [28]. It consists of a simply supported multilayered 3D plate with a sinusoidal loading on top and a loading-free bottom face (see Figure 2). This problem can be easily parameterized which allows to analyze many cases in terms of layer number and distribution. In the following study, a few examples are considered using different numbers of layers (i.e., 3, 4, 11, and 34). In all these cases, the loading conditions are the same (namely, a two dimensional sinus with a period equal to twice the length of the plate), while the thickness of every single layer is set to $1\,\text{mm}$, and the length of the plate is chosen to be $S$ times larger than the total thickness $t$ of the laminate. Figure 2 shows the summary of the test problem in the case of three layers.

The laminate is composed by orthotropic layers placed orthogonally on top of each other (thus creating a 90/0/90/… laminate). For each layer, we have:

[TABLE]

where the elasticity tensor $\mathbb{C}$ can be expressed, using Voigt notation, as:

[TABLE]

and the material parameters for these layers are the following:

[TABLE]

The pressure field used as the top boundary condition is

[TABLE]

In our test, we always select $E_{1}=25$ GPa, $G_{23}=0.5$ GPa and $\sigma_{0}=1$.

3.1 Analytical results

Each stress component presents a different distribution along thickness. In particular, in-plane normal components of the stress tensor (i.e., $\sigma_{11}$ and $\sigma_{22}$) are discontinuous along the thickness although the displacement solution is continuous in all cases. This is obviously due to the differences in term of material parameters between the layers. On the other hand, equilibrium implies that through-the-thickness stress components have to be C0-continuous (i.e., $\sigma_{13}$, $\sigma_{23}$ and $\sigma_{33}$) . Along the article, stress components are normalized, according to [28], using the following formulas:

[TABLE]

Figure 3 reports the analytical solution (computed using the equations given in [28]) for the $\sigma_{13}$, $\sigma_{11}$ and $\sigma_{12}$ components of the stress tensor in the case of 3 and 11 layers (blue solid lines). We can clearly identify the different behavior of the layers, with a “soft” and a “hard” direction.

3.2 First results

In this section, we comment the results obtained using the layerwise and the single-element approaches, which are compared with the analytical solutions for $\sigma_{13}$, $\sigma_{11}$ and $\sigma_{12}$ in Figure 3 for the cases of 3 and 11 layers at the position $X=0.25L$ and $Y=0.25L$. Results are given for an in-plane mesh composed of $9\times 9$ quartic elements, with cubic shape functions through the thickness. All the results shown in this work have been obtained using an in-house code based on the igatools library (see [29] for further details).

The layerwise approach captures accurately the solution for both displacements and stresses, no matter if they are discontinuous or not (see Figure 3). This is due to the layerwise displacement $C^{0}$-continuity which allows a proper description of the stress state even between two layers.

On the contrary, using a poor description (i.e., a single high order element through the thickness) does not allow a correct description of the stress state through the thickness (see Figure 3). The high order continuity of the shape functions leads in fact to the continuity of strains and the jump in the material parameters from one layer to another leads to the discontinuity of the stress components (even when they should be continuous through the thickness). Therefore, the differences between the analytical solution and the computed results are very significant.

Figure 4 gathers the relative computation time for these two approaches with different numbers of layers. As expected, the computational time of the layerwise approach is higher because the number of degrees of freedom is higher (1,183, 5,746 and 17,407 versus 676) as well as the number of integration points required per layer. Note that the computational times of the layerwise approach are, in this case, about 2.5 times higher than the single-element approach with 2 Gauss points per layer.

The layerwise approach seems to be suitable for small stacks (i.e., with a small amount of layers), while it is not suitable for high number of layers as it becomes extremely time and memory consuming.

On the other hand, the single-element approach seems to be more scalable, as only one element through-the-thickness is used, but only the global solution in terms of displacement and in-plane stresses is accurately captured by the simulation. The idea of the paper is to obtain a better approximation of the out-of-plane stresses by post-processing the accurately computed displacements and in-plane stresses, without significantly increasing the overall simulation cost.

4 Reconstruction from Equilibrium

Although in-plane stress components are almost correctly captured by the single-element approach (as it can be seen in Figures 3(c) and 3(d) for the component $\sigma_{11}$ and 3(e) and 3(f) for the component $\sigma_{12}$), the out-of-plane components are not. As interlaminar delamination and other fracture processes rely mostly on such out-of-plane components, a proper through-the-thickness stress description is required. In order to compute a more accurate stress state, we choose to use the following post-processing approach based on the equilibrium equations (as it has been proposed in [20, 21, 23, 25] in a finite element framework), relying on the higher regularity granted by IGA shape functions.

In an equilibrium state, the stresses inside the material should satisfy the equilibrium equation

[TABLE]

where ${\bm{\nabla}}\cdot$ is the divergence operator. In a 3D case, this reads in components as

[TABLE]

with

[TABLE]

By integrating along the thickness, we recover the out-of-plane stresses as

[TABLE]

It is important to note that, by introducing equations (8a) and (8b) into (8c), the component $\sigma_{33}$ can then be computed as

[TABLE]

The integral constants should be chosen to fulfill the boundary conditions at the top or bottom surfaces.

Remark 1**.**

It should be noted that the integration could be performed from both surfaces and then averaged in order to divide the resulting error by two. Although it is easy to compute the boundary conditions in this case, we chose to avoid this as it relies on the perfect knowledge of the stress boundary conditions on both surfaces (including the loaded one) which are not always available.

Assuming that the elasticity tensor $\mathbb{C}$ of equation (1) is constant, the derivatives of the in-plane components of the stress are computed from displacements as

[TABLE]

where

[TABLE]

Equation (11) clearly demonstrates the necessity of a highly regular displacement solution in order to recover a proper stress state. Such condition can be easily achieved using IGA with B-spline shape functions of order $p>2$ while in a standard finite element framework, continuity of the solution must be enhanced (see [23, 25]). One could thus retrieve a good approximation of the out-of-plane stress state once an accurate description of the in-plane stress state is available.

The displacement solution computed using an IGA method is a spline, which means that its post-processing and further handling of its derivatives can be efficiently performed. Splines can be indeed easily differentiated and many dedicated tools are available [30].

4.1 Numerical results

Figures 5, 6 and 7 show that applying the post-processing allows to greatly improve the quality of the results. Indeed, we observe that the post-processed stress state is very close to the analytical one. The correct stress state is recovered on the whole plate, and a great improvement is observed everywhere. The results were obtained using an in-plane mesh composed of $9\times 9$ quartic elements with cubic shape functions through the thickness.

In order to validate such an approach in a wider variety of cases, computations with several sets of number of layers for the laminate or a different ratio between the thickness of the plate and its length can be performed. Figure 8 gathers the results for different cases (varying the number of layers) with respect to the length-thickness ratio $S$ used in the computations. In each layer, a standard Gaussian integration with $p+1$ points (i.e., 5, in our case) along each in-plane direction has been used. The error is computed as:

[TABLE]

Several conclusions can be drawn from Figure 8. First, the post-processing method provides better results when the thickness ratio is larger. The method should therefore be used on slender plates rather than on thick stacks. Concerning the integration scheme, $p+1$ points per layer provide basically the same accuracy as using $p-1$ points per layer (i.e., 2 points, given the cubic through-the-thickness approximation). $p-2$ integration points per layer (corresponding to just 1 point in our case) lead in general to a loss of accuracy except, in our example, in the case of 34 layers (Figure 8(c)) where some super-convergence effect is present. Based on our tests, this “super-convergence” seems to be a special case happening only when a large unsymmetrical stack is involved.

Finally, the more layers in the stack, the better the results, since a stack with a large number of thin layers is indeed closer to a plate with average properties, whose solution would be perfectly captured by a standard isogeometric simulation with one high-order element through the thickness.

Therefore, the post-processing approach seems to be particularly suited to handle large and thin plates with a significant number of layers (whereas the first approach with one element per layer is far too costly to handle these kinds of problems, due to the huge number of elements and degrees of freedom). In such situations, this method provides results very close to the correct solution, at a fraction of the cost, representing a far better choice than a full 3D computation in terms of accuracy-to-cost ratio.

Although the results in this paper are obtained using a 9x9 elements in-plane mesh, Figure 9 shows that even using a single element mesh, the error remains acceptable in this case.

5 Conclusion

In this paper, we have presented a very simple approach for an accurate simulation of laminates combining a standard 3D coarse isogeometric analysis with a layerwise integration rule and a post-processing based on equilibrium equations. Using this approach, results close to the analytical solution are obtained, while using only few degrees of freedom through the thickness in contrast with a standard 3D layerwise approach.

Multiple numbers of layers (both even and odd) have been considered in our numerical tests and different thicknesses and numbers of integration points per layer have been studied. Whatever the number of layers, the method gives better results the thinner the layers are. Based on our results, a through-the-thickness integration rule involving $p-1$ points per layer can be considered as the default choice, since it grants a good accuracy even when the plate becomes thicker.

This method is easy to implement. Indeed, the post-processing is only based on the integration of equilibrium equations, and all the required components can be easily computed from the displacement solution. Furthermore, a standard 3D isogeometric computation is used, with a special set of integration points, but without any modification in the model. As the IGA code provides spline fields, its derivation and handling are straightforward and the whole process (coarse simulation plus post-processing) is far less time consuming than a full layerwise 3D approach.

Further research topics currently under investigation involve the extension of this approach to more complex problems involving curved geometries and/or large deformations. Moreover, more efficient through-the-thickness integration strategies or the combination with bivariate shells will be also considered in future studies.

Acknowledgment

P. Antolin was partially supported by the European Research Council through the H2020 ERC Advanced Grant 2015 n.694515 CHANGE. G. Sangalli was partially supported by the European Research Council through the FP7 ERC Consolidator Grant n.616563 HIGEOM.