High-order Shape Functions Research Articles

A quasi-3D SinZZ model-driven multi-field Chebyshev FEM for nonlinear vibration control in multilayer multiferroic composite plates

This paper presents an advanced numerical method for modeling and mitigating nonlinear vibrations in thin multilayered fiber-reinforced multiferroic composite plates, addressing complex multi-physical interactions. The proposed methodology leverages a multi-physical coupling Chebyshev finite element formulation, utilizing high-order shape functions derived from Chebyshev polynomials. By integrating the strengths of spectral element methods and Legendre spectral finite element methods, this approach effectively overcomes challenges such as shear locking and spurious zero energy modes while ensuring high convergence rates in multi-physical problems. A quasi-3D refined model, incorporating the Murakami zig-zag model and sinusoidal shear deformation theory, is employed to accurately capture the nonlinear Von Kármán strain–displacement relationship and the magneto-electro-elastic coupling in multilayer structures with thickness-dependent material properties. To suppress nonlinear vibrations, the study utilizes a closed-loop multiphysical Chebyshev finite element model for time-domain analysis of viscoelastically damped systems, employing the Golla–Hughes–McTavish model. The results underscore the significant influence of multiferroic properties and the strategic distribution of ferroelectric fibers within the substrate on the dynamic behavior of the plate. The numerical validation, supported by rigorous verification, demonstrates the robustness of the proposed method in effectively simulating and controlling multilayer fiber-reinforced multiferroic composite plates. Additionally, this research highlights the potential for significant vibration reduction through semi-active damping mechanisms, offering valuable insights for practical applications in industries where precision and stability are critical.

A semi-analytical spectral element model for guided wave propagation in composite laminated conical shells

The conical shell is a typical non-uniform waveguide with its circumferential radius varies linearly. This work aims to provide physical insights for the composite laminated conical shell in terms of waves. A semi-analytical spectral element model is presented, which adopts the Fourier series in the circumferential direction and the higher-order shape functions in the axial direction. The first-order shear deformation theory and Hamilton principle are employed to derive the energy functions of the conical shell, and they are discretized using the spectral elements to achieve a weak-form variational problem. The dispersion equation and the group velocity expression are derived via the piecewise approximation concept and Bloch's periodic theorem. The accuracy of the developed spectral element is verified by comparison with published results. The deformation evolution mechanisms of the composite conical shell are revealed, and the wave propagation characteristics of uniform and ply drop-off composite conical shells are investigated. The results demonstrated that wavenumber distribution at 0–5 kHz is spatially dependent, and wavenumber becomes cut off and bifurcate near the cone vertex. When the frequency exceeds 2 kHz, the deformation of conical shell changes from rigid to flexible, and the structure tends to be a uniform waveguide. Compared to the geometric curvature effect, the asymmetric lay-up leads to more significant mode coupling.

Boolean finite cell method for multi-material problems including local enrichment of the Ansatz space

AbstractThe Finite Cell Method (FCM) allows for an efficient and accurate simulation of complex geometries by utilizing an unfitted discretization based on rectangular elements equipped with higher-order shape functions. Since the mesh is not aligned to the geometric features, cut elements arise that are intersected by domain boundaries or internal material interfaces. Hence, for an accurate simulation of multi-material problems, several challenges have to be solved to handle cut elements. On the one hand, special integration schemes have to be used for computing the discontinuous integrands and on the other hand, the weak discontinuity of the displacement field along the material interfaces has to be captured accurately. While for the first issue, a space-tree decomposition is often employed, the latter issue can be solved by utilizing a local enrichment approach, adopted from the extended finite element method. In our contribution, a novel integration scheme for multi-material problems is introduced that, based on the B-FCM formulation for porous media, originally proposed by Abedian and Düster (Comput Mech 59(5):877–886, 2017), extends the standard space-tree decomposition by Boolean operations yielding a significantly reduced computational effort. The proposed multi-material B-FCM approach is combined with the local enrichment technique and tested for several problems involving material interfaces in 2D and 3D. The results show that the number of integration points and the computational time can be reduced by a significant amount, while maintaining the same accuracy as the standard FCM.

Analysis of the sound insulation performance of periodic wall structures by a virtual acoustic laboratory

The propagation of acoustic waves in periodic structures, also known as phononic crystals or PCs, is prohibited in certain frequency ranges, which are referred to as the frequency band-gaps. The existence, the location and the width of the frequency band-gaps are mainly determined by the geometrical parameters and the material properties of the PCs. In this work, a virtual acoustic laboratory based on a fully coupled fluid-structure interaction (FSI) model is developed to determine the sound insulation capacity of one-dimensional (1D) and two-dimensional (2D) periodic walls. The FSI model is discretized using the frequency-domain spectral element method (FDSEM), which is an advanced finite element method (FEM) using special high-order shape functions. Following the guidelines of the ISO10140, the setup of the developed FSI model allows us to take into account the essential physical phenomena, especially the interaction of the wall structure with the fluid domains (air). The FSI model based on the FDSEM increases the computational efficiency and accuracy in comparison with the standard FEM. Several numerical examples will be presented and discussed to show that the designed periodic walls in certain frequency ranges within the band-gaps may exhibit much better sound insulation capabilities than monolithic walls.

A nodal immersed finite element-finite difference method

The immersed finite element-finite difference (IFED) method is a computational approach to modeling interactions between a fluid and an immersed structure. The IFED method uses a finite element (FE) method to approximate the stresses, forces, and structural deformations on a structural mesh and a finite difference (FD) method to approximate the momentum and enforce the incompressibility of the entire fluid-structure system on a Cartesian grid. The fundamental approach used by this method follows the immersed boundary framework for modeling fluid-structure interaction (FSI), in which a force spreading operator prolongs structural forces to a Cartesian grid, and a velocity interpolation operator restricts a velocity field defined on that grid back onto the structural mesh. With an FE structural mechanics framework, force spreading first requires that the force itself be projected onto the finite element space. Similarly, velocity interpolation requires projecting velocity data onto the FE basis functions. Consequently, evaluating either coupling operator requires solving a matrix equation at every time step. Mass lumping, in which the projection matrices are replaced by diagonal approximations, has the potential to accelerate this method considerably. This paper provides both numerical and computational analyses of the effects of this replacement for evaluating the force projection and for the IFED coupling operators. Constructing the coupling operators also requires determining the locations on the structure mesh where the forces and velocities are sampled. Here we show that sampling the forces and velocities at the nodes of the structural mesh is equivalent to using lumped mass matrices in the IFED coupling operators. A key theoretical result of our analysis is that if both of these approaches are used together, the IFED method permits the use of lumped mass matrices derived from nodal quadrature rules for any standard interpolatory element. This is different from standard FE methods, which require specialized treatments to accommodate mass lumping with higher-order shape functions. Our theoretical results are confirmed by numerical benchmarks, including standard solid mechanics tests and examination of a dynamic model of a bioprosthetic heart valve.

A 3D finite-strain beam model for thermo-mechanical deformations of 2D shape memory alloys in 3D space

In this paper, a robust constitutive model adapted to a 3D beam theory in the finite strain regime is proposed to simulate thermo-mechanical behaviors of shape memory alloys (SMAs) under triaxial normal-shear loadings. For the first time, four internal parameters are considered to capture martensitic transformation, orientation/reorientation of martensite variants, normal-shear deformation coupling, and asymmetric/anisotropic strain generation in polycrystalline SMAs. Solution algorithms of the model are presented for stress- and strain-control cases based on the elastic-predictor inelastic-corrector return mapping process. To develop the Cauchy-Green deformation tensor, a 3D exact displacement field is proposed based on the centroid movements and rotations of the cross-section. A finite element formulation along with the Newton-Raphson and Riks techniques is also established to trace the non-linear equilibrium path of 3D SMA beam structures in the large defamation regime. Higher-order shape functions are also considered to avoid shear and membrane locking issues. To explore and demonstrate the capabilities of the proposed model, the results of the triaxial and 3D SMA constitutive models are compared with existing experimental results on uniaxial tension, and combined tension-torsion tests. Afterwards, the results of 3D solid element and proposed 3D beam element are compared with experimental results of deep stretch SMA helical springs under three loading/unloading cycles. An excellent qualitative and quantitative correlation is observed between numerical and experimental results verifying the model accuracy and the solution procedure. The combination of the proposed triaxial SMA constitutive model with 3D beam theory leads to a significant reduction in computational time and storage requirements without any accuracy loss. Extra simulations are then conducted for a simplified model of a 3D biomedical vascular SMA stent under various loading/unloading conditions. The simulations reveal that the new 3D finite-strain SMA-beam element is well appropriate for modeling of 3D spatial lattice-based stents. Due to simplicity and accuracy, the 3D SMA beam model could serve in future efforts for a computationally efficient modeling of large defamations of complicated 3D SMA structures with 2D members (e.g., SMA lattices). • A robust constitutive model is developed to simulate behaviors of SMAs under triaxial normal-shears loadings. • The solution algorithms combined is extended to trace the non-linear equilibrium path of 3D SMA beam structures. • The results of present 3D Beam model has been compared with existing experimental results as well as 3D solid elements. • A good qualitative and quantitative correlation is observed between numerical and experimental results. • The proposed triaxial SMA beam model leads to a significant reduction in computational time without loss inaccuracy.

Octree-based integration scheme with merged sub-cells for the finite cell method: Application to non-linear problems in 3D

Fictitious domain methods, such as the Finite Cell Method (FCM), allow for an efficient and accurate simulation of complex geometries by utilizing higher-order shape functions and an unfitted discretization based on rectangular elements. Since the mesh does not conform to the geometry, cut elements arise that are intersected by domain boundaries. For optimal convergence rates and the efficiency of the simulation in general, special integration schemes have to be used in such elements. In this contribution, the often used, robust octree-decomposition-based integration scheme is enhanced by a novel approach reducing the computational effort when evaluating the discontinuous integrals. This is realized by introducing an additional step, in which the local integration mesh is simplified using data compression techniques leading to fewer integration domains/points. An important advantage of the proposed method is that it can be added in a modular fashion to already existing codes. While it inherits all desired properties of the octree-decomposition-based integration scheme, it significantly reduces the number of integration points and has hardly any negative effect on the simulation accuracy. In this paper, the proposed integration scheme is introduced in detail, and investigated by means of numerical examples in the context of 3D non-linear problems.

Thermal Buckling of Laminated Plates using Modified Mantari Function

Critical buckling temperature of laminated plate under thermal load varied linearly along the thickness, is developed using a higher-order shape function which depends on a parameter ‘‘m’’, which is improved to obtain results for thin and thick plates. Laminated plates’ equations of motion are obtained using virtual work principle and solved for simply supported boundary conditions. Angle and cross laminates thermal buckled mode shapes with different E1/E2 proportion, number of plies, (α2/α1) proportion, aspect ratios, are investigated. It is observed that this shape function gives thermal buckling for thin and thick plates but with m = 0.05 that agree well with other theories and linear distribution of temperature gives a rise to critical temperature approach to 50% than those caused by uniform thermal distribution.

On the Use of High-Order Shape Functions in the SAFE Method and Their Performance in Wave Propagation Problems

In this research, high-order shape functions commonly used in different finite element implementations are investigated with a special focus on their applicability in the semi-analytical finite element (SAFE) method being applied to wave propagation problems. Hierarchical shape functions (p-version of the finite element method), Lagrange polynomials defined over non-equidistant nodes (spectral element method), and non-uniform rational B-splines (isogeometric analysis) are implemented in an in-house SAFE code, along with different refinement strategies such as h-, p-, and k-refinement. Since the numerical analysis of wave propagation is computationally quite challenging, high-order shape functions and local mesh refinement techniques are required to increase the accuracy of the solution, while at the same time decreasing the computational costs. The obtained results reveal that employing a suitable high-order basis in combination with one of the mentioned mesh refinement techniques has a notable effect on the performance of the SAFE method. This point becomes especially beneficial when dealing with applications in the areas of structural health monitoring or material property identification, where a model problem has to be solved repeatedly.

Efficient simulations for the democratization of duct transmission loss analysis

As in many industries, NVH reduction in the automotive market requires car manufacturers to adapt their designs and development methods. In particular duct system design necessitates accurate methods to capture the resonance frequencies and address noise phenomena. To do so, engineers evaluate the quality of the components through numerical simulations and validate against real-life measurements. Two main challenges arise: the simulations are often ran by non-numerical experts, who might have difficulties with setting up the models. Secondly, the computation techniques, together with the model preparation, need to be efficient to improve productivity. In this paper, an automated workflow is presented to democratize numerical techniques in duct design and sound transmission loss computation. It gives the possibility for non-experts to easily set up and solve their models with few input parameters. Furthermore, effective numerical techniques based on high-order shape functions with a-priori error estimator are used to optimize performance and accuracy: FEMAO and BEMAO (Finite/Boundary Element Methods with Adaptive Order). Two models are exposed: a simple duct line with three Helmholtz resonators, and a more complex intake system. The numerical results of the different methods are compared, as well as measurements.

Eliminating finite-grid instabilities in gyrokinetic particle-in-cell simulations

Finite-grid (aliasing) instabilities are known to place severe limitations on momentum-conserving particle-in-cell (PIC) methods applied to models including charge separation effects by requiring the resolution of the Debye length. Gyrokinetic models, on the other hand, generally enforce quasi-neutrality thereby removing the Debye length analytically. Recent studies with momentum-conserving PIC applied to gyrokinetic models, however, show that a manifestation of this instability exists in certain physical parameter regimes for arbitrary spatial resolution. In this paper, we show that a simple reformulation of the discrete equations, making use of a co-located discretization of the continuity equation, eliminates this instability for all practical purposes. We perform numerical dispersion analyses for both the original and reformulated schemes, including the effects of finite mean parallel velocity and finite β. We demonstrate that the reformulated scheme is numerically stable for stationary plasmas at any spatial resolution and for plasmas in which the electron mean parallel velocity is smaller than the electron thermal velocity (generally ensured by ambipolarity in plasmas of interest). This reformulation may be particularly useful for codes with complicated meshes, where high-order shape functions or energy-conserving schemes are difficult to implement. The analysis presented here may also be useful for explaining the success of previously considered methods to remove particle instabilities, such as the split-weight scheme.

Three dimensional extended finite element simulation of cracked functionally graded pipe and pipe bend

In this work, extended finite element method (XFEM) is used to simulate the fracture mechanics problems of part-through semi-elliptical circumferential crack in functionally graded pipe and pipe bend. Part-through semi-elliptical circumferential crack present in functionally graded pipe and pipe bend at different location has been subjected to internal pressure. Higher order shape function is used to model the crack with level set functions. The inner surface of the functionally graded pipe/pipe bend contains 100% steel alloy and outer surface contains 100% alumina ceramic. The material property gradation is followed by exponential law from inner surface to outer surface of the pipe and pipe bend in radial direction. Domain-based interaction integral approach has been used to determine the stress intensity factor at different location of the crack front.

Buckling analysis of orthotropic nanoplates based on nonlocal strain gradient theory using the higher-order finite strip method (H-FSM)

In this paper, a semi-analytical higher-order finite strip method is developed based on the nonlocal strain gradient theory (NSGT) for buckling analysis of orthotropic nanoplates. NSGT contains two material length scale parameters related to the nonlocal and strain gradient effects. To consider the effect of the strain gradient in the governing equation of the plate in the transverse direction, the higher-order polynomial shape functions (higher-order Hermitian shape functions) are utilized to assess the second derivatives, in addition to the displacement and first derivatives. Also, some numerical study is presented on the effects of different factors such as boundary conditions, nonlocal and strain gradient parameters, aspect ratio, and different types of in-plane loading for the isotropic and orthotropic rectangular nanoplate to verify the proposed formulation. In the following, a relation for the nanoplates based on the nonlocal strain gradient theory using the Navier method is extracted and presented for preliminary comparisons. According to the proposed relation, it is shown that in simply supported nanoplates, if the non-local parameter is equal to the second power of the strain gradient parameter, the responses obtained for the nanoplates are equal to the locally available responses.

Inviscid flow passing a lifting body with a higher order boundary element method

A higher order boundary element method (HOBEM) is presented for inviscid flow passing a lifting body. In addition to the boundary integral equation used for the velocity potential, similar integral equation is derived and used for the tangential velocity on the body surface. Higher order elements are used to discretize the body surface, which ensures the continuity of slope at the element nodes. The velocity potential is also expanded using higher order shape function, in which the unknown coefficients involve the tangential velocity. The expansion then ensures the continuity of the velocity at element nodes and it also allows the Kutta condition to be imposed directly through the velocity. A particular shape function is also derived and used near the trailing edge to account for the continuity of the velocity and its sharp variations there. The unknown potential and tangential velocity are then found through solving their integral equations simultaneously. Through extensive comparison of the results for a Karman-Trefftz (KT) foil, it is shown that the present HOBEM is much more accurate than the conventional BEM, in particular for the velocity and local results near the trailing edge.

