Next Article in Journal
Revealing the Corrosion Resistance Mechanism of Plain Carbon Steel Micro-Alloyed by La in Simulated Industrial Atmosphere
Next Article in Special Issue
Recent Advancements in Fabrication of Metal Matrix Composites: A Systematic Review
Previous Article in Journal
Evaluation of the Influence of the Tool Set Overhang on the Tool Wear and Surface Quality in the Process of Finish Turning of the Inconel 718 Alloy
Previous Article in Special Issue
Effects of Rubber Core on the Mechanical Behaviour of the Carbon–Aramid Composite Materials Subjected to Low-Velocity Impact Loading Considering Water Absorption
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Multiscale Analysis of Sandwich Beams with Polyurethane Foam Core: A Comparative Study of Finite Element Methods and Radial Point Interpolation Method

1
Department of Mechanical Engineering, ISEP, Polytechnic of Porto, Rua Dr. António Bernardino de Almeida, n. 431, 4200-072 Porto, Portugal
2
Institute of Science and Innovation in Mechanical and Industrial Engineering (INEGI), Rua Dr. Roberto Frias, 4200-465 Porto, Portugal
Materials 2024, 17(18), 4466; https://doi.org/10.3390/ma17184466
Submission received: 23 May 2024 / Revised: 8 July 2024 / Accepted: 31 July 2024 / Published: 11 September 2024

Abstract

:
This study presents a comprehensive multiscale analysis of sandwich beams with a polyurethane foam (PUF) core, delivering a numerical comparison between finite element methods (FEMs) and a meshless method: the radial point interpolation method (RPIM). This work aims to combine RPIM with homogenisation techniques for multiscale analysis, being divided in two phases. In the first phase, bulk PUF material was modified by incorporating circular holes to create PUFs with varying volume fractions. Then, using a homogenisation technique coupled with FEM and four versions of RPIM, the homogenised mechanical properties of distinct PUF with different volume fractions were determined. It was observed that RPIM formulations, with higher-order integration schemes, are capable of approximating the solution and field smoothness of high-order FEM formulations. However, seeking a comparable field smoothness represents prohibitive computational costs for RPIM formulations. In a second phase, the obtained homogenised mechanical properties were applied to large-scale sandwich beam problems with homogeneous and approximately functionally graded cores, showing RPIM’s capability to closely approximate FEM results. The analysis of stress distributions along the thickness of the beam highlighted RPIM’s tendency to yield lower stress values near domain edges, albeit with convergence towards agreement among different formulations. It was found that RPIM formulations with lower nodal connectivity are very efficient, balancing computational cost and accuracy. Overall, this study shows RPIM’s viability as an alternative to FEM for addressing practical elasticity applications.

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.

2. Numerical Method

2.1. Meshless Methods

Meshless methods employ a well-established discretisation approach distinct from mesh-based methods [20,21]. Unlike predefined meshes with element connectivity, for domain discretisation, meshless methods rely on a set of nodes, which can be distributed either regularly or irregularly. The accuracy of the meshless solution hinges on the nodal density, which is similar to how mesh density impacts accuracy in mesh-dependent methods using weak formulations for elasticity [21]. Generally, as the nodal density increases, the corresponding solution becomes more accurate. Furthermore, for problems with specific geometric features or concentrated loads that significantly influence stress distribution, a refined nodal density in those regions becomes crucial to improve the meshless method’s accuracy. Following the nodal discretisation, a set of background integration points is defined across the domain, which is known as the integration mesh.
Unlike finite element methods (FEMs) that rely on predefined element connectivity, meshless methods require an alternative approach to establish relationships between nodes and integration points. This is achieved through the concept of the influence domain [21]. In the influence domain approach, each integration point, denoted by x I , identifies a set of nodes within a predefined radius. These nodes will be used to construct the interpolation function of its corresponding integration point, x I , forming its influence domain. This concept offers a simple technique to establish the nodal connectivity in meshless methods, and it is employed by various meshless formulations (e.g., the RPIM, EFGM, MLPG, and RKPM). The two main approaches for defining influence domains are fixed or variable radial search. Within a fixed radial search, a constant radius for all integration points is used. While simple to implement, it can lead to inconsistencies near boundaries. Integration points close to physical boundaries may have fewer nodes within their influence domain compared to those located within the domain, which can affect the accuracy of the solution. Alternatively, a variable radial search aims to ensure a consistent number of nodes influencing each integration point. This is achieved by adjusting for each integration point x I the searching radius, which will make it possible to maintain the same number of nodes within the influence domain of all integration points. Although this technique enforces that all integration points will construct shape functions with the same order, it can be computationally more expensive due to the dynamic radius calculation.
After establishing the nodal connectivity of each integration point, shape functions functions are constructed, making it possible to approximate/interpolate the variable field at an integration point x I with
u ( x I ) = i = 1 n φ i ( x I ) · u ( x i ) = Φ ( x I ) T · u s = { φ 1 ( x I ) , φ 2 ( x I ) , , φ n ( x I ) } · u ( x 1 ) u ( x 2 ) u ( x n )
where the number of nodes within the influence domain of the integration point x I is defined as n, the field variable components of each node x i inside the influence domain of x I are contained in u s , and φ i ( x I ) is the ith component of the shape function constructed for x I . These functions, and their partial derivatives, are then incorporated into the chosen governing equation (strong or weak form) to establish a system of algebraic equations. Solving this system yields the solution variable fields (pressure, displacement, velocity, etc.) across the domain. Details regarding the RPIM (nodal connectivity, background integration mesh, shape functions, and discretisation) and the elastostatic formulations can be found in Appendix A.

2.2. Material Homogenisation Technique

The representative volume element (RVE) concept, introduced by Hill [13], offers an efficient tool for the multiscale analysis of materials. An RVE represents a statistically relevant subvolume of the material, capturing the essential microstructural features that govern its overall behaviour. A key aspect is ensuring statistical representativeness: the RVE must encompass a sufficient sampling of all the microstructural features within the material, even for materials with varying constituent volume fractions. The appropriate RVE size remains a critical parameter. It must be large enough to capture the relevant heterogeneities but remain significantly smaller than the macroscopic domain of interest. Appendix B introduces the RVE concept to present a micromechanical approach for determining the effective elastic properties of materials. This approach makes it possible to bridge scales by numerically connecting the microscopic behaviour of the heterogeneous constituents and the homogenised effective behaviour observed at the macroscopic level.

3. Numerical Results

In this section, microscale and macroscale examples are analysed, and the obtained results are discussed. First, the homogenised mechanical properties of a polyurethane foam (PUF) sample with circular voids were obtained using the homogenisation procedure described in Section 2.2. Then, using the homogenised mechanical properties, a macroscale application was analysed: a sandwich cantilever beam with aluminium face sheet layers and PUF core. Several PUF densities were analysed, making it possible to observe the effects of considering a constant density for the PUF core or a functionally graded density.
In order to compare the obtained RPIM results, finite element formulations were applied [38]. Thus, in this work, the following FEM formulations were considered:
  • FEM-3n: constant strain linear triangular elements (three-node elements), with one integration point per element.
  • FEM-6n: quadratic triangular elements (six-node elements), with full integration.
  • FEM-4n: Lagrangian four-node elements, with full integration.
Regarding the RPIM formulations, in this work, four RPIM formulations were studied:
  • CRPIM: a classic RPIM considering 16 nodes inside the influence domain, shape parameters c = 1.42 and p = 1.03 , and a linear polynomial basis.
  • MRPIM16: a modified version of the RPIM considering 16 nodes inside the influence domain, shape parameters c = 0.0001 and p = 0.9999 , and a constant polynomial basis.
  • MRPIM9: a modified version of the RPIM considering nine nodes inside the influence domain, shape parameters c = 0.0001 and p = 0.9999 , and a constant polynomial basis.
  • MRPIM4: a modified version of the RPIM considering four nodes inside the influence domain, shape parameters c = 0.0001 and p = 0.9999 , and a constant polynomial basis.
Regardless the RPIM formulation, the following integration schemes were assumed:
  • If the background integration mesh is made with triangular cells, one integration point is considered for each triangular integration cell, located at the triangular cell centre, and with an integration weight ω I equal to its volume ( ω I = A I · h I , with A I being the triangle area and h I being the element thickness).
  • If quadrilaterals are used to build the background integration mesh, the Gauss–Legendre quadrature scheme is adopted, assuming 2 × 2 integration points per quadrilateral cell (as shown in Figure A1).
All the finite element and meshless methods codes were completely written by the author using Matlab R2017bTM. The nodal and element meshes were built with Simcenter Femap Student Edition SoftwareTM. All the analyses were performed with a standard laptop computer, with an Intel© CoreTM i5-3230M CPU (Samsung, Suwon, Republic of Korea), 2.60 GHz, and 8 GB of RAM.

3.1. Microscale Analysis

In this section, the homogenised mechanical properties of PUF with microscale cylindrical holes are investigated. Thus, assuming the plane strain deformation theory and a unitary thickness, a 2D parametric RVE was considered, Figure 1, with L = D = 10 mm, x c = { L / 2 , D / 2 } , and radius R of the central hole defined as
R ( v f ) = 1 v f π
in which v f represents the volume fraction of the RVE,
v f = v m v d
with v m being the volume occupied by effective PUF mass and v d being the volume of the RVE: v d = L · D · 1 mm3.
For the bulk PUF ( v f = 1 ), the isotropic elastic mechanical properties documented in the literature [39] were considered: E = 171.43 MPa, G = 57.810 MPa, and ν = 0.30 .

3.1.1. Convergence Study

