## 3 Meshless Simulations  
Several computational difficulties arose in the analysis of missile dynamic effects on the foam pad. The non-linear behaviour and the extremely large strain on the polymeric foam, together with the presence of an impulsive load, required a very small time step to be fixed to ensure the numerical solution of the pad equilibrium equations. For this reason, the central differencing scheme [18] was adopted to avoid the inversion of the pad stiffness matrix in each calculation step.  
Moreover, the large deformations of the foam made the Lagrangian mesh extremely skewed and compressed. The remarkable distortion of the constitutive elements led to a reduction in the accuracy of the solution algorithm. To avoid accuracy reduction, periodic remeshing of the pad volume during the evolution of the simulation was needed. However, this solution produced a very high computational cost and the remeshing advantages were partially reduced from the degradation of accuracy due to the projection operation of the field variables from the old meshes to the new ones.  
To overcome these difficulties, the Element-Free Galerkin (EFG) method of [PERSON] et al. [16] was adopted to analyse the pad behaviour under missile impact. Using this method, the solution is constructed entirely in terms of a set of nodes and no elements are needed. The unknown displacement field \(u^{(h)}\) is searched in the form given by  
\[u^{(h)}=\mathbf{p}^{T}(x)*\mathbf{a}(x), \tag{1}\]  
where \(\mathbf{p}(x)\) is the vector of the monomial basis functions and \(\mathbf{a}(x)\) is a vector of unknown parameters. These are determined at any point \(x\) by solving the following optimization problem:  
\[minimize\ J=\sum_{l=1}w(x-x_{l})\left[\mathbf{p}^{T}(x)*\mathbf{a}(x)-u_{l}\right], \tag{2}\]  
where \(u_{l}\) are the unknown values assumed by the displacement field at node I having coordinate \(x_{l}\), \(w(x-x_{l})\) is a weighting function with compact support, and NP is the number of nodes within the support of \(w\). When the solution to Equation (2) is substituted into Equation (1), we obtain for the approximation \(u^{(h)}\) the expression  
\[u^{(h)}_{i}=\sum_{l=1}\varphi_{l}(x)u_{l},\]  
with \(\varphi_{l}(x)\) being the shape functions of the EFG method given by  
\[\varphi=\left[\varphi_{1}(x)\cdots\varphi_{n}(x)\right]=\ p^{T}(x)A^{-1}(x)B( x),\]  
where \(A(x)=P^{T}WP\) and \(B(x)=P^{T}W\), with  
\[P=\left[\begin{array}{c}p(x_{1})\\ \cdots\\ \cdots\\ p^{T}(x_{n})\end{array}\right]\text{and }W=\left[\begin{array}{cccc}w(x-x_{1})&0& \ldots&0\\ 0&w(x-x_{2})&\ldots&0\\ \ldots&\ldots&\ldots&\ldots\\ 0&0&\ldots&w(x-x_{n})\end{array}\right]\]  
The discretized equilibrium equations were derived by substituting the approximating function in the weak form of equilibrium conditions and carrying out numerical integration [15].  
It is worth noting that the convergence of the Galerkin method chosen to discretize the problem at hand depended both on the approximation adopted for the unknown functions and the numerical integration of the weak form. To achieve the linear exactness solution in the Galerkin approximation, some integration constraints had to be met [19].  
Therefore, the Lagrangian strain smoothing strategy described in [20] was used for the domain integration by Gauss quadrature. Finally, since the EFG shape functions are not interpolation functions, the essential boundary conditions were enforced by the transformation method proposed in [21].