Fast and accurate elastic analysis of laminated composite plates via isogeometric collocation and an equilibrium-based stress recovery approach

TL;DR

This paper presents a fast, accurate method for analyzing laminated composite plates using isogeometric collocation combined with an equilibrium-based stress recovery, achieving high accuracy with reduced computational effort.

Contribution

It introduces a novel combined approach of isogeometric collocation and equilibrium-based stress recovery for laminated plates, improving efficiency and accuracy over traditional methods.

Findings

High order convergence rates achieved

Accurate displacement and in-plane stress results

Effective out-of-plane stress recovery with minimal collocation points

Abstract

A novel approach which combines isogeometric collocation and an equilibrium-based stress recovery technique is applied to analyze laminated composite plates. Isogeometric collocation is an appealing strong form alternative to standard Galerkin approaches, able to achieve high order convergence rates coupled with a significantly reduced computational cost. Laminated composite plates are herein conveniently modeled considering only one element through the thickness with homogenized material properties. This guarantees accurate results in terms of displacements and in-plane stress components. To recover an accurate out-of-plane stress state, equilibrium is imposed in strong form as a post-processing correction step, which requires the shape functions to be highly continuous. This continuity demand is fully granted by isogeometric analysis properties, and excellent results are obtained using a minimal number of collocation points per direction, particularly for increasing values of length-to-thickness plate ratio and number of layers.

Figures

Tables

Table 1. Table 1: Material properties for 0°-oriented layers employed in the numerical tests.

$E_1$	$E_2$	$E_3$	$G_{23}$	$G_{13}$	$G_{12}$	${\nu }_{23}$	${\nu }_{13}$	${\nu }_{12}$
[GPa]	[GPa]	[GPa]	[GPa]	[GPa]	[GPa]	[-]	[-]	[-]
25000	1000	1000	200	500	500	0.25	0.25	0.25

Table 2. Table 2: Simply supported composite plate under sinusoidal load with a number of layers equal to 11. Out-of-plane stress state maximum relative error with respect to Pagano’s solution  [ 60 ] at x 1 = x 2 = 0.25 ​ L subscript 𝑥 1 subscript 𝑥 2 0.25 𝐿 x_{1}=x_{2}=0.25L . Comparing the isogeometric collocation-based homogenized single element approach (IGA-C) and the coupled post-processing technique (IGA-C+PP) for different approximation degrees.

Degree			$p=q=6$, $r=4$			$p=q=r=6$
S	Method	$e\left({\sigma }_{13}\right)$	$e\left({\sigma }_{23}\right)$	$e\left({\sigma }_{33}\right)$	$e\left({\sigma }_{13}\right)$	$e\left({\sigma }_{23}\right)$	$e\left({\sigma }_{33}\right)$
		[%]	[%]	[%]	[%]	[%]	[%]
20	IGA-C	97.6	56.7	6.34	96.6	56.1	6.31
	IGA-C+PP	0.31	2.94	0.90	1.97	1.20	0.05
30	IGA-C	98.7	55.6	6.36	98.3	55.4	6.34
	IGA-C+PP	0.16	1.34	0.47	0.91	0.57	0.08
40	IGA-C	99.2	55.3	6.37	98.9	55.2	6.36
	IGA-C+PP	0.07	0.78	0.34	0.50	0.35	0.12
50	IGA-C	99.4	55.1	6.38	99.2	55.1	6.38
	IGA-C+PP	0.03	0.52	0.29	0.30	0.25	0.15

Table 3. Table 3: Simply supported composite plate under sinusoidal load with a number of layers equal to 3. Out-of-plane stress state maximum relative error with respect to Pagano’s solution  [ 60 ] at x 1 = x 2 = 0.25 ​ L subscript 𝑥 1 subscript 𝑥 2 0.25 𝐿 x_{1}=x_{2}=0.25L . Comparing the isogeometric collocation-based homogenized single element approach (IGA-C) and the coupled post-processing technique (IGA-C+PP) for different approximation degrees.

Degree			$p=q=6$, $r=4$			$p=q=r=6$
S	Method	$e\left({\sigma }_{13}\right)$	$e\left({\sigma }_{23}\right)$	$e\left({\sigma }_{33}\right)$	$e\left({\sigma }_{13}\right)$	$e\left({\sigma }_{23}\right)$	$e\left({\sigma }_{33}\right)$
		[%]	[%]	[%]	[%]	[%]	[%]
20	IGA-C	292	57.2	5.80	291	57.2	5.79
	IGA-C+PP	10.4	3.16	0.54	11.9	1.41	0.33
30	IGA-C	311	57.5	5.77	311	57.5	5.77
	IGA-C+PP	5.05	1.40	0.28	5.75	0.63	0.11
40	IGA-C	319	57.6	5.76	319	57.6	5.76
	IGA-C+PP	2.91	0.81	0.21	3.32	0.38	0.02
50	IGA-C	323	57.6	5.76	322	57.6	5.76
	IGA-C+PP	1.87	0.54	0.20	2.14	0.26	0.07

Table 4. Table 4: Simply supported composite plate under sinusoidal load with a number of layers equal to 33. Out-of-plane stress state maximum relative error with respect to Pagano’s solution  [ 60 ] at x 1 = x 2 = 0.25 ​ L subscript 𝑥 1 subscript 𝑥 2 0.25 𝐿 x_{1}=x_{2}=0.25L . Comparing the isogeometric collocation-based homogenized single element approach (IGA-C) and the coupled post-processing technique (IGA-C+PP) for different approximation degrees.