Since this work proposes a discrete meshless method for a multiscale analysis, it is necessary to verify the convergence of the meshless technique at the microscale. Therefore, a volume fraction of v f = 0.75 was considered, corresponding to a radius R = 2.821 mm. Next, five distinct meshes were constructed, as shown in Figure 2.
For each mesh indicated in Figure 2, the procedure described in Section 2.2 was applied to obtain the elastic components of the homogenised constitutive material matrix, C i j , and the homogenised material properties of the RVE: the Young modulus in direction O x and O y E x and E y , respectively—and the in-plane elastic shear modulus and Poisson’s ratio— G x y and ν x y , respectively. The obtained results are presented in Figure 3.
Although the final converged value was not the same for the all the techniques, the results of Figure 3 show that all the discretisation techniques converged. Both of the FEM formulations converged to the same value, with FEM-6n being the one showing a higher convergence rate. Similarly, regarding the elastic mechanical property, all the RPIM formulations converged to the same value as well, with the MRPIM4 being the one showing the highest convergence rate and the CRPIM being the formulation with the lower one. It is possible to notice that the FEM formulations presented an upper-bound convergence path, and the RPIM formulations showed a lower-bound path. Some formulations presented very close results. Thus, their convergence curves apparently overlapped. As Figure 3 shows, the CRPIM and MRPIM16 overlapped their solutions in all presented results, and in Figure 3b, the solution of all the RPIM formulations overlapped, being visible in only one line. In order to show how close the mentioned solutions truly are and allow future comparisons, the numerical values of the convergence study are shown Table 1 and Table 2.
Concerning the homogenised elastic mechanical properties obtained for the most refined mesh, in Table 1, it is possible to observe that the difference between the RPIM and FEM formulations is around 3.5 % for the Young modulus E x and E y , 1.5 % for the Young’s modulus E z , 6.5 % for the elastic shear modulus G x y , and 1.2 % for the Poisson’s ratio. Since the FEM and RPIM techniques are fundamentally different in their formulations, it is expected to find such differences in the obtained results.
Observing the elastic constitutive constants obtained for the most refined mesh, in Table 2, the output can be seen to follow the same trend. The results of the FEM formulations are very similar between each other, and the RPIMs’ results are all very close to each other as well. Comparing the FEM results with the RPIM solutions, it is possible to observe the following differences: 3.3 % for C 11 and C 22 , 3.0 % for C 12 , 1.5 % for C 33 , and 6.5 % for C 44 .
In order to understand the stress distribution obtained with each formulation, the von Mises equivalent stress field ( σ e f = 3 · J 2 ) obtained with each formulation is presented in Figure 4 and Figure 5.
Observing Figure 4 and Figure 5, it is possible to observe that the FEM-6n formulation delivered a very smooth stress distribution. Such a result was expected, because its shape functions possess quadratic terms. Then, the CRPIM and MRPIM16 presented stress distributions very similar with the FEM-3n formulation, which was not expected.
Generally, the literature shows that meshless methods are capable of delivering smoother stress fields [21,23]. However, the obtained RPIM effective stress fields, presented in Figure 4 and Figure 5, apparently contradict such a notion. The stress distributions of the CRPIM and MRPIM16 formulations present the same granulated distribution observed in the FEM-3n formulation. The results became even worst when a lower number of nodes was considered inside the influence domains, as for the MRPIM9 and MRPIM4 formulations (Figure 5d–i).
The reason behind this lack of smoothness comes from the adopted integration scheme. In order to deliver a faster and more flexible RPIM formulation, when triangular integration cells are used, only one integration point is inserted inside the triangular integration cell. Although such a procedure does not significantly affect the accuracy of the variable fields, it affects the smoothness distribution of such fields. When the integration order was increased, using for example four integration points inside each integration cell (following the Gauss–Legendre integration scheme described in [21]), the smoothness of the stress fields significantly increased, as shown in Figure 6.
Comparing the results obtained with the CRPIM and the CRPIM* version with an higher integration order (comparing Figure 4g–i with Figure 5a–c), it is possible to visualise that the stress distribution improved significantly, then becoming more smooth and continuous (the granulated distributions have disappeared).
Nevertheless, although the stress field (and all the other variable fields as well) became smoother, the homogenised values did not change. Table 3 shows the results obtained with the CRPIM* version with an higher integration order. Comparing such results with the solutions presented in Table 1 and Table 2, it is possible to observe that the homogenised values were not affected by the integration scheme.
The main disadvantage of using the higher integration scheme (for any of the RPIM formulations) is the computational cost. Using one integration point inside each triangular integration cell makes it possible to perform analyses about 10 times faster than using four integration points, which becomes relevant when nonlinear analyses are being considered. Thus, the next examples will continue to consider only one integration point inside each triangular integration cell.
Regarding the overall computational cost of the FEM and RPIM formulations, Figure 7 is presented. The displayed data concers the complete computational cost, including the nodal connectivity, construction of the background integration mesh, and construction of the shape functions, thus establishing the global system of equations and obtaining all the variable fields (displacements, strains, and stresses). As Figure 7 shows, the FEM presented a very low computational cost. The fastest RPIM formulation was the MRPIM4. Even so, it was four times slower than the FEM. The CRPIM and MRPIM16 presented the highest computational costs, being almost 10 times slower than the FEM. Nevertheless, it is relevant to clarify that the code written by the author was not written for speed and to achieve computational efficiency. Several programming improvements could be applied to the RPIM code to significantly reduce the overall computational cost. Most of the RPIMs’ computational burden is due to the preprossessing phase (nodal connectivity, construction of the background integration mesh, and construction of the shape functions). In the FEM, the element mesh (containing the information of which node belongs to each element) is already available. In the RPIM, it is necessary to first find the influence domains. Then, in the FEM, the shape functions and their partial derivatives are already available (using isoparametric transformation), so it is not necessary to invert any geometric matrix for each integration point as in the RPIM. Moreover, the FEMs’ shape functions possess three or six nodes, and the RPIMs’ shape functions can possess 16 nodes, which will delay the shape function construction step, as well as the stiffness matrix construction phase. Also, the background integration mesh is already established in the FEM (again using isoparametric transformation), and the RPIM formulation requires the explicit construction of such a background mesh. All these reasons significantly increase the computational cost of RPIM formulations.

3.1.2. Homogenised Material Properties

Next, since the aim was to investigate the influence of the PUF volume fraction ( v f ) on the final homogenised mechanical properties, several decreasing volume fractions were considered. Thus, according to Equation (2) and considering the volume fractions presented in Table 4, distinct corresponding radii for the circular hole were obtained, as seen in Table 4.
The results of the convergence test, Figure 3, show that a nodal mesh with a nodal density corresponding to the nodal mesh of about 4000 nodes for v f = 0.75 was enough, providing results very close to the theoretical converged value. Thus, the following analysed RVEs were all discretised with nodal meshes containing a proportional mesh density, as can be observed in Figure 8.
The same geometry presented in Figure 1, with L = D = 10 mm and bulk PUF mechanical properties, was considered: E = 171.43 MPa, G = 57.810 MPa, and ν = 0.30 . The homogenised mechanical properties for each volume fraction v f were once again estimated using FEM and RPIM formulations, and the obtained results are presented in Figure 9.
As Figure 9 shows, the homogenised elastic mechanical properties are dependent on the PUC volume fraction. From the graphs, it is perceptible that it is possible to overlap a polynomial curve for each homogenised elastic mechanical property as a function of the volume fraction. The results show that the FEM and RPIM formulations provided close results and followed the same trend. The most significant difference can be observed for the Poisson’s ratio, particularly for lower volume fractions. With these results, it is possible to build macroscale sandwich beam or plate models, with a PUF core, and analyse such large-scale applications using the homogenised mechanical properties of PUF with distinct volume fractions (due to the cylindrical holes).
As in the previous nodal mesh convergence study, the results obtained for the homogenised elastic mechanical properties presented in Figure 9 show that all the RPIM formulations led to very close results. They were so close that apparently all the RPIM curves overlapped each other. Such an observation can be confirmed with Table 5, in which the homogenised mechanical properties of the PUF, with distinct volume fractions, are presented for the FEM-3n, CRPIM, and MRPIM16 formulations.
Notice that for the MRPIM16 formulation, the homogenised Poisson’s ratio for v f = 0.5 was higher than ν = 0.5 . Such a result, also observed in other kind of foams [40] for sandwich plates, hinders the possibility to use a plane strain assumption in a macroscale application using such homogenised mechanical properties. However, 2D plane stress deformation and complete 3D deformation theories can be used. Thus, in order to use the same numerical starting point, in the macroscale analysis, the homogenised mechanical properties were only obtained with the FEM-3n.

3.2. Macroscale Analysis

In this section, the homogenised mechanical properties shown in Table 5 are applied to analyse a sandwich cantilever beam. In the first numerical example, homogeneous PUF cores bounded by aluminium face sheets were numerically analysed. Next, the mechanical properties of the PUF core were modified along the beam thickness, leading to an approximated functionally graded sandwich cantilever beam. In both examples, the following linear elastic isotropic mechanical properties for aluminium were considered [41]: a Young modulus of E = 69.6 GPa and a Poisson’s ratio of ν = 0.33 . The problems here presented were analysed with five distinct formulations—FEM-4n, CRPIM and MRPIM4, MRPIM9, and MRPIM16—considering quadrilateral integration cells for the construction of the background integration mesh. Regarding the discretisation, both analysed problems used a regular mesh with 30 × 120 divisions, leading to a uniform regular mesh of 3751 nodes and 3600 quadrilateral elements (or, in the case of the CRPIM and MRPIM, 3600 integration cells).

3.2.1. Sandwich Cantilever Beam

In the first example, the sandwich cantilever beam presented in Figure 10 was analysed. The beam, clamped in the left side, possesses the following dimensions: L = 4 m, D 1 = D 3 = 0.1 m, and D 2 = 0.8 m. In the O z direction, an unitary thickness ( H = 1 m) was considered. In the right-top corner of the beam, a localized load was applied: P = 10 6 N. Aiming to understand the influence on the beam’s displacement and stress fields of using PUF cores with distinct volume fractions, several analyses were performed considering for the PUF core the following six volume fractions: v f = { 0.5 , 0.6 , 0.7 , 0.8 , 0.9 , 1.0 } . Then, from Table 5, the corresponding homogenised mechanical properties (obtained using the FEM-3n) were considered.
The displacement components along O x (u) and along O y (v) obtained in point A (represented in Figure 10) are presented in Table 6 and Table 7, respectively.
It is possible to observe that the MRPIM4 consistently presented results very close compared with the FEM-4n, with both results separated by only a 1% difference. On the other hand, the CRPIM and MRPIM16 also presented results that were extremely close, showing differences below 0.01%. Comparing the CRPIM (or the MRPIM16) with the FEM-4n, the difference between solutions increased to around 8%.
Regarding the stress field, the normal stress σ x x at point B (in the top aluminium layer) and the shear stress τ x y at point C (in the PUF core at y = 0.5 m) were documented, which are both represented in Figure 10. The obtained results are presented in Table 8 and Table 9.
In accordance with the results of u A and v A already discussed, it is possible to observe that once again the MRPIM4 and FEM-4n produced results very close to each other. Also, the CRPIM and MRPIM16 obtained very similar solutions. One possible reason for this results is the similarity between the nodal connectivity of the MRPIM4 and the FEM-4n, as well as the CRPIM and the MRPIM16. Both the MRPIM and FEM-4n used only four nodes to construct their shape functions at each integration point. Moreover, due to geometric coincidences, since uniform regular nodal meshes with regular quadrilateral elements (or integration cells, in the case of the CRPIM or MRPIM formulaitons) were being used, the four nodes closer to each integration point were the same, regardless of being an FEM-4n or an MRPIM formulation. For the CRPIM and MRPIM16 formulations, the background integration mesh is the same, and the shape functions of each integration point will be constructed using exactly the same 16 nodes in both formulations. Therefore, the only difference between the CRPIM and MRPIM16 is the shape parameters c and p, explaining why both formulations produced such close results.
In order to show how the normal stress σ x x varies on the clamped edge (along the beam thickness), Figure 11, Figure 12 and Figure 13 are presented. Since the magnitude of the normal stress σ x x changes significantly from the aluminium face sheets to the PUF core, in Figure 12 presents the normal stress σ x x on the PUF core, and in Figure 11 and Figure 13, the normal stress σ x x on the top and bottom aluminium layers are shown, respectively.
It is also relevant to understand how the shear stress τ x y varies along the beam thickness on the clamped edge. Therefore, in Figure 14, Figure 15 and Figure 16, such variation is presented. Likewise, the normal stress σ x x and the magnitude of the shear stress τ x y observed in the top and bottom aluminium layers are very different from the magnitudes in the PUF core. Therefore, once again, their representation is separated: Figure 15 presents the shear stress τ x y on the PUF core, and Figure 14 and Figure 16 correspond to the shear stress τ x y on the top and bottom aluminium layers, respectively.
With Figure 11 (normal stress σ x x at the top aluminium layer), it is possible to observe that the FEM-4n and MRPIM4 produced very close solutions. Another interesting observation is how close the CRPIM and MRPIM16 solutions were with the FEM-4n and MRPIM4 solutions. The main difference between these four formulations was observed only at the ends of the top aluminium layer ( y [ 0.98 , 1.00 ] and y [ 0.90 , 0.92 ] ). This localised difference is explained by the number of nodes inside the influence domain of both the CRPIM and MRPIM16. By assuming a fixed number of 16 nodes inside the influence domain of each integration point, corner integration points (such as the ones close to point B) will construct shape functions with larger support domain radii than other inner integration points. This effect will smooth the approximated variable (decreasing its magnitude) at these peripheral integration points [21]. The same observation can be made with Figure 13 (normal stress σ x x at the bottom aluminium layer).
Regarding the PUF core, Figure 12 shows that the CRPIM and MRPIM16 followed a similar trend. Again, the FEM-4n and MRPIM4 solutions were almost overlapped. Due to the bending effect, positive values (tensile stress) for σ x x when y > 0.5 m were expected and negative values (compressive stress) when y > 0.5 m. However, such an effect was not clearly observed, particularly for the CRPIM and MRPIM16 formulations, showing σ x x < 0 for y > 0.5 m, and σ x x > 0 for y < 0.5 m, which featured the opposite of the expected effect. Such results can be explained (once again) by the number of nodes inside the influence domain. The integration points in the PUF core, near the interface between the aluminium sheet and the PUF core, have the strong influence of the strain/stress field of the aluminium sheet. Thus, since some nodes belonging to the aluminium sheet are contained within the influence domain of those integration points (that actually belong to the PUF core), the integration points on the PUF core near the interface will show a localised perturbation on the strain/stress field. Such an effect is not visible when a lower connectivity was considered: the MRPIM4. Both the CRPIM and MRPIM16 formulations showed σ x x > 0 for y > 0.5 m and σ x x < 0 for y < 0.5 m, which was the expected effect. Only at the interface, due to the vicinity of the aluminium sheet, the signal was inverted.
The normal stress σ x x level in the PUF core was very low, showing the highest magnitudes observed in the top and bottom aluminium layers.
For the shear stress τ x y , the results show that the highest values could be found, again, at the top and bottom aluminium layers (Figure 14 and Figure 16). Concerning the PUF core, Figure 15, the distribution was approximately parabolic, as expected. The shear stress τ x y maintained the trend observed for the displacement and normal stress σ x x , the FEM-4n and MRPIM4s solutions, were always very close between each other, and the CRPIM and MRPIM16 results were also very close to each other. Regarding the perturbation effect observed for the normal stress σ x x near the aluminium/PUF interface, the same was also visible here; however, it was highly attenuated.