An energy- and charge-conserving electrostatic implicit particle-in-cell algorithm for simulations of collisional bounded plasmas

The present work investigates the applicability of the implicit energy- and charge-conserving particle-in-cell method to collisional bounded plasmas. The possible modifications might require modeling of the charge dynamics associated with a particle impinging on a bounding surface of different materials and coupling to an external RLC circuit. The former requires no change if the lowest order shape function is used for the current density assignment, whereas the method needs to be altered when the shape function is linear and particle recombination at the surface of a perfect conductor is assumed. The proposed fix consists of adding a virtual oppositely charged particle symmetrically with respect to the conductor boundary. A newly derived numerical form of the power balance equation contains only physical terms and demonstrates that the numerical heating is controlled solely by the accuracy of the energy-conserving iteration procedure. The external circuit's contribution is incorporated into the field iterations through Kirchhoff's laws for a system consisting of the plasma reactor and the external circuit. Due to the implicit nature of the proposed algorithm modification, it does not lead to numerical instabilites, and thus is more robust than explicit analogs. The energy-conserving PIC method can be employed for numerical modeling of a wide range of technological bounded plasmas, including those with high densities, which is demonstrated with 1D simulations of a capacitively coupled plasma discharge in Cartesian geometry. It is found that the algorithm remains stable even if the iterative procedure needed for the consistent update of the particle orbits and the electric field in the framework of an implicit method is stopped after the first iteration and the lowest order shape function is used for the electric field interpolation along with the conduction current assignment. The resulting solutions prove reasonably close to the more precise counterparts obtained with the energy conservation ensured to much better accuracy and with higher order shape functions. Since this makes the computational complexity of the new algorithm comparable to that of the explicit PIC code whereas utilizing much smaller grids is possible, the new algorithm should be preferred over the latter when simulating high-density plasmas.

A Nodal Immersed Finite Element-Finite Difference Method

The immersed finite element-finite difference (IFED) method is a computational approach to modeling interactions between a fluid and an immersed structure. This method uses a finite element (FE) method to approximate the stresses and forces on a structural mesh and a finite difference (FD) method to approximate the momentum of the entire fluid-structure system on a Cartesian grid. The fundamental approach used by this method follows the immersed boundary framework for modeling fluid-structure interaction (FSI), in which a force spreading operator prolongs structural forces to a Cartesian grid, and a velocity interpolation operator restricts a velocity field defined on that grid back onto the structural mesh. Force spreading and velocity interpolation both require projecting data onto the finite element space. Consequently, evaluating either coupling operator requires solving a matrix equation at every time step. Mass lumping, in which the projection matrices are replaced by diagonal approximations, has the potential to accelerate this method considerably. Constructing the coupling operators also requires determining the locations on the structure mesh where the forces and velocities are sampled. Here we show that sampling the forces and velocities at the nodes of the structural mesh is equivalent to using lumped mass matrices in the coupling operators. A key theoretical result of our analysis is that if both of these approaches are used together, the IFED method permits the use of lumped mass matrices derived from nodal quadrature rules for any standard interpolatory element. This is different from standard FE methods, which require specialized treatments to accommodate mass lumping with higher-order shape functions. Our theoretical results are confirmed by numerical benchmarks, including standard solid mechanics tests and examination of a dynamic model of a bioprosthetic heart valve.

An hp-version adaptive finite element algorithm for eigensolutions of moderately thick circular cylindrical shells via error homogenisation and higher-order interpolation

PurposeThis study presents a novel hp-version adaptive finite element method (FEM) to investigate the high-precision eigensolutions of the free vibration of moderately thick circular cylindrical shells, involving the issues of variable geometrical factors, such as the thickness, circumferential wave number, radius and length.Design/methodology/approachAn hp-version adaptive finite element (FE) algorithm is proposed for determining the eigensolutions of the free vibration of moderately thick circular cylindrical shells via error homogenisation and higher-order interpolation. This algorithm first develops the established h-version mesh refinement method for detecting the non-uniform distributed optimised meshes, where the error estimation and element subdivision approaches based on the superconvergent patch recovery displacement method are introduced to obtain high-precision solutions. The errors in the vibration mode solutions in the global space domain are homogenised and approximately the same. Subsequently, on the refined meshes, the algorithm uses higher-order shape functions for the interpolation of trial displacement functions to reduce the errors quickly, until the solution meets a pre-specified error tolerance condition. In this algorithm, the non-uniform mesh generation and higher-order interpolation of shape functions are suitable for addressing the problem of complex frequencies and modes caused by variable structural geometries.FindingsNumerical results are presented for moderately thick circular cylindrical shells with different geometrical factors (circumferential wave number, thickness-to-radius ratio, thickness-to-length ratio) to demonstrate the effectiveness, accuracy and reliability of the proposed method. The hp-version refinement uses fewer optimised meshes than h-version mesh refinement, and only one-step interpolation of the higher-order shape function yields the eigensolutions satisfying the accuracy requirement.Originality/valueThe proposed combination of methodologies provides a complete hp-version adaptive FEM for analysing the free vibration of moderately thick circular cylindrical shells. This algorithm can be extended to general eigenproblems and geometric forms of structures to solve for the frequency and mode quickly and efficiently.

An adaptive phase-field model based on bilinear elements for tensile-compressive-shear fracture

An efficient adaptive phase-field method based on bilinear elements for tensile-compressive-shear fracture is developed. On the one hand, refined meshes are needed near the crack path to obtain accurate results, and the computational efficiency is very low if a prior mesh refinement is applied, especially when the crack paths are unknown. On the other hand, conventional phase-field models are not applicable to tensile-compressive-shear fracture under complex stress states at present. In this paper, an adaptive scheme is proposed to improve the computational efficiency, in which a bilinear multi-node element is adopted to avoid using high order quadrature and shape functions when new nodes are inserted in elements, and a new multi-node triangular element is outlined to expand the scope of application of the proposed adaptive method. Another important aspect of the contribution is the development of the phase-field model for tensile-compressive-shear fracture under complex stress states, in which a universal fracture criterion is embedded. To determine the optimal parameters in the adaptive scheme and demonstrate the advantages of the proposed bilinear adaptive method, a series of numerical examples are performed for sensitivity analysis. Comparison with experimental results is also conducted to validate the proposed phase-field model for tensile-compressive-shear fracture under complex stress states.

Higher Order Multiscale Finite Element Method for Heat Transfer Modeling.

In this paper, we present a new approach to model the steady-state heat transfer in heterogeneous materials. The multiscale finite element method (MsFEM) is improved and used to solve this problem. MsFEM is a fast and flexible method for upscaling. Its numerical efficiency is based on the natural parallelization of the main computations and their further simplifications due to the numerical nature of the problem. The approach does not require the distinct separation of scales, which makes its applicability to the numerical modeling of the composites very broad. Our novelty relies on modifications to the standard higher-order shape functions, which are then applied to the steady-state heat transfer problem. To the best of our knowledge, MsFEM (based on the special shape function assessment) has not been previously used for an approximation order higher than p = 2, with the hierarchical shape functions applied and non-periodic domains, in this problem. Some numerical results are presented and compared with the standard direct finite-element solutions. The first test shows the performance of higher-order MsFEM for the asphalt concrete sample which is subject to heating. The second test is the challenging problem of metal foam analysis. The thermal conductivity of air and aluminum differ by several orders of magnitude, which is typically very difficult for the upscaling methods. A very good agreement between our upscaled and reference results was observed, together with a significant reduction in the number of degrees of freedom. The error analysis and the p-convergence of the method are also presented. The latter is studied in terms of both the number of degrees of freedom and the computational time.

A cubic spline layerwise spectral finite element for robust stress predictions in laminated composite and sandwich strips

A new theoretical and numerical tool, comprising of a novel cubic Hermite spline layerwise theory and a spectral strip finite element with integration points collocated to its nodes and high-order polynomial shape functions, is evaluated regarding its capability to predict intralaminar and interlaminar stresses in laminated and sandwich composite strips. Stress predictions in various laminated and thick sandwich beams are presented and compared with exact solutions, high-fidelity plane-strain FE models and reported zig-zag layerwise models. Emphasis is placed on the capability of the present model to predict interlaminar stresses on the free-edges, at ply interfaces, in resin-rich layers and in laminates with artificial delaminations. All numerical results illustrate the capabilities of the present computational tool to provide accurate and robust interlaminar stress predictions implementing coarse discretization, thus ensuring the validity and rapid convergence of the combined higher-order layerwise laminate theory and spectral finite element model.