Degree			$p=q=6$, $r=4$			$p=q=r=6$
S	Method	$e\left({\sigma }_{13}\right)$	$e\left({\sigma }_{23}\right)$	$e\left({\sigma }_{33}\right)$	$e\left({\sigma }_{13}\right)$	$e\left({\sigma }_{23}\right)$	$e\left({\sigma }_{33}\right)$
		[%]	[%]	[%]	[%]	[%]	[%]
20	IGA-C	81.6	69.7	6.33	80.7	68.9	6.33
	IGA-C+PP	1.16	2.21	0.93	0.54	0.50	0.07
30	IGA-C	81.5	69.0	6.34	81.2	68.7	6.34
	IGA-C+PP	0.53	1.01	0.48	0.23	0.25	0.09
40	IGA-C	81.5	68.7	6.35	81.3	68.6	6.34
	IGA-C+PP	0.32	0.59	0.34	0.11	0.16	0.12
50	IGA-C	81.6	68.6	6.35	81.4	68.5	6.35
	IGA-C+PP	0.23	0.40	0.30	0.05	0.13	0.16

Equations

$$\Theta =\{{\theta }_1^d,...,{\theta }_{m_d+p_d+1}^d\}d=1,...,d_p,$$

$$B_{\text{i},\text{p}}(\mathbf{\xi })=\prod _{d=1}^{d_p}{N}_{i_d,p_d}({\xi }^d),$$

$$\text{S}(\mathbf{\xi })=\sum _{\text{i}}{B}_{\text{i},\text{p}}(\mathbf{\xi }){\text{P}}_{\text{i}},$$

$$R_{\text{i},\text{p}}(\mathbf{\xi })=\frac{B_{\text{i},\text{p}}(\mathbf{\xi })w_{\text{i}}}{{\sum }_{\text{j}}{B}_{\text{j},\text{p}}(\mathbf{\xi })w_{\text{j}}}$$

$$\text{S}(\mathbf{\xi })=\sum _{\text{i}}{R}_{\text{i},\text{p}}(\mathbf{\xi }){\text{P}}_{\text{i}}.$$

$${\overline{\theta }}_i^d=\frac{{\theta }_{i+1}^d+{\theta }_{i+2}^d+...+{\theta }_{i+p}^d}{p_d}i=1,...,m^d,$$

$${\sigma }_{ij,j}+b_i=0\text{in}\Omega$$

$${\sigma }_{ij}{n}_j=t_i\text{on}{\Gamma }_N$$

$$u_i={\overline{u}}_i\text{on}{\Gamma }_D$$

$${\epsilon }_{ij}=\frac{u_{i,j}+u_{j,i}}{2},$$

$${\sigma }_{ij}={\mathbb{C}}_{ijkm}{\epsilon }_{km},$$

$$\mathbb{C}=\left[\begin{array}{cccccc}{\mathbb{C}}_{11}& {\mathbb{C}}_{12}& {\mathbb{C}}_{13}& 0& 0& 0\\ & {\mathbb{C}}_{22}& {\mathbb{C}}_{23}& 0& 0& 0\\ & & {\mathbb{C}}_{33}& 0& 0& 0\\ & symm& & {\mathbb{C}}_{44}& 0& 0\\ & & & & {\mathbb{C}}_{55}& 0\\ & & & & & {\mathbb{C}}_{66}\end{array}\right]=\text{bBigg@}17.5[\begin{array}{cccccc}\frac{1}{E_1}& -\frac{{\nu }_{12}}{E_1}& -\frac{{\nu }_{13}}{E_1}& 0& 0& 0\\ & \frac{1}{E_2}& -\frac{{\nu }_{23}}{E_2}& 0& 0& 0\\ & & \frac{1}{E_3}& 0& 0& 0\\ & symm& & \frac{1}{G_{23}}& 0& 0\\ & & & & \frac{1}{G_{13}}& 0\\ & & & & & \frac{1}{G_{12}}\end{array}\text{bBigg@}17.5{]}^{-1}.$$

$$\text{u}(\mathbf{\xi })=R_{\text{i},\text{p}}(\mathbf{\xi }){\hat{\text{u}}}_{\text{i}},$$

$$\text{v}(\mathbf{\xi })=R_{\text{i},\text{p}}(\mathbf{\xi }){\hat{\text{v}}}_{\text{i}},$$

$$\text{w}(\mathbf{\xi })=R_{\text{i},\text{p}}(\mathbf{\xi }){\hat{\text{w}}}_{\text{i}}.$$

$$\left[\begin{array}{ccc}{\text{K}}_{11}(\mathbf{\tau })& {\text{K}}_{12}(\mathbf{\tau })& {\text{K}}_{13}(\mathbf{\tau })\\ & {\text{K}}_{22}(\mathbf{\tau })& {\text{K}}_{23}(\mathbf{\tau })\\ symm& & {\text{K}}_{33}(\mathbf{\tau })\end{array}\right]\cdot \left(\begin{array}{c}{\hat{\text{u}}}_{\text{i}}\\ {\hat{\text{v}}}_{\text{i}}\\ {\hat{\text{w}}}_{\text{i}}\end{array}\right)=-\text{b}(\mathbf{\tau }),\forall \mathbf{\tau }\in \Omega ,$$

$${\text{K}}_{11}(\mathbf{\tau })={\mathbb{C}}_{11}\frac{{\partial }^2{R}_{\text{i},\text{p}}(\mathbf{\tau })}{\partial {x_1}^2}+{\mathbb{C}}_{66}\frac{{\partial }^2{R}_{\text{i},\text{p}}(\mathbf{\tau })}{\partial {x_2}^2}+{\mathbb{C}}_{55}\frac{{\partial }^2{R}_{\text{i},\text{p}}(\mathbf{\tau })}{\partial {x_3}^2},$$

$${\text{K}}_{22}(\mathbf{\tau })={\mathbb{C}}_{66}\frac{{\partial }^2{R}_{\text{i},\text{p}}(\mathbf{\tau })}{\partial {x_1}^2}+{\mathbb{C}}_{22}\frac{{\partial }^2{R}_{\text{i},\text{p}}(\mathbf{\tau })}{\partial {x_2}^2}+{\mathbb{C}}_{44}\frac{{\partial }^2{R}_{\text{i},\text{p}}(\mathbf{\tau })}{\partial {x_3}^2},$$

$${\text{K}}_{33}(\mathbf{\tau })={\mathbb{C}}_{55}\frac{{\partial }^2{R}_{\text{i},\text{p}}(\mathbf{\tau })}{\partial {x_1}^2}+{\mathbb{C}}_{44}\frac{{\partial }^2{R}_{\text{i},\text{p}}(\mathbf{\tau })}{\partial {x_2}^2}+{\mathbb{C}}_{33}\frac{{\partial }^2{R}_{\text{i},\text{p}}(\mathbf{\tau })}{\partial {x_3}^2},$$

$${\text{K}}_{23}(\mathbf{\tau })=({\mathbb{C}}_{23}+{\mathbb{C}}_{44})\frac{{\partial }^2{R}_{\text{i},\text{p}}(\mathbf{\tau })}{\partial x_2\partial x_3},$$

$${\text{K}}_{13}(\mathbf{\tau })=({\mathbb{C}}_{13}+{\mathbb{C}}_{55})\frac{{\partial }^2{R}_{\text{i},\text{p}}(\mathbf{\tau })}{\partial x_1\partial x_3},$$

$${\text{K}}_{12}(\mathbf{\tau })=({\mathbb{C}}_{12}+{\mathbb{C}}_{66})\frac{{\partial }^2{R}_{\text{i},\text{p}}(\mathbf{\tau })}{\partial x_1\partial x_2},$$

$$\left[\begin{array}{ccc}{\tilde{\text{K}}}_{11}(\mathbf{\tau })& {\tilde{\text{K}}}_{12}(\mathbf{\tau })& {\tilde{\text{K}}}_{13}(\mathbf{\tau })\\ & {\tilde{\text{K}}}_{22}(\mathbf{\tau })& {\tilde{\text{K}}}_{23}(\mathbf{\tau })\\ symm& & {\tilde{\text{K}}}_{33}(\mathbf{\tau })\end{array}\right]\cdot \left(\begin{array}{c}{\hat{\text{u}}}_{\text{i}}\\ {\hat{\text{v}}}_{\text{i}}\\ {\hat{\text{w}}}_{\text{i}}\end{array}\right)=\text{t}(\mathbf{\tau }),\forall \mathbf{\tau }\in {\Gamma }_N$$

$${\tilde{\text{K}}}_{11}(\mathbf{\tau })={\mathbb{C}}_{11}\frac{\partial R_{\text{i},\text{p}}(\mathbf{\tau })}{\partial x_1}{n}_1+{\mathbb{C}}_{66}\frac{\partial R_{\text{i},\text{p}}(\mathbf{\tau })}{\partial x_2}{n}_2+{\mathbb{C}}_{55}\frac{\partial R_{\text{i},\text{p}}(\mathbf{\tau })}{\partial x_3}{n}_3,$$

$${\tilde{\text{K}}}_{22}(\mathbf{\tau })={\mathbb{C}}_{66}\frac{\partial R_{\text{i},\text{p}}(\mathbf{\tau })}{\partial x_1}{n}_1+{\mathbb{C}}_{22}\frac{\partial R_{\text{i},\text{p}}(\mathbf{\tau })}{\partial x_2}{n}_2+{\mathbb{C}}_{44}\frac{\partial R_{\text{i},\text{p}}(\mathbf{\tau })}{\partial x_3}{n}_3,$$

$${\tilde{\text{K}}}_{33}(\mathbf{\tau })={\mathbb{C}}_{55}\frac{\partial R_{\text{i},\text{p}}(\mathbf{\tau })}{\partial x_1}{n}_1+{\mathbb{C}}_{44}\frac{\partial R_{\text{i},\text{p}}(\mathbf{\tau })}{\partial x_2}{n}_2+{\mathbb{C}}_{33}\frac{\partial R_{\text{i},\text{p}}(\mathbf{\tau })}{\partial x_3}{n}_3,$$

$${\tilde{\text{K}}}_{23}(\mathbf{\tau })={\mathbb{C}}_{23}\frac{\partial R_{\text{i},\text{p}}(\mathbf{\tau })}{\partial x_3}{n}_2+{\mathbb{C}}_{44}\frac{\partial R_{\text{i},\text{p}}(\mathbf{\tau })}{\partial x_2}{n}_3,$$

$${\tilde{\text{K}}}_{13}(\mathbf{\tau })={\mathbb{C}}_{13}\frac{\partial R_{\text{i},\text{p}}(\mathbf{\tau })}{\partial x_3}{n}_1+{\mathbb{C}}_{55}\frac{\partial R_{\text{i},\text{p}}(\mathbf{\tau })}{\partial x_1}{n}_3,$$