3.2.2. Functionally Graded Sandwich Cantilever Beam

In the final macroscale example, two cantilever beams with aluminium face sheets and approximately functionally graded PUF cores were analysed. Thus, as shown in Figure 17, the PUF core was divided into eight layers, and the beam possessed the following dimensions: L = 4 m and D i = 0.1 m for i = { 1 , 2 , , 10 } . The beam was clamped on the left side, and a localised load P = 10 6 N was applied at the top-right corner.
Two distinct approximately functionally graded PUF cores were analysed (FG1 and FG2), with the mechanical properties of each layer indicated in Table 10.
Both beams were analysed with the following formulations: the FEM-4n, CRPIM, MRPIM4, MRPIM9, and MRPIM16. The O x displacements obtained in point A (represented in Figure 17) are presented in Table 11, and the corresponding O y displacements are shown in Table 12. As already observed in previous examples, the FEM-4n and MRPIM4, as well as the CRPIM and MRPIM16, produced similar results between each other. Thus, the difference in the magnitude of the u A or v A displacement components between the FEM-4n and MRPIM4 was about 0.04%, and between the CRPIM and MRPIM16, it was about 0.02%. However, the difference between the pair FEM-4n/MRPIM4 with the pair CRPIM/MRPIM16 was around 4%. These results show that all formulations produced approximated solutions.
Regarding the normal stress at point B and the shear stress at point C, both represented in Figure 17, the obtained results are shown in Table 13 and Table 14, respectively.
The results again show a similarity between the FEM-4n and MRPIM4 solutions, around 1% for σ x x and 2.5% for τ x y , and also for the solutions of the CRPIM and MRPIM16, around 0.5% for σ x x and 0.1% for τ x y . Such differences are in accordance with previous studies. Comparing the solutions of the FEM-4n with CRPIM, higher differences were found: about 10% for σ x x and only 5% for τ x y . Such a difference corroborates the previous idea that (in meshless methods with a large number of nodes inside the influence domain) integration points near a boundary (such as point B) have larger support domains, locally lowering the magnitude of the variable field (due to smoothing effects), and integration points inside the physical domain (such as point C) possess smaller support domains, allowing for sharper values.
In order to visualise the distribution of the stress field on the clamped edge, in Figure 18 and Figure 19 are shown the distribution of the normal stress σ x x and shear stress τ x y along the thickness of the beam.
In Figure 18a,b,e,f, it is possible to visualise that the CRPIM and MRPIM16 values near the edge surface of the beam ( y 0.0 and y 1.0 ) were smoothed. This is why the localised results of the CRPIM and MRPIM16, presented in Table 13 and Table 14, are different from the ones obtained with the FEM-4n and MRPIM4. Observing the normal stress distribution between y ] 0.9 , 1.0 [ m and y ] 0.0 , 0.1 [ m, it is possible to understand that all solutions were very close (with the exception of the MRPIM9). Regarding the normal stress along y [ 0.1 , 0.9 ] m, despite the similarity of pairs FEM-4n/MRPIM4 and CRPIM/MRPIM16, all solutions were acceptably close to each other.
For the shear stress τ x y , the results are in accordance with previous studies. Notice that at the free edges of the beam ( y 0.0 and y 1.0 ), in Figure 19a,b,e,f, the solutions of the CRPIM and MRPIM16 significantly differed from the solutions of the FEM-4n and MRPIM4 (due to the size of the support domain). However, as the solution moved to the PUF core, as seen in Figure 19c,d, the results started to match each other. All solutions seemed to produce similar results.
Figure 18 and Figure 19 also allow us to understand that, regardless of the volume fraction distribution along the PUF thickness, the distinct formulations were capable of providing very similar stress distributions.

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 σ x x , 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.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The original contributions presented in the study are included in the article, further inquiries can be directed to the corresponding author.

Acknowledgments

The author acknowledges the support provided by the Ministério da Ciência, Tecnologia e Ensino Superior—Fundação para a Ciência e a Tecnologia (Portugal) and LAETA under internal project UIDB/50022/2020.

Conflicts of Interest

The author declares no conflicts of interest.

Appendix A. The Radial Point Interpolation Method

Appendix A.1. Shape Functions

In order to obtain the background integration points and to allow the numerical integration of the integrodifferential equations governing the studied phenomenon, RPIM classic formulation uses the Gauss–Legendre quadrature integration scheme. Therefore, the problem domain can be divided with regular integration cells, Figure A1, and then each integration cell is transformed into a isoparametric square in which the integration points are included following the Gauss–Legendre integration scheme, as seen in Figure A1.
The classic formulation of the RPIM employs the Gauss–Legendre quadrature integration scheme for numerical integration of the governing integrodifferential equations [21,28], which facilitates the efficient evaluation of integrals within the problem domain. To achieve this, the domain can be discretised into a set of regular integration cells, as seen in Figure A1. Each cell is then mapped (transformed) into a unit square using a parametric transformation (isoparametric mapping). This transformation allows for the distribution of integration points within the isoparametric square, thus respecting the Gauss–Legendre integration scheme [21,28].
In the end, as Figure A1 illustrates, it is possible to retrieve the Cartesian coordinates of the integration points from their isoparametric coordinates with Equation (A1).
x I = x I y I = N ( ξ I , η I ) T N ( ξ I , η I ) T · x 1 T x 2 T x 3 T x 4 T = N 1 N 2 N 3 N 4 N 1 N 2 N 3 N 4 · x 1 y 1 x 2 y 2 x 3 y 3 x 4 y 4
where N ( ξ , η ) = { N 1 ( ξ , η ) , N 2 ( ξ , η ) , N 3 ( ξ , η ) , N 4 ( ξ , η ) } T = { N 1 , N 2 , N 3 , N 4 } T . The generic equation of N i ( ξ , η ) is defined as
N i ( ξ , η ) = 1 4 · ( 1 + ξ ^ i · ξ ) · ( 1 + η ^ i · η )
and assuming the following ξ ^ i and η ^ i ,
Ξ = { ξ ^ 1 , ξ ^ 2 , ξ ^ 3 , ξ ^ 4 } = { 1 , 1 , 1 , 1 } H = { η ^ 1 , η ^ 2 , η ^ 3 , η ^ 4 } = { 1 , 1 , 1 , 1 }
The ratio of between the number of integration points and nodes discretising the problem domains is a recurrent topic in the literature [21,28,42]. If the nodes are coincident with the quadrilateral vertices, 2 × 2 integration points per integration cell suffice [21]. Thus, in this work, 2 × 2 integration points are assumed for the quadrilateral integration cell, with the following isoparametric coordinates and weights:
ξ I η I ω I = 1 3 1 3 1 3 1 3 1 3 1 3 1 3 1 3 1 1 1 1
Then, the Cartesian integration weight of each integration point can be obtained with
ω ^ I = A c e l A i s o · ω I
where A c e l is the total area of the integration cell, and A i s o the total area of the isoparametric cell, which is A i s o = 2 × 2 = 4 .
Therefore, it is possible to integrate a function f ( x , y ) , defined inside an quadrilateral domain Ω h , using a Gauss–Legendre integration scheme with 2 × 2 = 4 integration points:
F = Ω h f ( x , y ) d Ω h = I = 1 4 f ( x I , y I ) · ω ^ I
where we assume that the Cartesian coordinates x I = { x I , y I } T and integration weights ω ^ I of the integration points are calculated using Equations (A1) and (A5), respectively.
Repeating this process for each integration cell, the problem’s domain is discretised with a cloud of n Q integration points, forming the background integration mesh. Since the background integration cell are coincident with the problem’s domain, it is possible to calculate the domain’s area with A Ω = I = 1 n Q ω ^ I .
As previously established, the RPIM relies on overlapping influence domains to enforce nodal connectivity. In this study, each integration point x I will identify its n closest nodes, establishing its influence domain. Studies by Wang et al. suggest that for two-dimensional (2D) problems, an appropriate number of nodes per influence domain falls within the range of n = 9 to n = 16 [28]. Maintaining a constant number of nodes within each influence domain ensures that the RPIM shape functions constructed for all integration points share a similar level of complexity.
Figure A1. Integration scheme for the RPIM using a rectangular domain and a regular background integration grid.
Figure A1. Integration scheme for the RPIM using a rectangular domain and a regular background integration grid.
Materials 17 04466 g0a1
As the name of the RPIM suggests, the RPIM shape functions are constructed using the radial point interpolating technique. In order to understand its construction, consider a 2D domain, Ω R 2 , in which a field function u ( x ) is discretised with N nodes: X = { x 1 , x 2 , , x N } R 2 . Then, for a given integration point x I R 2 x I X , its interpolated value u ( x I ) can be obtained with
u ( x I ) = { r ( x I ) T , p ( x I ) T } a b
where r ( x I ) = { r 1 ( x I ) , , r n ( x I ) } T is a radial basis function (RBF), p ( x I ) is a polynomial basis function p ( x I ) = { p 1 ( x I ) , , p m ( x I ) } T with m monomials, and a i ( x I ) and b j ( x I ) are nonconstant coefficients of r ( x I ) and p ( x I ) , respectively, which are defined as a ( x I ) = { a 1 ( x I ) , , a n ( x I ) } T and b ( x I ) = { b 1 ( x I ) , , b m ( x I ) } T .
The inclusion of the polynomial basis function confers robustness and stability to the RPIM shape functions. For instance, the inclusion of a polynomial of order k ensures C k consistency and allows the RPIM to pass the standard patch test. Nevertheless, due to the computational cost of adding high-order polynomials, commonly, only low-order polynomial basis are included, such as constant polynomial basis ( p ( x I ) = { 1 } , m = 1 ), linear polynomial basis ( p ( x I ) = { 1 , x I , y I } T , m = 3 ), or quadratic polynomial basis ( p ( x I ) = { 1 , x I , y I , x I 2 , x I · y I , y I 2 } T , m = 6 ).
The preliminary works on the RPIM [28,42] suggest for the RBF the Multiquadratic Radial Basis Function (MQ-RBF), which was initially proposed by Hardy [43]. Later, Belinha [21] proposed a slight modification to the MQ-RBF to account for the spatial dimension of the problem’s domain:
r j ( x i ) = r i j = ( d i j 2 + ( ω ^ I · c ) 2 ) p = ( x j x i ) 2 + ( y j y i ) 2 2 + ( ω ^ I · c ) 2 p
where the MQ-RBF shape parameters are defined by c and p. For 2D analyses, the original RPIM literature recommends c = 1.42 and p = 1.03 . However, posterior works have shown that there are other suitable values for c and p. In the work of Belinha [21], through the optimisation of the Kronecker delta property of the 2D RPI shape functions, it was found that c should be close to zero, but not zero, and p should be close to one, but not one [21], where the following MQ-RBF shape parameters were found to be suitable: c = 10 4 and p = 1 10 4 . The results documented in the literature [21] show that these alternative values for the shape parameters make it possible to reproduce highly complex variable fields and to achieve accurate and stable results.
In the seminal works on the RPIM [28,42], it was found that to obtain a unique solution, the following equality has to be imposed in order to allow a unique solution:
i = 1 n p j ( x i ) a i ( x i ) = 0
with j = { 1 , 2 , , m } . Consequently, a new equation matrix can be established by combining Equation (A7) with Equation (A9):
R P P T 0 a b = G a b = u s 0
where the vector of the nodal parameters u s is defined as
u s = { u 1 , u 2 , , u n } T
The RBF matrix R can be computed with
R = r 11 r 12 r 1 n r 21 r 22 r 2 n r n 1 r n 2 r n n
and the polynomial basis matrix P can be obtained with
p = p 1 ( x 1 ) p 2 ( x 1 ) p m ( x 1 ) p 1 ( x 2 ) p 2 ( x 2 ) p m ( x 2 ) p 1 ( x n ) p 2 ( x n ) p m ( x n )
Because R is a symmetric matrix, moment G is symmetric as well. Resolving Equation (A10) makes it possible to obtain
a b = G 1 · u s 0
Finally, replacing the nonconstant values, { a , b } T , of Equation (A14) in Equation (A7), the interpolation of x I is obtained:
u ( x I ) = { r ( x I ) T , p ( x I ) T } · G 1 · u s 0 = { Φ ( x I ) , Ψ ( x I ) } · u s 0
where the interpolation function of x I is represented by Φ ( x I ) :
{ Φ ( x I ) , Ψ ( x I ) } = { r ( x I ) T , p ( x I ) T } · G 1 = { ϕ 1 ( x I ) , , ϕ n ( x I ) , ψ 1 ( x I ) , , ψ m ( x I ) }
The integrodifferential equations ruling elasticity require the calculation of the partial derivatives of Φ ( x I ) , which can be computed (for general direction ξ ) with
Φ ξ ( x I ) = { r ( x I ) ξ T , p ( x I ) ξ T } · G 1
where the corresponding MQ-RBF partial derivatives are obtained with
( r i j ) ξ = 2 · p · ( d i j 2 + c 2 ) p 1 · ( ξ j ξ i )
Radial point interpolating functions satisfy the partition of unity and possess the Kronecker delta property, permitting the possibility to directly impose boundary conditions [21].

Appendix A.2. Discrete System of Equations

Within the RPIM, the virtual work principle can be applied to obtain the global system of equations. Therefore, consider a body with a solid domain Ω , bounded by a surface boundary Γ , containing natural and essential boundaries Γ t and Γ u , respectively. The solid domain has its movement constrained at Γ u , and it experiences a body force b and external forces t ¯ on Γ t . Establishing that the work produced by external forces is equal to the work produced by the internal forces, W i n t = W e x t , it is possible to write
Ω δ ε T · σ · d Ω = Ω δ u ( x I ) T · b · d Ω + Γ t δ u ( x I ) T · t ¯ · d Γ
Expanding Equation (1) makes it possible to interpolate both displacement components { u , v } simultaneously:
u ( x I ) = u ( x I ) v ( x I ) = H ( x I ) · u = ϕ 1 ( x I ) 0 ϕ n ( x I ) 0 0 ϕ 1 ( x I ) 0 ϕ n ( x I ) · u 1 v 1 u n v n
Consequently, the deformation vector can be represented as
ε ( x I ) = L · u ( x I ) = x 0 0 y y x · Φ ( x I ) · u = B ( x I ) · u = B ( x I ) · u 1 v 1 u n v n
making it possible to obtain the deformation matrix B ( x I ) :
B ( x I ) = ϕ 1 ( x I ) x 0 ϕ n ( x I ) x 0 0 ϕ 1 ( x I ) y 0 ϕ n ( x I ) y ϕ 1 ( x I ) y ϕ 1 ( x I ) x ϕ n ( x I ) y ϕ n ( x I ) x
Next, assuming Hooke’s law, the stress can be calculated at the integration point x I :
σ ( x I ) = c ( x I ) · ε ( x I ) = c ( x I ) · B ( x I ) · u
where c ( x I ) is the material constitutive matrix of integration point x I , which for a 2D problem can be defined for plane stress or plane stress assumptions [38]. Thus, considering plane stress, it can be defined as
c = E 2 1 α 1 · ν 21 2 α 1 α 1 · ν 21 0 α 1 · ν 21 1 0 0 0 α 2 · ( 1 α 1 · ν 21 2 )
and assuming plane strain, the constitutive matrix can be written as
c = E 2 α 3 · α 4 α 1 · ( 1 α 1 · ν 21 2 ) α 1 · ν 21 · α 3 0 α 1 · ν 21 · α 3 1 ν 12 2 0 0 0 α 2 · α 3 · α 4
where α 1 = E 1 / E 2 , α 2 = G 12 / E 2 , α 3 = 1 + ν 12 , and α 4 = 1 ν 12 2 · α 1 ν 21 2 . The directional mechanical properties E i , G i j , and ν i j are defined as follows: E i is the Young modulus in material direction i, and G i j and ν i j are the elastic shear modulus and Poisson’s ratio in the material distortion plane O i j , respectively.
Thus, with the substitution of stress and strain vectors in Equation (A19), the following expression is obtained:
Ω δ B ( x I ) · u T · c ( x I ) · B ( x I ) · u · d Ω = Ω δ H ( x I ) · u T · b · d Ω + Γ t δ H ( x I ) · u T · t ¯ · d Γ
which can be further simplified because small strains ( δ B ( x I ) = 0 and δ H ( x I ) = 0 ) are being considered:
Ω δ u T · B ( x I ) T · c ( x I ) · B ( x I ) · u · d Ω = Ω δ u T · H ( x I ) T · b · d Ω + Γ t δ u T · H ( x I ) T · t ¯ · d Γ
δ u T Ω B ( x I ) T · c ( x I ) · B ( x I ) · d Ω · u = δ u T Ω H ( x I ) T · b · d Ω + δ u T Γ t H ( x I ) T · t ¯ · d Γ
leading to the final system of equations,
K · u = f b + f t
which makes it possible to calculate the unknown global displacement field:
u = K 1 · f b + f t
where K is the global stiffness matrix
K = Ω B ( x I ) T · c ( x I ) · B ( x I ) · d Ω = I = 1 N Q B ( x I ) T · c ( x I ) · B ( x I ) · ω ^ I · h ( x I )
and f t and f b are the external force and global body vectors, respectively:
f t = Γ t H ( x I ) T · t ¯ · d Γ = 0 = J = 1 n Q H ( x J ) T · t ¯ · ω ^ J
f b = Ω H ( x I ) T · b · d Ω = I = 1 N Q H ( x I ) T · b · ω ^ I · h ( x I )
Notice that h ( x I ) is the thickness of the 2D solid at the position of integration point x I . Given that the RPIM is an interpolating technique, penalty methods or the direct imposition method can be applied to impose essential and natural boundary conditions.

Appendix B. Homogenisation Technique

The RVE, denoted by a volume domain Ω ˜ , is associated with a specific material point x I at the macroscale (defined by its Cartesian coordinates x I ). Within the RVE, at a microscale domain, infinitesimal points x ˜ Ω ˜ collectively represent the macroscale material point x I . A local deformation enforced at the macroscale point, x I , induces a corresponding perturbation of the RVE’s equilibrium state. Considering the motion of x represented by function ϕ , the macroscale deformed configuration of x can be represented as X = ϕ ( x ) , making it possible to define the displacement of material point x as
u ( x ) = ϕ ( x ) x = X x
leading to,
X = x + u ( x )
Thus, the deformation gradient, F ( x ) , can be defined as a function of the displacement:
F ( x ) = X x = ( x + u ( x ) ) x = I + u ( x ) x
where the second-order identity tensor is denoted as I .
If a macroscale deformation gradient at material point x (i.e., to an RVE) is applied, and assuming X ˜ as the microscale coordinates of a deformed point belonging to the RVE, a microscopic displacement field u ˜ ( X ˜ ) is obtained [44]:
u ˜ ( X ˜ ) = ( F ( X ) I ) · x ˜ + u ¯ ( X ˜ )
where ( F ( X ) I ) · x ˜ is the linear displacement field, and u ¯ ( X ˜ ) is the displacement fluctuation field.
With the equilibrium equations, and assuming a specific macroscale deformation gradient on the RVE, it is possible to calculate the RVE’s displacement field u ˜ ( x ˜ ) . Afterwards, for each RVE’s integration point, x ˜ I Ω ˜ , the corresponding RVE’s strain and stress fields are obtained using u ˜ ( x ˜ ) .
Applying the homogenisation principle and using the volume average of stresses, σ ˜ , and strains, ε ˜ , on the RVE’s volume Ω ˜ , the macroscale stress, σ ¯ , and strain, ε ¯ , at material point x I can be calculated:
σ ¯ ( x I ) = 1 Ω ˜ Ω ˜ σ ˜ ( x ˜ I ) · d Ω ˜ = 1 I = 1 n Q ˜ ω ˜ I I = 1 n Q ˜ σ ˜ ( x ˜ I ) · ω ˜ I
ε ¯ ( x I ) = 1 Ω ˜ Ω ˜ ε ˜ ( x ˜ I ) · d Ω ˜ = 1 I = 1 n Q ˜ ω ˜ I I = 1 n Q ˜ ε ˜ ( x ˜ I ) · ω ˜ I
where ω ˜ I is the microscale integration weight of the integration point x ˜ I , and n Q ˜ is the total number of the microscale integration points discretising the RVE.
Subscale modelling relies on the fundamental principle of energy equivalence embodied by the Hill–Mandel principle [13]. This principle ensures consistency between the macroscopic and microscopic energy states within a model. In simpler terms, for a model to be energetically valid, the total deformation energy observed at the macroscopic level must be equal to the average work performed by the stresses acting at the microscopic level. This relationship can be mathematically written as
σ ¯ ( x I ) · ε ¯ ( x I ) = 1 Ω ˜ Ω ˜ σ ˜ ( x ˜ I ) · ε ˜ ( x ˜ I ) · d Ω ˜
Thus, the virtual work expressed in Equation (A28) is also valid for the RVE, leading to the following expression for the microscale domain (neglecting body forces), which means that the virtual work expressed in Equation (A28) is valid for the RVE. Thus, considering a microscopic displacement δ u ˜ and neglecting body forces, applying the virtual work principle to the RVE leads to
δ u ˜ T · Ω ˜ B ( x ˜ I ) T · c ( x ˜ I ) · B ( x ˜ I ) · u ˜ · d Ω ˜ = δ u ˜ T · Γ ˜ H ( x ˜ I ) T · t ˜ · d Γ ˜
where δ u ˜ is the virtual microscopic displacement, and t ˜ is the traction force applied in the RVE’s boundary surface Γ ˜ . Considering ε ˜ ( x ˜ I ) = B ( x ˜ I ) · u ˜ makes it possible to rewrite Equation (A41) as
Ω ˜ B ( x ˜ I ) T · c ( x ˜ I ) · ε ˜ ( x ˜ I ) · d Ω ˜ = Γ ˜ H ( x ˜ I ) T · t ˜ · d Γ ˜ = f ˜
This expression makes it possible to impose the RVE macroscale deformation fields and obtain corresponding equivalent forces.
To determine the homogenised elastic properties of the 2D RVE, plane strain simplification will be considered. This approach is valid when one dimension (in this case, the O z direction) is significantly larger than the other two ( O x or z directions). Consequently, the RVE is assumed to have a theoretical infinite length in the O z direction. For instance, any voids within the RVE, such as a circular void, would be extended infinitely along O z , resulting in an infinitely long cylinder. The plain strain deformation theory assumes that ε z z = 0 , γ x z = 0 and γ y z = 0 , with σ z z 0 . The generalised Hooke’s law for the plane strain condition can be written using Voigt notation as [35]
σ x x σ y y σ z z τ x y = c [ 4 × 4 ] · ε x x ε y y 0 γ x y = C 11 C 12 C 13 0 C 21 C 22 C 23 0 C 31 C 32 C 33 0 0 0 0 C 44 · ε x x ε y y 0 γ x y
or,
σ x x σ y y σ z z τ x y = 1 ν 23 ν 32 E 2 E 3 Δ ν 21 + ν 23 ν 31 E 2 E 3 Δ ν 31 + ν 21 ν 32 E 2 E 3 Δ 0 ν 12 + ν 32 ν 13 E 1 E 3 Δ 1 ν 13 ν 31 E 1 E 3 Δ ν 32 + ν 12 ν 31 E 1 E 3 Δ 0 ν 31 + ν 12 ν 23 E 1 E 2 Δ ν 23 + ν 21 ν 13 E 1 E 2 Δ 1 ν 12 ν 21 E 1 E 2 Δ 0 0 0 0 G 12 · ε x x ε y y 0 γ x y
in which parameter Δ can be defined as
Δ = 1 ν 23 ν 32 ν 13 ν 31 ν 12 ν 21 2 ν 32 ν 13 ν 21 E 1 E 2 E 3
Parameters E 1 and E 2 are the elasticity modulus in plane material directions 1 and 2, and E 3 is the elasticity modulus along direction 3 (i.e., O z ). The Poisson ratio, representing the ratio between the deformation observed in direction j when a force is applied in direction i, is defined with ν i j , and the G 12 stands for the plane shear elasticity modulus. Although Equation (A44) is useful to obtain the stress field, its matrix size [ 4 × 4 ] is not compatible with the matrix size of the deformation matrix B , hindering the computational application of (A30). Thus, at the microscale or at the macroscale, the stiffness matrix is constructed using the following material constitutive matrix [ 3 × 3 ] :
c = C 11 C 12 0 C 21 C 22 0 0 0 C 44
In this work, the homogenised elastic properties of the RVE are obtained with the following procedure:
  • First, because only 2D examples are analysed in this work, three deformation fields are defined:
    ε 100 = 1 0 0 ; ε 010 = 0 1 0 ; ε 001 = 0 0 1
  • Then, three load cases are obtained by substituting each of the deformation fields of Equation (A43) into Equation (A42):
    Ω ˜ B ( x ˜ I ) T · c ( x ˜ I ) · ε 100 · d Ω ˜ = f ˜ 100 Ω ˜ B ( x ˜ I ) T · c ( x ˜ I ) · ε 010 · d Ω ˜ = f ˜ 010 Ω ˜ B ( x ˜ I ) T · c ( x ˜ I ) · ε 001 · d Ω ˜ = f ˜ 001
  • Afterwards, the following essential boundary conditions are imposed at the RVE for equivalent force obtained with Equation (A48):
    for x ˜ i = 0 0 , u ˜ ¯ = 0 and v ˜ ¯ = 0 for x ˜ i = L ˜ 0 , v ˜ ¯ = 0
    where u ˜ ¯ and v ˜ ¯ are the imposed displacements along the O x ˜ and O y ˜ directions, respectively, and L ˜ is the length of the RVE.
  • Next, with Equation (A30), the system of equations is solved, and the RVE’s displacement field for each considered load case is obtained: u ˜ 100 , u ˜ 010 and u ˜ 001
  • Consequently, with Equations (A21) and (A23) the RVE’s strain fields are obtained—( ε ˜ 100 , ε ˜ 010 , and ε ˜ 001 )—as well as the stress fields—( σ ˜ 100 , σ ˜ 010 and σ ˜ 001 ).
  • Thereafter, it is possible to obtain the RVE’s homogenised stress ( σ ¯ 100 , σ ¯ 010 , and σ ¯ 001 ) and strain ( ε ¯ 100 , ε ¯ 010 , and ε ¯ 001 ) with Equations (A38) and (A39), respectively.
  • Then, with Equation (A43), a system of 3 × 4 equations is established:
    σ ¯ 100 = c [ 4 × 4 ] · ε ¯ 100 σ x x [ 100 ] σ y y [ 100 ] σ z z [ 100 ] τ x y [ 100 ] = C 11 C 12 C 13 0 C 21 C 22 C 23 0 C 31 C 32 C 33 0 0 0 0 C 44 · ε x x [ 100 ] ε y y [ 100 ] 0 γ x y [ 100 ] σ ¯ 010 = c [ 4 × 4 ] · ε ¯ 010 σ x x [ 010 ] σ y y [ 010 ] σ z z [ 010 ] τ x y [ 010 ] = C 11 C 12 C 13 0 C 21 C 22 C 23 0 C 31 C 32 C 33 0 0 0 0 C 44 · ε x x [ 010 ] ε y y [ 010 ] 0 γ x y [ 010 ] σ ¯ 001 = c [ 4 × 4 ] · ε ¯ 001 σ x x [ 001 ] σ y y [ 001 ] σ z z [ 001 ] τ x y [ 001 ] = C 11 C 12 C 13 0 C 21 C 22 C 23 0 C 31 C 32 C 33 0 0 0 0 C 44 · ε x x [ 001 ] ε y y [ 001 ] 0 γ x y [ 001 ]
    Notice that 3 of these 12 equations are linear dependent and lead to the same result: τ = C 44 · γ x y —making it possible to obtain G 12 = C 44 . The remaining nine components C i j , with { i , j } = 1 , 2 , 3 , are obtained solving the other 9 equations.
  • Finally, after calculating C i j , Equation (A44) can be applied to obtain the elastic mechanical properties:
    E 1 = ζ C 22 C 33 ( C 23 ) 2 E 2 = ζ C 11 C 33 ( C 13 ) 2 E 3 = ζ C 11 C 22 ( C 12 ) 2 G 12 = C 44 , ν 12 = C 12 C 33 C 13 C 23 C 22 C 33 ( C 23 ) 2 ν 13 = C 13 C 22 C 12 C 23 C 22 C 33 ( C 23 ) 2 ν 23 = C 23 C 11 C 12 C 13 C 11 C 33 ( C 13 ) 2
    where
    ζ = C 11 C 22 C 33 + 2 C 23 C 13 C 12 C 11 ( C 23 ) 2 C 22 ( C 13 ) 2 C 33 ( C 12 ) 2
The homogenised mechanical properties will be determined under 2D plane strain assumptions. This necessitates an initial estimation of the elasticity modulus in the direction perpendicular to the plane strain (i.e., the O z direction), because term C 33 in Equation (A50) is always multiplied by zero due to the plane strain condition. Therefore, the rule of mixtures (ROMs) is employed to estimate this out-of-plane modulus [35]
E 3 R O M = i = 1 n v f ( i ) E 3 ( i )
where the number of material fractions is defined by n. In this work, only two fractions are considered: material and absence of material (i.e., the void). The volume fraction of fraction i is represented by v f ( i ) , and the Young’s modulus in material direction 3 of fraction i is denoted as E 3 ( i ) .

References

  1. Castanie, B.; Bouvet, C.; Ginot, M. Review of composite sandwich structure in aeronautic applications. Compos. Part C Open Access 2020, 1, 100004. [Google Scholar] [CrossRef]
  2. Garg, A.; Belarbi, M.; Chalak, H.; Chakrabarti, A. A review of the analysis of sandwich FGM structures. Compos. Struct. 2021, 258, 113427. [Google Scholar] [CrossRef]
  3. Khan, T.; Acar, V.; Aydin, M.; Hülagü, B.; Akbulut, H.; Seydibeyoğlu, M. A review on recent advances in sandwich structures based on polyurethane foam cores. Polym. Compos. 2020, 41, 2355–2400. [Google Scholar] [CrossRef]
  4. Jamil, A.; Guan, Z.; Cantwell, W. The static and dynamic response of CFRP tube reinforced polyurethane. Compos. Struct. 2020, 161, 85–92. [Google Scholar] [CrossRef]
  5. Caliri, M.; Ferreira, A.; Tita, V. A review on plate and shell theories for laminated and sandwich structures highlighting the Finite Element Method. Compos. Struct. 2016, 156, 63–77. [Google Scholar] [CrossRef]
  6. Vaghefi, R. Thermo-elastoplastic analysis of functionally graded sandwich plates using a three-dimensional meshless model. Compos. Struct. 2020, 242, 112144. [Google Scholar] [CrossRef]
  7. Pagano, N. Exact Solutions for Rectangular Bidirectional Composites and Sandwich Plates. J. Compos. Mater. 1970, 4, 20–34. [Google Scholar] [CrossRef]
  8. Zenkour, A. Three-dimensional Elasticity Solution for Uniformly Loaded Cross-ply Laminates and Sandwich Plates. J. Sandw. Struct. Mater. 2007, 9, 213–238. [Google Scholar] [CrossRef]
  9. Kashtalyan, M.; Menshykova, M. Three-dimensional elasticity solution for sandwich panels with a functionally graded core. Compos. Struct. 2009, 87, 36–43. [Google Scholar] [CrossRef]
  10. Woodward, B.; Kashtalyan, M. Bending response of sandwich panels with graded core: 3D elasticity analysis. Mech. Adv. Mater. Struct. 2010, 17, 586–594. [Google Scholar] [CrossRef]
  11. Woodward, B.; Kashtalyan, M. 3D elasticity analysis of sandwich panels with graded core under distributed and concentrated loadings. Int. J. Mech. Sci. 2011, 53, 872–885. [Google Scholar] [CrossRef]
  12. Ghosh, S. Micromechanical Analysis and Multi-Scale Modeling Using the Voronoi Cell Finite Element Method, 1st ed.; CRC Press: Boca Raton, FL, USA, 2011. [Google Scholar] [CrossRef]
  13. Hill, R. Elastic properties of reinforced solids: Some theoretical principles. J. Mech. Phys. Solids 1963, 11, 357–372. [Google Scholar] [CrossRef]
  14. Hill, R. The essential structure of constitutive laws for metal composites and polycrystals. J. Mech. Phys. Solids 1967, 15, 79–95. [Google Scholar] [CrossRef]
  15. Kheyabani, A.; Massarwa, E.; Kefal, A. Multiscale structural analysis of thick sandwich structures using parametric HFGMC micromechanics and isogeometric plate formulation based on refined zigzag theory. Compos. Struct. 2022, 297, 115988. [Google Scholar] [CrossRef]
  16. Rao, S. The Finite Element Method in Engineering, 6th ed.; Butterworth-Heinemann: Oxford, UK, 2018. [Google Scholar] [CrossRef]
  17. Patel, V.; Rachchh, N. Meshless method—Review on recent developments. Mater. Today Proc. 2020, 26, 1598–1603. [Google Scholar] [CrossRef]
  18. Zhang, M.; Abidin, A.; Tan, C. State-of-the-art review on meshless methods in the application of crack problems. Theor. Appl. Fract. Mech. 2024, 131, 104348. [Google Scholar] [CrossRef]
  19. Gu, Y. Meshfree methods and their comparisons. Int. J. Comput. Methods 2005, 2, 477–515. [Google Scholar] [CrossRef]
  20. Nguyen, V.; Rabczuk, T.; Bordas, S.; Duflot, M. Meshless methods: A review and computer implementation aspects. Math. Comput. Simul. 2008, 79, 763–813. [Google Scholar] [CrossRef]
  21. Belinha, J. Meshless Methods in Biomechanics—Bone Tissue Remodelling Analysis, 1st ed.; Springer: Berlin/Heidelberg, Germany, 2014. [Google Scholar] [CrossRef]
  22. Nayroles, B.; Touzot, G.; Villon, P. Generalizing the finite element method: Diffuse approximation and diffuse elements. Comput. Mech. 1992, 10, 307–318. [Google Scholar] [CrossRef]
  23. Belytschko, T.; Lu, Y.; Gu, L. Element-free Galerkin methods. Int. J. Numer. Methods Eng. 1994, 37, 229–256. [Google Scholar] [CrossRef]
  24. Liu, W.; Jun, S.; Zhang, F. Reproducing kernel particle methods. Int. J. Numer. Methods Fluids 1995, 20, 1081–1106. [Google Scholar] [CrossRef]
  25. Atluri, S.; Zhu, T. A new Meshless Local Petrov-Galerkin (MLPG) approach in computational mechanics. Comput. Mech. 1998, 22, 117–127. [Google Scholar] [CrossRef]
  26. Sukumar, N.; Moran, B.; Belytschko, T. The natural element method in solid mechanics. Int. J. Numer. Methods Eng. 1998, 43, 839–887. [Google Scholar] [CrossRef]
  27. Liu, G.R.; Gu, Y.T. A point interpolation method for two-dimensional solids. Int. J. Numer. Methods Eng. 2001, 50, 937–951. [Google Scholar] [CrossRef]
  28. Wang, J.G.; Liu, G.R. A point interpolation meshless method based on radial basis functions. Int. J. Numer. Methods Eng. 2002, 54, 1623–1648. [Google Scholar] [CrossRef]
  29. Idelsohn, S.; Onate, E.; Calvo, N.; Del Pin, F. The meshless finite element method. Int. J. Numer. Methods Eng. 2003, 58, 893–912. [Google Scholar] [CrossRef]
  30. Belinha, J.; Dinis, L.; Natal Jorge, R. The natural radial element method. Int. J. Numer. Methods Eng. 2013, 93, 1286–1313. [Google Scholar] [CrossRef]
  31. Dinis, L.; Natal Jorge, R.; Belinha, J. Analysis of 3D solids using the natural neighbour radial point interpolation method. Comput. Methods Appl. Mech. Eng. 2007, 196, 2009–2028. [Google Scholar] [CrossRef]
  32. Neves, A.; Ferreira, A.; Carrera, E.; Cinefra, M.; Roque, C.; Jorge, R.; Soares, C. Static, free vibration and buckling analysis of isotropic and sandwich functionally graded plates using a quasi-3D higher-order shear deformation theory and a meshless technique. Compos. Part B 2013, 44, 657–674. [Google Scholar] [CrossRef]
  33. Dinis, L.; Natal Jorge, R.; Belinha, J. A 3D shell-like approach using a natural neighbour meshless method: Isotropic and orthotropic thin structures. Compos. Struct. 2010, 92, 1132–1142. [Google Scholar] [CrossRef]
  34. Dinis, L.; Natal Jorge, R.; Belinha, J. Composite Laminated Plates: A 3D Natural Neighbor Radial Point Interpolation Method Approach. J. Sandw. Struct. Mater. 2010, 12, 119–138. [Google Scholar] [CrossRef]
  35. Rodrigues, D.; Belinha, J.; Pires, F.; Dinis, L.; Natal Jorge, R. Material homogenization technique for composites: A meshless formulation. Sci. Technol. Mater. 2018, 30, 2603–6363. [Google Scholar] [CrossRef]
  36. Rodrigues, D.; Belinha, J.; Pires, F.; Dinis, L.; Natal Jorge, R. Homogenization technique for heterogeneous composite materials using meshless methods. Eng. Anal. Bound. Elem. 2018, 92, 73–89. [Google Scholar] [CrossRef]
  37. Wang, J.; Zhang, L.; Liew, K. A multiscale modeling of CNT-reinforced cement composites. Comput. Methods Appl. Mech. Eng. 2016, 309, 411–433. [Google Scholar] [CrossRef]
  38. Zienkiewicz, O.; Taylor, R. The Finite Element Method, 5th ed.; Butterworth-Heinemann: Oxford, UK, 2000. [Google Scholar]
  39. Dhaliwal, G.; Newaz, G. Flexural Response of Degraded Polyurethane Foam Core Sandwich Beam with Initial Crack between Facesheet and Core. Materials 2020, 13, 5399. [Google Scholar] [CrossRef]
  40. Luo, Y.; Yuan, K.; Shen, L.; Liu, J. Sandwich panel with in-plane honeycombs in different Poisson’s ratio under low to medium impact loads. Rev. Adv. Mater. Sci. 2021, 60, 145–157. [Google Scholar] [CrossRef]
  41. Comte, C.; von Stebut, J. Microprobe-type measurement of Young’s modulus and Poisson coefficient by means of depth sensing indentation and acoustic microscopy. Surf. Coatings Technol. 2020, 154, 42. [Google Scholar] [CrossRef]
  42. Wang, J.G.; Liu, G.R. On the optimal shape parameters of radial basis functions used for 2-D meshless methods. Comput. Methods Appl. Mech. Eng. 2002, 191, 2611–2630. [Google Scholar] [CrossRef]
  43. Hardy, R. Theory and applications of the multiquadric-biharmonic method 20 years of discovery 1968–1988. Comput. Math. Appl. 1990, 19, 163–208. [Google Scholar] [CrossRef]
  44. Reis, F.; Pires, F. A mortar based approach for the enforcement of periodic boundary conditions on arbitrarily generated meshes. Comput. Methods Appl. Mech. Eng. 2014, 274, 169–191. [Google Scholar] [CrossRef]
Figure 1. Parametric representation of the analysed RVE.
Figure 1. Parametric representation of the analysed RVE.
Materials 17 04466 g001
Figure 2. Discretisation meshes used with the FEM and RPIM analyses. (a) A total of 138 nodes and 234 triangular elements. (b) A total of483 nodes and 884 triangular elements. (c) A total of 1785 nodes and 3408 triangular elements. (d) A total of 4005 nodes and 7768 triangular elements. (e) A total of 6993 nodes and 13,664 triangular elements.
Figure 2. Discretisation meshes used with the FEM and RPIM analyses. (a) A total of 138 nodes and 234 triangular elements. (b) A total of483 nodes and 884 triangular elements. (c) A total of 1785 nodes and 3408 triangular elements. (d) A total of 4005 nodes and 7768 triangular elements. (e) A total of 6993 nodes and 13,664 triangular elements.
Materials 17 04466 g002
Figure 3. Discretisation convergence study of the following elastic mechanical properties: (a) E x , (b) E z , (c) G x y , and (d) ν x y .
Figure 3. Discretisation convergence study of the following elastic mechanical properties: (a) E x , (b) E z , (c) G x y , and (d) ν x y .
Materials 17 04466 g003
Figure 4. von Mises equivalent stress distribution for the following formulations: FEM-6n (ac), FEM-3n (df), and CRPIM (gi). The results corresponding to f ˜ 100 are shown in (a,d,g), f ˜ 010 results are shown in (b,e,h), and f ˜ 001 results are shown in (c,f,i).
Figure 4. von Mises equivalent stress distribution for the following formulations: FEM-6n (ac), FEM-3n (df), and CRPIM (gi). The results corresponding to f ˜ 100 are shown in (a,d,g), f ˜ 010 results are shown in (b,e,h), and f ˜ 001 results are shown in (c,f,i).
Materials 17 04466 g004
Figure 5. von Mises equivalent stress distribution for the following formulations: MRPIM16 (ac), MRPIM9 (df), and MRPIM4 (gi). The results corresponding to f ˜ 100 are shown in (a,d,g), f ˜ 010 results are shown in (b,e,h), and f ˜ 001 results are shown in (c,f,i).
Figure 5. von Mises equivalent stress distribution for the following formulations: MRPIM16 (ac), MRPIM9 (df), and MRPIM4 (gi). The results corresponding to f ˜ 100 are shown in (a,d,g), f ˜ 010 results are shown in (b,e,h), and f ˜ 001 results are shown in (c,f,i).
Materials 17 04466 g005
Figure 6. von Mises equivalent stress distribution obtained with CRPIM considering a higher-order integration scheme: (a) f ˜ 100 , (b) f ˜ 010 , and (c) f ˜ 001 .
Figure 6. von Mises equivalent stress distribution obtained with CRPIM considering a higher-order integration scheme: (a) f ˜ 100 , (b) f ˜ 010 , and (c) f ˜ 001 .
Materials 17 04466 g006
Figure 7. Overall computational cost of FEM and RPIM formulations.
Figure 7. Overall computational cost of FEM and RPIM formulations.
Materials 17 04466 g007
Figure 8. Discretisation meshes used with the FEM and RPIM analyses. (a) v f = 0.95 ; 5073 nodes and 9839 triangular elements. (b) v f = 0.90 ; 4807 nodes and 9321 triangular elements. (c) v f = 0.85 ; 4539 nodes and 8803 triangular elements. (d) v f = 0.80 ; 4273 nodes and 8285 triangular elements. (e) v f = 0.75 ; 4005 nodes and 7767 triangular elements. (f) v f = 0.70 ; 3739 nodes and 7249 triangular elements. (g) v f = 0.65 ; 3472 nodes and 6732 triangular elements. (h) v f = 0.60 ; 3205 nodes and 6214 triangular elements. (i) v f = 0.55 ; 2937 nodes and 5696 triangular elements. (j) v f = 0.50 ; 2671 nodes and 5178 triangular elements.
Figure 8. Discretisation meshes used with the FEM and RPIM analyses. (a) v f = 0.95 ; 5073 nodes and 9839 triangular elements. (b) v f = 0.90 ; 4807 nodes and 9321 triangular elements. (c) v f = 0.85 ; 4539 nodes and 8803 triangular elements. (d) v f = 0.80 ; 4273 nodes and 8285 triangular elements. (e) v f = 0.75 ; 4005 nodes and 7767 triangular elements. (f) v f = 0.70 ; 3739 nodes and 7249 triangular elements. (g) v f = 0.65 ; 3472 nodes and 6732 triangular elements. (h) v f = 0.60 ; 3205 nodes and 6214 triangular elements. (i) v f = 0.55 ; 2937 nodes and 5696 triangular elements. (j) v f = 0.50 ; 2671 nodes and 5178 triangular elements.
Materials 17 04466 g008
Figure 9. Influence of the volume fraction v f on the homogenised elastic mechanical properties: (a) E x , (b) E y , (c) G x y , and (d) ν x y .
Figure 9. Influence of the volume fraction v f on the homogenised elastic mechanical properties: (a) E x , (b) E y , (c) G x y , and (d) ν x y .
Materials 17 04466 g009
Figure 10. Sandwich cantilever beam with aluminium face sheets and PUF core. Three material domains were assumed: aluminium top-face sheet (3), PUF-core (2), and aluminium bottom face sheet (1).
Figure 10. Sandwich cantilever beam with aluminium face sheets and PUF core. Three material domains were assumed: aluminium top-face sheet (3), PUF-core (2), and aluminium bottom face sheet (1).
Materials 17 04466 g010
Figure 11. Variation in the normal stress σ x x along the aluminium top face sheet for distinct PUF cores. (a) v f = 0.5 . (b) v f = 0.6 . (c) v f = 0.7 . (d) v f = 0.8 . (e) v f = 0.9 . (f) v f = 1.0 .
Figure 11. Variation in the normal stress σ x x along the aluminium top face sheet for distinct PUF cores. (a) v f = 0.5 . (b) v f = 0.6 . (c) v f = 0.7 . (d) v f = 0.8 . (e) v f = 0.9 . (f) v f = 1.0 .
Materials 17 04466 g011
Figure 12. Variation in the normal stress σ x x along the PUF core for distinct PUF cores. (a) v f = 0.5 . (b) v f = 0.6 . (c) v f = 0.7 . (d) v f = 0.8 . (e) v f = 0.9 . (f) v f = 1.0 .
Figure 12. Variation in the normal stress σ x x along the PUF core for distinct PUF cores. (a) v f = 0.5 . (b) v f = 0.6 . (c) v f = 0.7 . (d) v f = 0.8 . (e) v f = 0.9 . (f) v f = 1.0 .
Materials 17 04466 g012
Figure 13. Variation in the normal stress σ x x along the aluminium bottom face sheet for distinct PUF cores. (a) v f = 0.5 . (b) v f = 0.6 . (c) v f = 0.7 . (d) v f = 0.8 . (e) v f = 0.9 . (f) v f = 1.0 .
Figure 13. Variation in the normal stress σ x x along the aluminium bottom face sheet for distinct PUF cores. (a) v f = 0.5 . (b) v f = 0.6 . (c) v f = 0.7 . (d) v f = 0.8 . (e) v f = 0.9 . (f) v f = 1.0 .
Materials 17 04466 g013
Figure 14. Variation in the shear stress τ x y along the aluminium top face sheet for distinct PUF cores. (a) v f = 0.5 . (b) v f = 0.6 . (c) v f = 0.7 . (d) v f = 0.8 . (e) v f = 0.9 . (f) v f = 1.0 .
Figure 14. Variation in the shear stress τ x y along the aluminium top face sheet for distinct PUF cores. (a) v f = 0.5 . (b) v f = 0.6 . (c) v f = 0.7 . (d) v f = 0.8 . (e) v f = 0.9 . (f) v f = 1.0 .
Materials 17 04466 g014
Figure 15. Variation in the shear stress τ x y along the PUF core for distinct PUF cores. (a) v f = 0.5 . (b) v f = 0.6 . (c) v f = 0.7 . (d) v f = 0.8 . (e) v f = 0.9 . (f) v f = 1.0 .
Figure 15. Variation in the shear stress τ x y along the PUF core for distinct PUF cores. (a) v f = 0.5 . (b) v f = 0.6 . (c) v f = 0.7 . (d) v f = 0.8 . (e) v f = 0.9 . (f) v f = 1.0 .
Materials 17 04466 g015
Figure 16. Variation in the shear stress τ x y along the aluminium bottom face sheet for distinct PUF cores. (a) v f = 0.5 . (b) v f = 0.6 . (c) v f = 0.7 . (d) v f = 0.8 . (e) v f = 0.9 . (f) v f = 1.0 .
Figure 16. Variation in the shear stress τ x y along the aluminium bottom face sheet for distinct PUF cores. (a) v f = 0.5 . (b) v f = 0.6 . (c) v f = 0.7 . (d) v f = 0.8 . (e) v f = 0.9 . (f) v f = 1.0 .
Materials 17 04466 g016
Figure 17. Sandwich cantilever beam with aluminium face sheets and a functionally graded PUF core. Ten material domains were assumed: aluminium top-face sheet (10), PUF-cores with potential distinct densities (2–9), and aluminium bottom face sheet (1).
Figure 17. Sandwich cantilever beam with aluminium face sheets and a functionally graded PUF core. Ten material domains were assumed: aluminium top-face sheet (10), PUF-cores with potential distinct densities (2–9), and aluminium bottom face sheet (1).
Materials 17 04466 g017
Figure 18. Variation in the normal stress σ x x along the beam thickness. (a) FG1 beam for y [ 0.9 , 1.0 ] m. (b) FG2 beam for y [ 0.9 , 1.0 ] m. (c) FG1 beam for y [ 0.1 , 0.9 ] m. (d) FG2 beam for y [ 0.1 , 0.9 ] m. (e) FG1 beam for y [ 0.0 , 0.1 ] m. (f) FG2 beam for y [ 0.0 , 0.1 ] m.
Figure 18. Variation in the normal stress σ x x along the beam thickness. (a) FG1 beam for y [ 0.9 , 1.0 ] m. (b) FG2 beam for y [ 0.9 , 1.0 ] m. (c) FG1 beam for y [ 0.1 , 0.9 ] m. (d) FG2 beam for y [ 0.1 , 0.9 ] m. (e) FG1 beam for y [ 0.0 , 0.1 ] m. (f) FG2 beam for y [ 0.0 , 0.1 ] m.
Materials 17 04466 g018
Figure 19. Variatio n in the shear stress τ x y along the beam thickness. (a) FG1 beam for y [ 0.9 , 1.0 ] m. (b) FG2 beam for y [ 0.9 , 1.0 ] m. (c) FG1 beam for y [ 0.1 , 0.9 ] m. (d) FG2 beam for y [ 0.1 , 0.9 ] m. (e) FG1 beam for y [ 0.0 , 0.1 ] m. (f) FG2 beam for y [ 0.0 , 0.1 ] m.
Figure 19. Variatio n in the shear stress τ x y along the beam thickness. (a) FG1 beam for y [ 0.9 , 1.0 ] m. (b) FG2 beam for y [ 0.9 , 1.0 ] m. (c) FG1 beam for y [ 0.1 , 0.9 ] m. (d) FG2 beam for y [ 0.1 , 0.9 ] m. (e) FG1 beam for y [ 0.0 , 0.1 ] m. (f) FG2 beam for y [ 0.0 , 0.1 ] m.
Materials 17 04466 g019
Table 1. Elastic mechanical properties of the RVE with v f = 0.75 .
Table 1. Elastic mechanical properties of the RVE with v f = 0.75 .
   Nodes E x [MPa] E y [MPa] E z [MPa] G xy [MPa] ν xy ν yx
FEM-3n13888.15588.251129.50824.7450.3130.313
48384.57184.575128.87422.2250.3260.326
178583.03383.044128.70021.2330.3320.332
400582.70282.702128.66921.0190.3330.333
699382.57982.576128.65820.9400.3340.334
FEM-6n50983.92483.922129.50821.6210.3320.332
184982.80882.808128.87421.0330.3340.334
697782.51982.519128.70020.8870.3350.335
CRPIM13863.73962.498115.22912.1470.3680.368
48371.28671.455121.17215.6020.3500.350
178576.76076.759124.84918.1150.3420.342
400578.57278.544126.07218.9880.3400.340
699379.52579.634126.75019.4540.3380.338
MRPIM1613863.78262.556115.22712.1750.3670.367
48371.32971.505121.17115.6240.3490.349
178576.76576.764124.84918.1190.3420.342
400578.57578.547126.07218.9910.3400.340
699379.52779.635126.75019.4560.3380.338
MRPIM913860.44362.784115.22011.7480.3590.359
48371.23870.952121.16315.2360.3530.353
178576.55376.448124.84917.9540.3440.344
400578.39978.462126.06918.8850.3400.340
699379.42579.626126.73519.4220.3370.337
MRPIM413869.31867.072115.09814.0460.3360.336
48372.49172.916121.17816.3230.3380.338
178576.91976.985124.85318.4450.3410.341
400578.62978.786126.05219.0680.3370.337
699379.60479.649126.75519.5020.3380.338
Table 2. Elastic constitutive constants of the RVE with v f = 0.75 . C i j in MPa.
Table 2. Elastic constitutive constants of the RVE with v f = 0.75 . C i j in MPa.
   Nodes C 11 C 12 C 21 C 22 C 33 C 44
FEM-3n138111.63244.50744.507111.762157.62424.745
483107.94044.15544.155107.946156.25122.225
1785106.38444.04644.046106.400155.77921.233
4005106.07244.04944.049106.072155.69121.019
6993105.95844.05444.054105.955155.66020.940
FEM-6n509107.60944.60144.601107.605156.90521.621
1849106.28444.21444.214106.283155.96321.033
6977105.93244.10044.100105.932155.70620.887
CRPIM13882.74235.56934.90980.928136.24312.146
48391.91339.15939.16692.144144.78715.603
178598.74441.55741.57098.746150.10618.115
4005100.99642.33942.343100.959151.86918.988
6993102.14142.74442.797102.301152.85319.454
MRPIM1613882.74535.52634.88780.957136.24012.173
48391.95139.15839.16992.192144.79415.625
178598.74941.55641.57098.750150.10618.119
4005100.99842.33842.343100.962151.87018.991
6993102.14442.74542.797102.302152.85419.456
MRPIM913878.23334.65834.09881.318135.71811.745
48392.08239.33439.43691.713144.80315.237
178598.64641.66441.68098.506150.09517.954
4005100.83442.37242.405100.927151.86018.885
6993101.96642.67442.889102.283152.83819.422
MRPIM413887.67234.97235.76484.823137.06114.043
48392.61638.60938.54693.179144.83816.325
178598.89241.56641.53498.976150.13718.445
4005100.88842.18342.400101.145151.86719.068
6993102.30342.84142.842102.364152.88619.502
Table 3. Elastic homogenised mechanical properties and constitutive constants of the RVE with v f = 0.75 obtained with CRPIM considering a higher-order integration scheme. C i j in MPa.
Table 3. Elastic homogenised mechanical properties and constitutive constants of the RVE with v f = 0.75 obtained with CRPIM considering a higher-order integration scheme. C i j in MPa.
    E x [MPa] E y [MPa] E z [MPa] G x y [MPa] ν x y ν y x
CRPIM*79.54179.636126.75119.4680.3380.338
    C 11 C 12 C 21 C 22 C 33 C 44
CRPIM*102.16642.74642.786102.301152.85519.468
Table 4. Relation between the volume fraction and the circular hole radius.
Table 4. Relation between the volume fraction and the circular hole radius.
v f 1.0000.9500.9000.8500.8000.7500.7000.6500.6000.5500.500
R [mm]0.0001.2621.7842.1852.5232.8213.0903.3383.5683.7853.989
Table 5. Homogenised mechanical properties with respect to volume fraction.
Table 5. Homogenised mechanical properties with respect to volume fraction.
    v f E x [MPa] E y [MPa] E z [MPa] G xy [MPa] ν xy ν yx
FEM-3n1.00171.430171.430171.43065.9350.3000.300
0.95149.211149.212162.91356.4810.3000.300
0.90129.549129.546154.40846.5240.3040.304
0.85111.759111.757145.65136.6800.3110.311
0.8096.67396.677137.28828.2900.3200.320
0.7582.70282.702128.66921.0190.3330.333
0.7070.09470.102120.08415.2350.3510.351
0.6558.51458.513111.43510.7570.3740.374
0.6048.07848.083102.8847.4580.4030.403
0.5538.87038.86594.5805.1100.4370.437
0.5030.10630.10685.8013.3320.4800.480
CRPIM1.00171.430171.430171.43065.9350.3000.300
0.95146.086146.025161.71754.9230.3000.300
0.90125.863125.848152.77444.4580.3040.304
0.85107.790107.747143.65134.3980.3110.311
0.8092.59092.555134.99226.0600.3200.320
0.7578.57278.544126.07218.9880.3330.333
0.7066.00065.923117.25213.5080.3510.351
0.6554.52954.436108.4209.3550.3740.374
0.6044.13744.07199.6486.3520.4030.403
0.5534.95234.96791.1134.2470.4370.437
0.5026.34226.33582.1512.6910.4800.480
MRPIM161.00171.430171.430171.43065.9350.3000.300
0.95146.087146.027161.71754.9260.3010.301
0.90125.864125.850152.77444.4600.3060.306
0.85107.793107.750143.65134.4020.3140.314
0.8092.59292.559134.99226.0630.3250.325
0.7578.57578.547126.07218.9910.3400.340
0.7066.00365.927117.25213.5100.3600.360
0.6554.53254.440108.4209.3570.3870.387
0.6044.14144.07599.6476.3540.4190.419
0.5534.95734.97291.1134.2480.4580.458
0.5026.34626.33982.1512.6920.5070.507
Table 6. Displacement obtained at point A along direction O x for sandwich cantilever beams with PUF cores with distinct volume fractions. u in [mm].
Table 6. Displacement obtained at point A along direction O x for sandwich cantilever beams with PUF cores with distinct volume fractions. u in [mm].
Vol. Frac.FEM-4nCRPIMMRPIM4MRPIM9MRPIM16
0.5−5.194−4.785−5.169−4.895−4.782
0.6−3.642−3.382−3.637−3.441−3.381
0.7−2.831−2.657−2.831−2.694−2.657
0.8−2.369−2.248−2.369−2.273−2.245
0.9−2.078−1.991−2.079−2.009−1.991
1.0−1.883−1.819−1.884−1.833−1.819
Table 7. Displacement obtained at point A along direction O y for sandwich cantilever beams with PUF cores with distinct volume fractions. v in [mm].
Table 7. Displacement obtained at point A along direction O y for sandwich cantilever beams with PUF cores with distinct volume fractions. v in [mm].
Vol. Frac.FEM-4nCRPIMMRPIM4MRPIM9MRPIM16
0.5−284.043−261.657−281.414−270.417−261.661
0.6−185.148−170.287−183.955−175.344−170.282
0.7−130.169−119.857−129.539−123.084−119.847
0.8−96.882−89.406−96.509−91.639−89.398
0.9−74.711−69.153−74.478−70.780−69.144
1.0−58.909−54.722−58.756−55.952−54.715
Table 8. Normal stress σ x x obtained at point B for sandwich cantilever beams with PUF cores with distinct volume fractions. σ x x in MPa.
Table 8. Normal stress σ x x obtained at point B for sandwich cantilever beams with PUF cores with distinct volume fractions. σ x x in MPa.
Vol. Frac.FEM-4nCRPIMMRPIM4MRPIM9MRPIM16
0.5−316.024−277.339−319.636−260.559−276.372
0.6−252.946−224.091−255.732−210.699−223.372
0.7−213.791−191.043−216.064−179.627−190.476
0.8−187.105−168.441−189.052−158.326−167.976
0.9−167.088−151.419−168.806−142.262−151.030
1.0−151.028−137.711−152.571−129.310−137.382
Table 9. Shear stress τ x y obtained at point C for sandwich cantilever beams with PUF cores with distinct volume fractions. τ x y in kPa.
Table 9. Shear stress τ x y obtained at point C for sandwich cantilever beams with PUF cores with distinct volume fractions. τ x y in kPa.
Vol. Frac.FEM-4nCRPIMMRPIM4MRPIM9MRPIM16
0.5−146.937−168.516−143.741−140.503−168.624
0.6−190.784−206.409−186.841−178.593−206.404
0.7−231.156−242.687−226.565−213.384−242.622
0.8−267.833−276.681−262.707−244.885−276.580
0.9−302.068−309.140−296.495−274.212−309.015
1.0−334.716−340.668−328.765−302.107−340.528
Table 10. Mechanical properties of each layer.
Table 10. Mechanical properties of each layer.
BeamLayerVol. Frac.MaterialE [MPa] ν
FG1D11aluminium696000.330
D20.9PUF129.5490.304
D30.8PUF96.6730.320
D40.7PUF70.0940.351
D50.6PUF48.0780.403
D60.6PUF48.0780.403
D70.7PUF70.0940.351
D80.8PUF96.6730.320
D90.9PUF129.5490.304
D101aluminium696000.330
FG2D11aluminium696000.330
D20.8PUF96.6730.320
D30.7PUF70.0940.351
D40.6PUF48.0780.403
D50.5PUF30.1060.480
D60.5PUF30.1060.480
D70.6PUF48.0780.403
D80.7PUF70.0940.351
D90.8PUF96.6730.320
D101aluminium696000.330
Table 11. Displacement obtained at point A along direction O x for sandwich cantilever beams with functionally graded PUF cores. u in [mm].
Table 11. Displacement obtained at point A along direction O x for sandwich cantilever beams with functionally graded PUF cores. u in [mm].
Vol. Frac.FEM-4nCRPIMMRPIM4MRPIM9MRPIM16
FG1−2.703−2.607−2.704−2.624−2.606
FG2−3.495−3.361−3.492−3.386−3.359
Table 12. Displacement obtained at point A along direction O y for sandwich cantilever beams with functionally graded PUF cores. v in [mm].
Table 12. Displacement obtained at point A along direction O y for sandwich cantilever beams with functionally graded PUF cores. v in [mm].
Vol. Frac.FEM-4nCRPIMMRPIM4MRPIM9MRPIM16
FG1−122.141−116.725−121.579−118.573−116.706
FG2−176.362−169.300−175.271−172.018−169.269
Table 13. Normal stress σ x x obtained at point B for sandwich cantilever beams with functionally graded PUF cores. σ x x in MPa.
Table 13. Normal stress σ x x obtained at point B for sandwich cantilever beams with functionally graded PUF cores. σ x x in MPa.
Vol. Frac.FEM-4nCRPIMMRPIM4MRPIM9MRPIM16
FG1−207.505−188.319−209.727−176.261−187.781
FG2−246.623−222.737−249.362−208.174−222.048
Table 14. Shear stress τ x y obtained at point B for sandwich cantilever beams with functionally graded PUF cores. τ x y in kPa.
Table 14. Shear stress τ x y obtained at point B for sandwich cantilever beams with functionally graded PUF cores. τ x y in kPa.
Vol. Frac.FEM-4nCRPIMMRPIM4MRPIM9MRPIM16
FG1−196.320−206.232−191.908−188.425−206.021
FG2−163.736−174.509−159.874−160.641−174.319
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Belinha, J. Multiscale Analysis of Sandwich Beams with Polyurethane Foam Core: A Comparative Study of Finite Element Methods and Radial Point Interpolation Method. Materials 2024, 17, 4466. https://doi.org/10.3390/ma17184466

AMA Style

Belinha J. Multiscale Analysis of Sandwich Beams with Polyurethane Foam Core: A Comparative Study of Finite Element Methods and Radial Point Interpolation Method. Materials. 2024; 17(18):4466. https://doi.org/10.3390/ma17184466

Chicago/Turabian Style

Belinha, Jorge. 2024. "Multiscale Analysis of Sandwich Beams with Polyurethane Foam Core: A Comparative Study of Finite Element Methods and Radial Point Interpolation Method" Materials 17, no. 18: 4466. https://doi.org/10.3390/ma17184466

APA Style

Belinha, J. (2024). Multiscale Analysis of Sandwich Beams with Polyurethane Foam Core: A Comparative Study of Finite Element Methods and Radial Point Interpolation Method. Materials, 17(18), 4466. https://doi.org/10.3390/ma17184466

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop