1. Introduction
In modern engineering, sandwich structures are extensively applied to construct lightweight laminated structures allowing for advanced mechanical solutions in naval, automotive, aerospace, and aeronautic applications [
1]. Among structural laminated solutions, sandwich structures stand out as a special type of laminate, which are characterised by two thin high-rigid face sheets (e.g., fibre-reinforced composite or aluminium) enclosing a soft thick core. Generally, its core is manufactured with low-stiffness and low-density materials. Nevertheless, being confined by two rigid materials allows the core to increase the bending rigidity of the structural set while keeping the overall weight of the sandwich structure low. Several factors influence the structural mechanical behaviour of sandwich structures, such as the microstructure of the core material; the thickness ratio of core and face sheets; the fibre volume ratio; and the material (and its mechanical orientation) of the face sheets [
1,
2]. Due to its mechanical complexity, the development of numerical methods capable of efficiently and accurately predicting the structural behaviour of sandwich structures is a never-ending task.
In this work, due to its versatility, polyurethane foam (PUF) has been selected as the core material. It is possible to find in the literature several research studies on PUF [
3], which is a widely used material in the polymer family due to its exceptional properties, such as lightweight, eco-friendliness, low density, shock absorption, processability, and elasticity. These foams have a cellular structure with polyhedral cells that fill three-dimensional space, which, combined with various manufacturing techniques, make it possible to design PUF with specific mechanical properties to be applied to particular applications [
4]. It is possible to produce PUF in multiple physical states and densities to meet specific needs, making it cost-effective for diverse uses. PUFs are commonly used as core materials in sandwich structures integrating aerospace, aeronautic, naval, automotive, and construction structures and mechanisms [
3].
Most commonly, sandwich panels and laminates are analysed using equivalent single-layer (ESL) theories [
5], in which the laminate thickness is modelled using a transverse deformation theory. Thus, the problem domain can be represented with a 2D domain. This dimensional simplification makes it possible to reduce the computational cost of the analysis when compared with a full 3D deformation elasticity theory for 3D discretisations. However, 3D solutions allow for more realistic physical simulations and predictions than 2D ESL theories, particularly for thick plates and shells [
6].
Thus, several analytical solutions using 3D elasticity theory have been developed to investigate the bending response of sandwich plates. In the pioneer research work of Pagano [
7], exact 3D solutions were derived for the stress analysis of simply supported rectangular sandwich plates. Later, Zenkour [
8] applied 3D elasticity equations to obtain the bending analytical solutions of rectangular multilayer plates under transversal distributed loading conditions. For sandwich plates with a functionally graded core, Kashtalyan and Menshykova [
9] proposed a 3D exact elasticity solution to predict their bending behaviour for sinusoidally distributed transverse loads. Woodward and Kashtalyan further developed the 3D exact elasticity solution to predict the bending response of sandwich plates under localised transverse loads [
10] and many other relevant transversal loading scenarios [
11].
Discretising with full detail the cellular material forming the sandwich core will lead to numerical analyses with high computational costs. Thus, multiscale homogenisation techniques are generally applied to reduce the time of the structural analyses, making it possible to obtain fair approximated solutions [
12]. Homogenisation techniques provide a multiscale analysis by assuming the existence of multiple spatial scales in materials and structures. Commonly, the analysis of heterogeneous materials has relied on effective properties derived from homogenising the response at microscopic scales, which are then transported to a macroscale analysis. Within the framework of small deformations and linear elasticity, numerous analytical models have been developed to predict the homogenised constitutive response of heterogeneous materials at the macroscopic level, incorporating the characteristics of their microstructure [
12]. These analytical models are based on the Hill–Mandel condition for homogeneity [
13]. This principle states that, when dealing with materials with vastly different microscopic and macroscopic length scales, the volume-averaged strain energy within a representative volume element (RVE) can be computed as the product of the volume-averaged stress and strain fields within that same RVE. The Hill–Mandel condition demonstrates the energy equivalency between the homogeneous and the heterogeneous material [
13,
14]. Micromechanical analysis of two-phase materials often utilises discretisation-based approaches to predict the overall response from their microstructure [
15].
In Computational Mechanics, the finite element method (FEM) is the most popular discretisation technique. Since its origin, the FEM has been successfully applied to several engineering fields [
16]. Nevertheless, the computational mechanics research community employs continuous efforts to discover and develop new discretisation techniques possessing higher efficiency and accuracy. Today, meshless methods (or meshfree methods) are alternative discretisation techniques, which are capable of replacing the FEM in most applications [
17,
18]. In meshless methods, the solid domain is discretised with an unstructured nodal set instead of the standard element mesh of FEM [
19,
20]. With the FEM, the nodal connectivity is obtained by means of the element concept. Differently, in meshless methods, such connectivity is achieved using the influence domain concept [
21].
The diffuse element method (DEM), proposed by Nayroles et al. as an generalisation of the FEM, was the first developed mature meshless method (or meshfree method) [
22] using moving least square (MLS) approximants. Later, Belytschko et al. slightly enhanced the DEM and extended it to elasticity problems [
23]. By introducing a background integration mesh based on the Gauss–Legendre quadrature scheme and applying the Lagrange multipliers to impose the boundary conditions, Belytschko et al. called their DEM version the Element-Free Galerkin Method (EFGM) [
23], becoming one of the most popular meshless method ever developed. After, in the same decade, other meshless methods appeared, such as the Reproducing Kernel Particle Method (RKPM) [
24] or the Meshless Local Petrov–Galerkin Method [
25]. By using approximant shape functions and higher nodal connectivity, these meshless methods are capable of delivering smoother variable fields and more accurate solutions [
21]. Although all these meshless methods were successfully applied to solid and fluid mechanics, their approximant shape functions lack the delta Kronecker property, hindering the direct imposition of natural and essential boundary conditions. In these approximant meshless methods, the most usual technique to impose the boundary conditions involes applying the Lagrange multipliers, which increases the size of the system of equations and, consequently, the overall cost of the analysis [
21]. Compared to the FEM, this is a significant disadvantage.
Thus, the research community started to focus their efforts in the development of meshless methods with interpolating shape functions. Proposed by Sukumar et al., the Natural Element Method (NEM) was one of the first interpolant meshless methods [
26] and a truly meshless method. With the NEM, the natural and essential boundary conditions could be exactly imposed with the same techniques applied to the FEM, reducing the overall cost of the meshless analysis. Afterwards, other ingenious interpolant meshless methods were developed, such as the Point Interpolation Method (PIM) [
27], the Radial Point Interpolation Method (RPIM) [
28], the Meshless Finite Element Method (MFEM) [
29], and the Natural Radial Element Method (NREM) [
30]. In the field of interpolant meshless methods, the RPIM is certainly the most popular technique [
21].
Later, combining the connectivity technique of the NEM with the interpolation functions of the RPIM, Dinis et al. developed a truly meshless method: the Natural Neighbour Radial Point Interpolation Method (NNRPIM) [
31]. The NNRPIM only requires the nodal discretisation of the problem’s domain. Afterwards, using that spatial information, it autonomously distributes the integration points and establishes nodal connectivity. There is no need for an independent background integration mesh, as in the EFGM or RPIM.
In the literature, it is possible to find several research works addressing the combination of meshless methods with ESL deformation theories for the analysis of sandwich structures [
32]. However, meshless formulations using a 3D deformation theory are not so common [
6,
33,
34]. Regarding multiscale homogenisation techniques combined with meshless methods, Rodrigues et al. extended the RPIM [
35] and the NNRPIM [
36] to the multiscale analysis of laminated composites, and Wang et al. developed a multiscale approach for the computational modelling of the mechanical behaviour of carbon nanotube-reinforced cement composites [
37].
Regarding the main novelty of this study, for the first time, this work combines the RPIM with a multiscale homogenisation technique to predict the macroscale behaviour of sandwich structures. This research aims to compare the performance and the efficiency of the RPIM with the FEM (computational cost, accuracy, stress field distribution) and to explore the meshless multiscale analysis of sandwich structures with a homogeneous or a functionally graded foam core. Thus, firstly, the homogenised mechanical properties of the cellular units are estimated and correlated with the foam density. Then, at a macroscale, distinct sandwich structures are modelled by varying the density distributions along the sandwich thickness. The obtained results are compared with FEM solutions. This work is a preliminary work on the structural analysis of sandwich structures using the RPIM, which is intended to understand if the RPIM is sufficiently accurate and robust in linear elastostatic applications of sandwich structures. If the accuracy and robustness of the RPIM is verified, it can be extended to more demanding applications, such as material and geometrical nonlinearities. Beyond the elastic limit, sandwich structures start to show elastoplastic effects (and damage) in the core and a consequent crushing. In computational analysis, such nonlinear behaviour will produce severe mesh distortions, which the FEM cannot handle efficiently. In opposition, meshless methods are capable of dealing with mesh distortions without requiring any remeshing.
This document is divided in four sections. In
Section 1, the state-of-the-art of sandwich structures are presented and their corresponding mathematical formulation, multiscale homogenisation approaches, and usual dicretization techniques. In
Section 2, the meshless method mathematical formulation for elasticity is presented, as well as the adopted homogenisation technique. Then, in
Section 3, the numerical results are presented and discussed. First, the obtained cellular foam homogenised mechanical properties with respect to their apparent density are presented. Afterwards, macroscale benchmark examples of sandwich structures are presented and their solutions discussed. In the end, the main conclusions and remarks are presented in
Section 4.
4. Conclusions
In this work, a full multiscale study involving sandwich beams with a polyurethane foam (PUF) core has been presented. First, the bulk PUF was modified through the inclusion of circular holes, making it possible to create PUFs with distinct volume fractions. Then, the homogenised mechanical properties of the PUF, with respect its volume fraction, were obtained using a homogenisation technique combined with finite element methods (FEMs) and four versions of the radial point interpolation method (RPIM): the classical RPIM (CRPIM) and the modified RPIM with 4, 9, and 16 nodes inside the influence domain (MRPIM4, MRPIM9, and MRPIM16, respectively). The obtained results show that RPIM formulations are capable of delivering solutions close to high-order FEM formulations, such as quadratic triangular elements (FEM-6n). However, since background triangular integration cells were considered for all RPIM formulations, to achieve variable fields with the same level of smoothness of FEM-6n, RPIM formulations require increasing the number of integration points per triangular integration cell, leading to a very high computational cost and, consequently, decreasing the overall RPIM numerical efficiency. Next, after calculating the homogenised mechanical properties of PUFs with respect their volume fractions, those mechanical properties were applied to large-scale problems: a sandwich cantilever beam with a homogeneous PUF core and a sandwich cantilever beam with an approximated functionally graded PUF core. In these macroscale examples, differently from the microscale examples, background quadrilateral integration cells were used for all RPIM formulations. Such modification made it possible to approximate significantly the RPIM and FEM results, particularly the results obtained with the MRPIM4 formulation. Regarding the vertical displacement field of the cantilever beam with a uniform core, it was observed that the solutions obtained with the CRPIM and MRPIM16 formulations had an 8% difference when compared with the FEM solution. When the comparison was made between the FEM and MRPIM4, the difference dropped to 1%. Regarding the cantilever beam with a functionally graded core, comparing the FEM solution with the CRPIM and MRPIM16 solutions (for the vertical displacement field), a 4% difference was found, and upon comparing the FEM and MRPIM4 solutions, the difference became lower: 0.04%. For the normal stress , similar differences were obtained. It was found that near the domain edge, RPIM formulations using influence domains with a large number of nodes led to local lower stress values. However, stress distributions obtained with all the distinct formulations studied tended to agree along the thickness of the beam. The results obtained for the macroscale examples consistently showed that the MRPIM4 is capable of producing results very close with the FEM-4n. Such findings are very interesting, since the MRPIM4 is the most efficient RPIM formulation, as well as the one with the lower computational cost. It was also possible to observe that all RPIM formulation possessed a higher computational cost when compared with ther FEM formulations. For instance, the MRPIM4 was four times slower than the FEM, and the CRPIM and MRPIM16 were almost ten times slower than the FEM. Such results can be explained by the lack of programming optimisation of the RPIM codes, heavy preprocessing routines, and the large nodal connectivity of RPIM formulations. The results documented in this work make it possible to verify that RPIM formulation, in its classical or modified version, is a solid alternative to the FEM. The results also show that if background triangular integration cells are being used, the number of integration points per cell should be increased, leading to a not very interesting RPIM version. Therefore, for the RPIM, background quadrilateral integration cells are much more suitable. Being a preliminary study on RPIM formulations combined with multiscale analysis, its main purpose was to assess the accuracy and robustness of RPIM in this context, paving the way for its potential application in more complex scenarios involving material and geometric nonlinearities. The results here documented will make it possible to apply the obtained mechanical properties to dynamic, buckling, large deformations, and contact analyses of sandwich structures with PUF cores. The information contained in this document make it possible to design/select PUF cores with circular voids suited for specific problems. The homogenised microscale results make it possible to correlate the density of the PUF with circular voids with macroscale mechanical properties. This output can be applied to construct more efficient structural components (maximising its stiffness/weight ratio), such as airplane wings and fuselage components, wind turbine blades, automobile impact parts (to absorb impact loads), and building construction materials, amongst many other applications.