$${\tilde{\text{K}}}_{12}(\mathbf{\tau })={\mathbb{C}}_{12}\frac{\partial R_{\text{i},\text{p}}(\mathbf{\tau })}{\partial x_2}{n}_1+{\mathbb{C}}_{66}\frac{\partial R_{\text{i},\text{p}}(\mathbf{\tau })}{\partial x_1}{n}_2.$$

$${\overline{\mathbb{C}}}_{11}=\sum _{k=1}^N{\overline{t}}_k{\mathbb{C}}_{11}^{(k)}+\sum _{k=2}^N({\mathbb{C}}_{13}^{(k)}-{\overline{\mathbb{C}}}_{13}){\overline{t}}_k\frac{({\mathbb{C}}_{13}^{(1)}-{\mathbb{C}}_{13}^{(k)})}{{\mathbb{C}}_{33}^{(k)}}$$

$${\overline{\mathbb{C}}}_{12}=\sum _{k=1}^N{\overline{t}}_k{\mathbb{C}}_{12}^{(k)}+\sum _{k=2}^N({\mathbb{C}}_{13}^{(k)}-{\overline{\mathbb{C}}}_{13}){\overline{t}}_k\frac{({\mathbb{C}}_{23}^{(1)}-{\mathbb{C}}_{23}^{(k)})}{{\mathbb{C}}_{33}^{(k)}}$$

$${\overline{\mathbb{C}}}_{13}=\sum _{k=1}^N{\overline{t}}_k{\mathbb{C}}_{13}^{(k)}+\sum _{k=2}^N({\mathbb{C}}_{33}^{(k)}-{\overline{\mathbb{C}}}_{33}){\overline{t}}_k\frac{({\mathbb{C}}_{13}^{(1)}-{\mathbb{C}}_{13}^{(k)})}{{\mathbb{C}}_{33}^{(k)}}$$

$${\overline{\mathbb{C}}}_{22}=\sum _{k=1}^N{\overline{t}}_k{\mathbb{C}}_{22}^{(k)}+\sum _{k=2}^N({\mathbb{C}}_{23}^{(k)}-{\overline{\mathbb{C}}}_{23}){\overline{t}}_k\frac{({\mathbb{C}}_{23}^{(1)}-{\mathbb{C}}_{23}^{(k)})}{{\mathbb{C}}_{33}^{(k)}}$$

$${\overline{\mathbb{C}}}_{23}=\sum _{k=1}^N{\overline{t}}_k{\mathbb{C}}_{23}^{(k)}+\sum _{k=2}^N({\mathbb{C}}_{33}^{(k)}-{\overline{\mathbb{C}}}_{33}){\overline{t}}_k\frac{({\mathbb{C}}_{23}^{(1)}-{\mathbb{C}}_{23}^{(k)})}{{\mathbb{C}}_{33}^{(k)}}$$

$$\displaystyle\overline{\mathbb{C}}_{33}=\frac{1}{\bigg{(}\sum_{k=1}^{N}\cfrac{\overline{t}_{k}}{\mathbb{C}_{33}^{(k)}}\bigg{)}}$$

$$\displaystyle\overline{\mathbb{C}}_{44}=\frac{\bigg{(}\sum_{k=1}^{N}\cfrac{\overline{t}_{k}\mathbb{C}_{44}^{(k)}}{\Delta_{k}}\bigg{)}}{\Delta},\hskip 14.22636pt\Delta=\bigg{(}\sum_{k=1}^{N}\frac{\overline{t}_{k}\mathbb{C}_{44}^{(k)}}{\Delta_{k}}\bigg{)}\bigg{(}\sum_{k=1}^{N}\frac{\overline{t}_{k}\mathbb{C}_{55}^{(k)}}{\Delta_{k}}\bigg{)}$$

$$\displaystyle\overline{\mathbb{C}}_{55}=\frac{\bigg{(}\sum_{k=1}^{N}\cfrac{\overline{t}_{k}\mathbb{C}_{55}^{(k)}}{\Delta_{k}}\bigg{)}}{\Delta},\hskip 14.22636pt\Delta_{k}=\mathbb{C}_{44}^{k}\mathbb{C}_{55}^{k}$$

$${\overline{\mathbb{C}}}_{66}=\sum _{k=1}^N{\overline{t}}_k{\mathbb{C}}_{66}^{(k)}$$

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

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

$${\sigma }_{13,1}+{\sigma }_{23,2}+{\sigma }_{33,3}=-b_3.$$

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

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

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

Full text

1 Introduction

Composite materials consist of two or more materials which combined present enhanced properties that could not be acquired employing any of the constituents alone (see, e.g., [28, 41, 65, 74] and references therein). The interest for composite structures in the engineering community has constantly grown in recent years due to their appealing mechanical properties such as increased stiffness and strength, reduced weight, improved corrosion and wear resistance, just to recall some of them. The majority of man-made composite materials consists of reinforced fibers embedded in a base material, called matrix (see, e.g., [34, 65]). The matrix material keeps the fibers together, acts as a load-transfer medium between fibers, process which takes place through shear stresses, and protects those elements from being exposed to the environment, while the resistence properties of composites are given by the fibers which are stiffer and stronger than the soft matrix. In this paper we focus on laminated composite materials, which are formed by a collection of building blocks or plies, stacked to achieve the desired stiffness and thickness. For this kind of structures also simple loading conditions, such as traction or bending, cause a complex 3D stress state because of the difference in the material properties between the layers, which may lead to delamination, consequently requiring an accurate stress evaluation through the thickness (see, e.g., [65, 70]). As an alternative to two-dimensional theories, often insufficiently accurate to depict delamination and interlaminar damage, and to layerwise theories, which typically show an high computational cost, a novel method combining an isogeometric analysis (IGA) Galerkin approach with a stress-recovery technique has been recently proposed in [23]. Introduced in 2005 by Hughes et al. [36], IGA aims at integrating design and analysis employing shape functions typically belonging to Computer Aided Design field (such as, B-Splines and NURBS). Using the same shape functions to approximate both geometry and field variables leads to a cost-saving simplification of expensive mesh generation and refinement processes required by standard finite element analyisis. One of the most important features of IGA is the high-regularity of its basis functions leading to superior approximation properties. IGA proved to be successful in a wide variety of problems ranging from solids and structures (see, e.g., [5, 11, 13, 22, 24, 37, 39, 49, 58]) to fluids (see, e.g., [1, 12, 30, 50]), fluid-structure interaction (see, e.g., [10, 35]), opening also the door to geometrically flexible discretizations of higher-order partial differential equations in primal form as in [3, 29, 43, 51]. However, a well-known important issue of IGA is related to the development of efficient integration rules when higher-order approximations are employed (see, e.g., [7, 26, 38, 67]). In attempt to address this problem taking full advantage of the special possibilities offered by IGA, isogeometric collocation (IGA-C) schemes have been proposed in [4]. The aim was to optimize the computational cost still relying on IGA geometrical flexibility and accuracy. Collocation main idea, in contrast to Galerkin-type formulations, consists in the discretization of the governing partial differential equations in strong form, evaluated at suitable points. Since integration is not required, isogeometric collocation results in a very fast method providing superior performance in terms of accuracy-to-computational effort ratio with respect to Galerkin formulations, in particular when higher-order approximation degrees are adopted (see [68]). Isogeometric collocation has been particularly successful in the context of structural elements, where isogeometric collocation has proven to be particularly stable in the context of mixed methods. In particular, Bernoulli-Euler beam and Kirchhoff plate elements have been proposed [64], while mixed formulations both for Timoshenko initially-straight planar [17] and non-prismatic [9] beams as well as for curved spatial rods [8] have been introduced and studied, and then effectively extended to the geometrically nonlinear case [47, 53, 54, 55, 75, 77]. Isogeometric collocation has been moreover successfully applied to the solution of Reissner-Mindlin plate problems in [44], and new formulations for shear-deformable beams [45, 47], as well as shells [46, 56] have been solved also via IGA collocation. Since its introduction, many promising significant works on isogeometric collocation methods have been published also in other fields, including phase-field modeling [31], contact [21, 48, 76] and poromechanics [59]. Moreover, combinations with different spline spaces, like hierarchical splines, generalized B-Splines, and T-Splines, have been successfully tested in [14, 52, 68], while alternative effective selection strategies for collocation points have been proposed in [2, 31, 57]. IGA-Galerkin methods have already been used to solve composite laminate problems, especially relying on high-order theories for enhanced plate and shell theories [27, 42, 66, 69, 72]. Recently an isogeometric collocation numerical formulation has been proposed [61] to study Reissner-Mindlin composite plates. Other Galerkin methods [32, 33, 66] compute instead a full 3D stress state using isogeometric analysis, applying a layerwise technique. This method can also be applied to isogeometric collocation adopting a multipatch approach, which models each layer as a patch (see Figure 1(a)), enforcing normal stress continuity at the inter-patches boundaries [6]. Clearly a layerwise method exploits a number of degrees of freedom directly proportional to the number of layers, inevitably leading to high computational costs. In this paper we apply a single patch 3D isogeometric collocation method to analyze the behavior of composite plates. We adopt a homogenized single-element approach (see Figure 1(b)), which conveniently uses one element through the thickness, coupled with a post-processing technique in order to recover a proper out-of-plane stress state. This method is significantly less expensive compared to a layerwise approach since employs a considerably lower number of degrees of freedom.

The post-processing approach, first proposed in [23], takes inspiration from recovery techniques which can be found in [16, 19, 25, 42, 63, 73] and is based on the direct integration of the equilibrium equations to compute the out-of-plane stress components from the in-plane ones directly derived from a coarse displacement solution.

The structure of the paper is organized as follows. In Section 2 the fundamental concepts of multivariate B-Splines and NURBS are presented, followed by an introduction to isogeometric collocation and a description of our IGA-C scheme for orthotropic elasticity. In Section 3 we define our isogeometric collocation strategy to study laminated plates, which combines a homogenized single-element approach with an equilibrium-based stress recovery technique. In Section 4 we present our reference test case and provide results for the single-element approach. Several numerical benchmarks are displayed, which show a significant improvement between non-treated and post-processed out-of-plane stress components. Finally we provide some mesh sensitivity tests considering an increasing length-to-thickness ratio and numbers of layers to show the effectiveness of the method. We draw our conclusions in Section 5.

2 Isogeometric Collocation: Basics and application to orthotropic elasticity

In this section we introduce the notions of multivariate B-Splines and NURBS, provide some details regarding isogeometric collocation and describe our collocation scheme in the context of linear orthotropic elasticity.

2.1 Multivariate B-Splines and NURBS

In the following, we introduce the basic definitions and notations about multivariate B-Splines and NURBS. For further details, readers may refer to [15, 36, 62], and references therein. Multivariate B-Splines are generated through the tensor product of univariate B-Splines. We denote with $d_{p}$ the dimension of the parametric space and therefore $d_{p}$ univariate knot vectors have to be introduced as

[TABLE]

where $p_{d}$ represents the polynomial degree in the parametric direction $d$, and $m_{d}$ is the associated number of basis functions. Given the univariate basis functions $N^{d}_{i_{d},p_{d}}$ associated to each parametric direction $\xi^{d}$, the multivariate basis functions $B_{\textbf{i},\textbf{p}}(\bm{\xi})$ are obtained as:

[TABLE]

where $\textbf{i}=\{i_{1},...,i_{d_{p}}\}$ plays the role of a multi-index which describes the considered position in the tensor product structure, $\textbf{p}=\{{p_{1},...,p_{d}}\}$ indicates the polynomial degrees, and $\bm{\xi}=\{\xi^{1},...,\xi^{d_{p}}\}$ represents the vector of the parametric coordinates in each parametric direction $d$. B-Spline multidimensional geometries are built from a linear combination of multivariate B-Spline basis functions as follows

[TABLE]

where the coefficients $\textbf{P}_{\textbf{i}}\in\mathbb{R}^{d_{s}}$ of the linear combination are the so-called control points ($d_{s}$ is the dimension of the physical space) and the summation is extended to all combinations of the multi-index i. NURBS geometries in $\mathbb{R}^{d_{s}}$ are instead obtained from a projective transformation of their B-Spline counterparts in $\mathbb{R}^{d_{s}+1}$. Defining $w_{\textbf{i}}$ as the collection of weights according to the multi-index i, multivariate NURBS basis functions are obtained as

[TABLE]

and NURBS multidimensional geometries are built as

[TABLE]

2.2 An Introduction to Isogeometric collocation

Collocation methods have been introduced within isogeometric analysis as an attempt to address a well-known important issue of early IGA-Galerkin formulations, related to the development of efficient integration rules for higher-order approximations. In fact, element-wise Gauss quadrature, typically used for finite elements and originally adopted for Galerkin-based IGA, does not properly take into account inter-element higher continuity leading to sub-optimal array formation and assembly costs, significantly affecting the performance of IGA methods. Isogeometric collocation aimed at optimizing computational cost, since it may be viewed as a variant of one-point quadrature numerical scheme, still taking advantage of IGA geometrical flexibility and accuracy. Collocation methods are based on the direct discretization in strong form of the differential equations governing the problem evaluated at suitable points. The isoparametric paradigm is adopted and the same basis functions are used to describe both geometry and problem unknowns. Once the approximations are carried out, as in a typical Galerkin-IGA context, by means of a linear combinations of IGA basis functions and control variables, the discrete differential equations are collocated at each collocation point. Consequently a delicate issue is represented by the determination of suitable collocation points. A widespread approach which is proposed in the engineering literature is to collocate at the images of Greville abscissae (see, e.g., [40]), but this represents just the simplest possible option (see, e.g., [18, 20] for alternative choices). Along each parametric direction $d$, Greville abscissae consist of a set of $m^{d}$ points, obtained from the knot vector components, $\theta^{d}_{i}$, as

[TABLE]

$p_{d}$ being the degree of approximation. Since the approximation is performed through direct collocation of the differential equations, no integrals need to be computed and consequently, evaluation and assembly operations lead to a significantly reduced computational cost.

2.3 Numerical formulation for orthotropic elasticity

Once a strategy to select collocation points and compute IGA shape functions is set, a proper description of the equations in strong form for the problem under examination is required, as mentioned in Section 2.2. We therefore recall the classical elasticity problem in strong form considering a small strain regime and detail equilibrium equations using Einstein’s notation (7). The following notations are used: $\Omega\subset\mathbb{R}^{3}$, is an open bounded domain, representing an elastic three-dimensional body, $\Gamma_{N}$ and $\Gamma_{D}$ are defined as boundary portions subjected respectively to Neumann and Dirichlet conditions such that $\Gamma_{N}\cup\Gamma_{D}=\partial\Omega$ and $\Gamma_{N}\cap\Gamma_{D}=\emptyset$. Accordingly, the equilibrium equations and the corresponding boundary conditions are:

[TABLE]

where $\sigma_{ij}$ and $u_{i}$ represent respectively the Cauchy stress and displacement components, while $b_{i}$ and $t_{i}$ the volume and traction forces, $n_{j}$ the outward normal, and $\overline{u}_{i}$ the prescribed displacements. The elasticity problem is finally completed by the kinematic relations in small strain

[TABLE]

as well as by the constitutive equations

[TABLE]

where $\mathbb{C}_{ijkm}$ is the fourth order elasicity tensor.

As we described in Section 1, the basic building block of a laminate is a lamina, i.e., a flat arrangement of unidirectional fibers, considering the simplest case, embedded in a matrix. In order to increase the composite resistance properties cross-ply laminates can be employed (i.e., all the plies used to form the composite stacking sequence are piled alternating different fiber layers orientations) in which all unidirectional layers are individually orthotropic. Since the proposed collocation approach uses one element through the thickness to model the composite plate as a homogenized single building block, we focus in this section on the collocation formulation for a plate formed by only one orthotropic elastic lamina. Considering three mutually orthogonal planes of material symmetry for each ply, the number of elastic coefficients of the fourth order elasticity tensor $\mathbb{C}_{ijkm}$ is reduced to 9 in Voigt notation, that can be expressed in terms of engineering constants as

[TABLE]

The displacement field is then approximate as a linear combination of NURBS multivariate shape functions and control points as follows

[TABLE]

Having defined $\bm{\tau}$ as the matrix of collocation points, we insert the approximations (11) into kinematics equations (8) and we combine the obtained expressions with the constitutive relations (9). Finally we substitute into equilibrium equations (7a) obtaining

[TABLE]

and substituting in (7b) we obtain:

[TABLE]

As we can see from equations (13), Neumann boundary conditions are directly imposed as strong equations at the collocation points belonging to the boundary surface (see, [6, 21]), with the usual physical meaning of prescribed boundary traction.

3 An IGA collocation approach to model 3D composite plates

In this section we describe our IGA 3D collocation strategy to model composite plates. The proposed method, known as single element approach, relies on a homogenization technique combined with a post-processing approach based on the imposition of equilibrium equations in strong form.

3.1 Single-element approach

The single-element approach considers the plate discretized by a single element through the thickness, which strongly reduces the number of degrees of freedom with respect to layerwise methods. The material matrix is therefore homogenized to account for the presence of the layers as Figure 1(b) clearly describes.

Remark 1**.**

Considering a single-element homogenized approach is effective only for through-the-thickness symmetric layer distributions, as for non-symmetric ply stacking sequences the plate middle plane is not balanced. In the case of non-symmetric layer distributions this technique is still applicable when the stacking sequence can be split into two symmetric piles, using one element per homogenized stack with a $C^{0}$ interface.

This method provides accurate results only in terms of displacements and in-plane stress components and, in order to recover a proper out-of-plane stress state, following [23], we propose to couple it with a post-processing technique. To characterize the variation of the material properties from layer to layer, we homogenize the constitutive behavior to create an equivalent single-layer laminate, referring to [71], where explicit expressions for the effective elastic constants of the equivalent laminate are given as

[TABLE]

where $\mathbb{C}_{ij}^{(k)}$ represents the $ij$-th component of the fourth order elasticity tensor in Voigt notation (10) for the $k$-th layer and $\overline{t}_{k}=\cfrac{t_{k}}{h}$ stands for the volume fraction of the $k$-th lamina, $h$ being the total thickness and $t_{k}$ the $k$-th thickness.

3.1.1 Post-processing step: Reconstruction from Equilibrium

As interlaminar delamination and other fracture processes rely mostly on out-of-plane components, a proper through-the-thickness stress description is required. In order to recover a more accurate stress state, we perform a post-processing step based on the equilibrium equations, following [23], relying on the higher regularity granted by IGA shape functions. This procedure, which takes its roots in [19, 25, 63, 73], has already been proven to be successful for IGA-Galerkin. Inside the plate the stresses should satisfy the equilibrium equation (7a) that can be expanded as

[TABLE]

Assuming the in-plane stress components to well approximate the laminate behaviour, as it will be shown in Section 4, we can integrate equation 15a and 15b along the thickness, recovering the out-of-plane shear stresses as

[TABLE]

where $\zeta$ represents the coordinate along the thickness direction.

Finally we can insert equations (16a) and (16b) into (15c), recovering the $\sigma_{33}$ component as

[TABLE]

Following [23], the integral constants are chosen to fulfil the boundary conditions at the top or bottom surfaces.

Recalling that

[TABLE]

where the homogenized elasticity tensor $\overline{\mathbb{C}}$ is constant, it is clear the necessity of a highly regular displacement solution in order to recover a proper stress state. Such a condition can be easily achieved using isogeometric collocation, due to the possibility to benefit from the high regularity of B-Splines or NURBS. We also remark that the proposed method strongly relies on the possibility to obtain an accurate description (with a relatively coarse mesh) of the in-plane stress state.

4 Numerical tests

In this section, to assess whether the proposed method can effectively reproduce composite plates behaviour, we consider a classical benchmark problem [60] and we address different aspects such as the effectiveness of the proposed post-processing step, the method sensitivity to parameters of interest (i.e., number of layers and length-to-thickness ratio), and its convergence.

4.1 Reference solution: the Pagano layered plate

A square laminated composite plate of total thickness $t$ made of $N$ orthotropic layers is considered. This structure is simply supported and a normal sinusoidal traction is applied on the upper surface, while the lower surface is traction-free, as shown in Figure 2.

In the proposed numerical tests we consider different numbers of layers, namely 3, 11 and 33. The thickness of every single layer is set to 1 mm, and the edge length, $L$, is chosen to be $S$ times larger than the total thickness $t$ of the laminate. Different choices of length-to-thickness ratio are considered (i.e., 20, 30, 40, and 50) which allow to draw interesting considerations about the laminate behaviour in the proposed convergence tests. For all examples we consider the same loading conditions proposed by Pagano, i.e., a double sinus with periodicity equal to twice the length of the plate. As depicted in Figure 2 the laminated plate is composed of layers organized in an alternated distribution of orthotropic plies (i.e., a 0°/90° stacking sequence in our case). Layer material parameters considered in the numerical tests are summarized in Table 1 for 0°-oriented plies.

The Neumann boundary conditions on the plate surfaces $x_{3}=\pm\dfrac{t}{2}$ are

[TABLE]

where $\sigma_{0}=$ 1 MPa.

The simple support edge conditions are taken as

[TABLE]

All results are then expressed in terms of the following normalized stress components

[TABLE]

4.2 Post-processed out-of-plane stresses

In this section, we comment the results obtained using the proposed IGA-collocation approach, as compared with Pagano’s analytical solution [60]. To give an idea of the improvement granted by the post-processing of out-of-plane stress components, in Figures 3 and 4 we compare the reference solution with non-treated and post-processed results for the cases with 3 and 11 layers, considering a length to thickness ratio $S=20$. All numerical simulations are carried out using an in-plane degree of approximation $p=q=6$ and 10 collocation points for each in-plane parametric direction, while we use an approximation degree $r=4$ and one element through the thickness (i.e., $r+1$ collocation points). The sampling point where we show results is located at $x_{1}=x_{2}=0.25L$. For both considered cases the in-plane stresses show a good behaviour, as expected, while the out-of-plane stress components, without a post-processing treatment, are erroneusly discontinuous. The proposed results clearly show the improvement granted by the post-process of out-of-plane components.

To show the effect of post-processing at different locations of the plate, in Figures 5-7 the out-of-plane stress state profile is recovered sampling the laminae every quarter of length in both in-plane directions, for the case of a length-to-thickness ratio equal to 20 and 11 layers.

4.3 Convergence behaviour

In order to validate the proposed approach in a wider variety of cases, computations with a different ratio between the thickness of the plate and its length are performed respectively for 3, 11, and 33 layers, considering an increasing number of knot spans. Figures 8 and 9 assess the convergence behaviour of the method, adopting the following error definition

[TABLE]

Note that relation (22) is used only to estimate the error inside the domain to avoid indeterminate forms. Different combinations of degree of approximations have been also considered. A poorer out-of-plane stress approximation is obtained using a degree equal to 4 in every direction, and, in addition, with this choice locking phenomena may occur for increasing values of length-to-thickness ratio. Therefore, we conclude that using a degree of approximation equal to 6 in-plane and equal to 4 through the thickness seems to be a reasonable choice to correctly reproduce the 3D stress state. Using instead uniform approximation degrees $p=q=r=6$ does not seem to significantly improve the results (see Figures 8, 9, and Table 2). The post-processing method provides better results for increasing values of length-to-thickness ratio and number of layers and therefore proves to be particularly convenient for very large and thin plates. This is clear since a laminae with a large number of thin layers resembles a plate with average properties. What really stands out looking at the displayed mesh sensitivity results, is the fact that collocation perfectly captures the plates behaviour not only using one element through the thickness but also employing only one knot span in the plane of the plate. A single element of degrees $p=q=6$ and $r=4$, comprising a total of 7x7x5 collocation points, is able to provide for this example maximum percentage errors of 5% or lower (and of 1% or lower in the cases of 11 and 33 layers) for $S=30$ or larger.

Quantitative results are presented in Table 2 for various plate cases, considering a number of layers equal to 11 and 10 collocation points for each in-plane parametric direction. Different number of layers (i.e., 3 and 33) are instead investigated in Appendix A. Increasing length-to-thickness ratios, namely 20, 30, 40, and 50 are considered and the maximum relative error results is reported for a reference point located at $x_{1}=x_{2}=0.25L$. Also different degrees of approximation are investigated. Given these results, we conclude that using an out-of-plane degree of approximation equal to 4 leads to a sufficiently accurate stress state. Furthermore the out-of-plane stress profile reconstruction shows a remarkable improvement for increasing values of number of layers and slenderness parameter $S$.

5 Conclusions

In this paper we present a new approach to simulate laminated plates characterized by a symmetric distribution of plies. This technique combines a 3D collocation isogeometric analysis with a post-processing step procedure based on equilibrium equations. Since we adopt a single-element appoach, to take into account variation through the plate thickness of the material properties, we average the constitutive behaviour of each layer considering an homogeneized response. Following this simple approach, we showed that acceptable results can be obtained only in terms of displacements and in-plain stresses. Therefore, we propose to perform a post-processing step which requires the shape functions to be highly continuous. This continuity demand is fully granted by typical IGA shape functions. After the post-processing correction is applied, good results are recovered also in terms of out-of-plane stresses, even for very coarse meshes. The post-processing stress-recovery technique is only based on the integration through the thickness of equilibrium equations, and all the required components can be easily computed differentiating the displacement solution. Several numerical tests are carried out to test the sensitivity of the proposed technique to different length-to-thickness ratios and number of layers. Regardless of the number of layers, the method gives better results the thinner the composites are. Multiple numbers of alternated layers and sequence of stacks (both even and odd) have been studied in our applications. Neverthless only tests which consider an odd number of layers or an odd disposition of an even number of stacks show good results as expected because considering a homogenized response of the material is effective only for symmetric distributions of plies. Further research topics currently under investigation consist in the extension of this approach to more complex problems involving curved geometries and large deformations.

Acknowledgments

This work was partially supported by Fondazione Cariplo – Regione Lombardia through the project "Verso nuovi strumenti di simulazione super veloci ed accurati basati sull’analisi isogeometrica", within the program RST – rafforzamento.

P. Antolin was partially supported by the European Research council through the H2020 ERC Advanced Grant 2015 n.694515 CHANGE.

Appendix A

Results in terms of maximum relative error considering a plate with a number of layers equal to 3 and 33 are herein presented for a reference point located at $x_{1}=x_{2}=0.25L$. Increasing length-to-thickness ratios, namely 20, 30, 40, and 50 are investigated for different degrees of approximations (i.e., $p=q=6$ and $r=4$, and $p=q=r=6$), using 10 collocation points for each in-plane parametric direction and one element through-the-thickness.