Abstract: In this study, a truly mesh-free method is implemented for the elastic analysis of discontinues media. The media is divided into separate parts which are connected to each other by the interface layers. The displacement field in each part is constructed by the radial basis point interpolation method, enriched by the polynomial terms. The final system of equations is derived by the substitution of the displacement field into the weak form of the governing equation. The dependency of the method on the background mesh has been taken care of by stabilized nodal integration. Based on derived equations a computer code is developed and some cases are investigated. The results of analyses show that the nodal integration improves the accuracy of results in comparison with the Gauss method. The stability of proposed method is also guaranteed by the suitable choice for shape parameters and the radius of supports.
INTRODUCTION
Contact problems are very common and hence important in mechanical and civil engineering. Anytime there is an interface between two materials, there is a discontinuity and the usual analytical methods used in continues solid mediums fail. Strong discontinuities can be solved easier by numerical methods than analytical ones. Conventional engineering numerical methods such as Finite Element (FEM) and Finite Difference (FDM) are based on element or mesh and hence may give rise to problems related to mesh distortion. To address these issues many mesh-free methods or element-free methods have been developed in the past few years. To name a few: the Smoothed Particle Hydrodynamics (SPH) method (Libersky and Petschek, 1991), Element Free Galerkin (EFG) method (Belytschko et al., 1994), Reproducing Kernel Particle Method (RKPM) (Liu et al., 1995), Local Petrov-Galerkin Method (LPGM) (Atluri and Zhu, 1998) and Point Interpolation Method (PIM) (Liu and Gu, 2001). Detailed descriptions and discussions on these mesh-free methods can be found in, for example (Liu, 2003).
There has been a great deal of research into the implementation of the Radial Basis Point Interpolation Method (RBPIM) in many engineering problem solutions. The RBPIM has the following advantages:
| • | The shape function has the Kronecker delta property |
| • | The moment matrix used in constructing shape functions is always invertible |
| • | The linear consistency is guaranteed by using the RBPIM basis functions enriched with the linear polynomials |
The generation of discrete system of equations in the RBPIM, like the majority of mesh-free methods, is based on Galerkin procedure. Gaussian integration is commonly used in Galerkin mesh-free methods for the integration of weak forms and thus a background mesh is needed to perform integration. However the misalignment of the local supports and the integration cells may deteriorate the accuracy of results. Dolbow and Belytschko (1999) investigated this problem and suggested the bounding box technique for the construction of the integration cells which match the shape functions supports. Due to the complexity involved in Gaussian integration, attempts have been devoted to explore stabilized nodal integration methods for Galerkin-based mesh-free methods.
The source of instability in nodal integration is the vanishing of the first derivatives of the shape functions at the nodes. Since the strain energy term in the weak form consists only the first derivatives of shape functions, nodal sampling does not adequately reflect the energy of the short-wavelength mode and its contribution to the stiffness matrix is severely underestimated. This instability is manifested by spurious short wavelength modes. Randles et al. (1999) used the idea of computing derivative away from the centroid of the support domain and proposed a stress point method to improve accuracy and to reduce spurious oscillations in SPH. Bonet and Kulasegaram (1999) introduced an integration correction to stabilize the nodal integration method and to eliminate the spurious modes in SPH. Beissel and Belytschko (1996) stabilized the nodal integration of Galerkin mesh-free methods by adding the square of the residual of the equilibrium equation to the potential energy. The stabilization term includes second derivatives of shape functions and thus the least squares term serves to stabilized nodal integration. Although, the addition of stabilization term improves the accuracy of solution for the problems with spurious near-singular modes, for problems that do not contain unstable modes in their original solution, this addition actually deteriorate the accuracy. A stabilized conforming nodal integration for Galerkin mesh-free methods was proposed by Chen et al. (2001). They identified that for linear exactness in the Galerkin approximation, the shape functions have to be linearly consistent and the nodal integration must satisfy integration constraints. They proposed a strain smoothing stabilization technique to stabilize the nodal integration. The basic idea is to replace the strain at a specific point with the average strain in a Voronoi cell that contains the point. Divergence theorem is then used to replace the area, or volume integration around the point by a contour integration of the Voronoi cell boundary.
In this study, the conventional RBPIM, enriched with polynomial basis functions is implemented for the elastic analysis of discontinues medium. The interface layers between different materials are defined by the concept of linkage element. The stabilized nodal integration technique has been used to integrate the weak form of equations.
RBPIM ENRICHED WITH POLYNOMIAL BASIS FUNCTIONS
A field function u(x) can be approximated using both radial and polynomial basis as:
(1) |
where, n is the number of field nodes in the local support domain for point x; the vector R(x) is defined as:
(2) |
Where:
(3) |
is the radial basis functions and
(4) |
is the distance between the point x and field node Xk; vector P(x), where:
(5) |
is the vector of polynomial basis functions in 2D space xT = [x, y] and m is the number of terms of polynomial basis functions. Vectors a and b, where:
(6) |
(7) |
are, respectively, coefficients for R(x) and P(x). The radial basis functions are used to guarantee the invertability of moment matrix and the polynomial basis functions are used to ensure the linear consistency of generated shape functions. The coefficient vectors a and b are determined by enforcing Eq. 1 to be satisfied at all the n nodes within the local support domain. By the lengthy but straightforward procedure given by Liu (2003) we have:
(8) |
Where:
(9) |
is a vector of nodal displacements of the nodes in local support domain and
(10) |
contains RBPIM shape functions for the n local nodes in which Φk(x) is as follows:
(11) |
where Saik is the (i, k) entry of matrix Sa where:
(12) |
and Sbjk is the (j, k) entry of matrix Sb where:
(13) |
The moment matrices RM and PM are consisted of row vectors RT(xi) and PT(xi) (i = 1, 2,…,n), respectively.
There are many types of radial basis functions. In this study the multi-quadratic form is used as:
(14) |
where, c and q are shape parameters, which are needed to be determined by numerical tests.
VARIATIONAL PRINCIPALS FOR DISCONTINUOUS PROBLEMS
Generally, the total potential energy functional for a discontinuous media is expressed as:
(15) |
where, Πe and Πi are the elastic strain energy of a solid media and interface layer, respectively and Πf is the potential energy related to external forces. These functionals can be defined as:
(16) |
(17) |
(18) |
where, ε and σ are, respectively, the strain and stress tensors
of solid domains and εi and σi are the
strain and stress tensors of interface layers between solid domains, respectively.
b is the body force vector, U is the vector of degrees of
freedom, Ω is the problem domain, β is the length parameter
along the interface layer,
The variational (weak) form of Eq. 15 can be written as:
(19) |
where, D and Di are the elasticity matrixes of medium and interface layer, respectively.
| Fig. 1: | Interface modeling |
Using the stress-strain relation, the discrete form of equations can be written as:
(20) |
Where:
(21) |
(22) |
where, Φ and B are, respectively, the shape functions and the gradient of shape functions matrices. Using the generated shape functions, Φ and B can be constructed as in conventional FEM. Matrices Bi and Di are related to the interface layer. The rest of this section is devoted to the expression of the interface modeling using mesh-free method.
The interface layer is considered as two parallel planes with an insignificant distance from each other. A series of springs are assumed between the planes. The springs are resistant to the relative normal and tangential displacement of planes. The stiffness of springs into the unit area (stiffness coefficient) along normal and tangential directions are kn and ks, respectively. As shown in Fig. 1, point P is considered between the planes of the interface layer to determine the strain-displacement relationship. The width of interface (h) is much less than its length (L), hence the strain in the n direction in the interface can be assumed to be constant. The relative displacement of point P can be related to the displacements of point A and B on the top and bottom planes as follows:
(23) |
where, δn and δs are the normal and shear
relative displacement of point P, respectively and
(24) |
On the other hand, the displacements of points A and B can be related to the displacements of the nodes in the support domain of those points (Fig. 1). So, we have:
(25) |
where, ΦA and ΦB are the
shape function matrices of the nodes in the support domains of point A
and B respectively and
Considering Eq. 23-25 we have:
(26) |
Where:
(27) |
(28) |
Defining the strain tensor in the interface layer as:
(29) |
where, γns and εn are the shear and normal strain, respectively, we can have:
(30) |
and by combination of Eq. 26 and 30 the strain-displacement relation in the interface layer can be written as follows:
(31) |
Considering the strain tensor definition in the interface layer, the stress-strain relation in elastic analysis can be written as:
(32) |
Where:
(33) |
(34) |
where, τ and σn are shear and normal stress in the interface layer, respectively.
STABILIZED NODAL INTEGRATION
As noticed before, the source of instability in nodal integration is the vanishing of the first derivatives of the shape functions at the nodes. Therefore attempts have been made either to calculate the derivatives away from the nodes or to derive the formulations in which the derivatives are not manifested. Chen et al. (2001) used the latter technique and proposed a strain smoothing as follows:
(35) |
where,
(36) |
and Ψ is a distribution function which is chosen as:
(37) |
where, AL is the area of the representative domain of node L obtained from the Voronoi diagram. Combination of Eq. 8 and 35-37 and applying the integration by parts gives:
(38) |
where, GL is a group of nodes in which their associated shape function supports cover node L. In two dimensional space we have:
(39) |
(40) |
(41) |
(42) |
(43) |
As shown in Fig. 2, ΓL is the boundary of Voronoi cell which contains node L. nx and ny are, respectively, the x and y components of vector n (outward vector which is normal to boundary of Voronoi cell). Note that in Eq. 41 no derivatives of shape functions are involved in evaluating the smoothed gradient matrix at the nodal points. Considering the Eq. 41 and the nodal integration method, Eq. 21 and 22 can be rewritten as:
(44) |
(45) |
where, NP is the number points in the local support domain of Node L, NPb is the number of points on the natural boundary, SL are the weights associated with the boundary points that can also be obtained from the Voronoi diagram and Ki is the stiffness matrix of the interface layer (second term in Eq. 21).
To calculate Eq. 42 and 43, any numerical integration method can be used, for instance by applying a two-point trapezoidal rule for each segment in Fig. 2, Eq. 42 and 43 can be written as:
(46) |
(47) |
where, Ns is the total number of segments of Voronoi cell
contains node L,
| Fig. 2: | The Voronoi cell contains node L (Chen et al., 2001) |
NUMERICAL STUDY
Example 1: Bimaterial bar: Consider a bimaterial bar (Ω =
Ω1UΩ2) as shown in Fig.
3, with length L = 2 m, height h = 8 m with material moduli E1
and Poisson ratio υ1 in Ω1 and E2,
υ2 in Ω2. The material constants were
chosen so that there is no singularity at the edge of the interface; this
condition
(48) |
As shown in Fig. 4 three mesh-free models with 30, 90 and 306 nodes in which the distances between nodes are 1, 0.5 and 0.25 meter respectively, are considered. The shear and normal stiffness coefficients of interface are determined by the following equations:
(49) |
| Fig. 3: | Bimaterial bar |
| Fig. 4: | Mesh-free models for bimaterial bar |
(50) |
where, ks and kn are, respectively, the shear and normal stiffness coefficients and t is the thickness of the interface layer which is chosen 0.005 m. For integration, both, the Gauss and the nodal integration method are used. In Gauss method a 4x16 background mesh with 16 Gauss points in each block and in nodal method the Voronoi diagrams shown in Fig. 4 are employed for integration. Each integration point is employed as the center of circular support domain. To guarantee sufficient and suitable nodes are covered in support domain, the radius of the support domain is specially devised with a slightly adjustable value in the program which is automatically self-tuned such that at least 15 nodes are selected for each local support domain.
In order to compare the results of mesh-free analysis with exact solution, the error indicator in energy is defined as follows:
(51) |
where, D is the elastic matrix, εexct and εnum are, respectively, the exact and the mesh-free (numerical) strain tensors. The variation of en with nodes distances is shown in Fig. 5 for both methods of integration. As it is obvious, the results of nodal integration are much better than the results of Gauss integration and, they are also in very good agreement with the results of exact solution (i.e., maximum value of en for nodal integration is 0.00228).
| Fig. 5: | Variation of en with the nodes distances |
| Fig. 6: | A column of reinforced elastic material; type 1: material;
type 2: two narrow reinforcement strips with the thickness of 0.1
m; Interface: layers with thickness of 0.05 m between materials |
The convergence of proposed mesh-free method is also confirmed by decreasing the nodal distance and consequently increasing the accuracy of results.
Example 2: Reinforced column: In this example, an elastic reinforced column is investigated. The dimensions and boundary conditions of the problem are shown in Fig. 6. The elastic modulus of material type1 is 3 GPa and its Poisson`s ratio is 0.3. The elastic modulus and Poisson`s ratio of material type 2 are 60 and 0.25 GPa, respectively. Interface layers are considered between different materials. The normal and shear stiffness of the interfaces are 104 and 106 kN m-3, respectively.
| Fig. 7: | Variation of η with respect to the number of elements |
Since, there is no exact solution for this problem, the FEM is used to compare the results with the mesh-free method. It is assumed that by increasing the number of elements the results of the FEM code gets closer to the exact solution. By increasing the number of elements, a parameter is calculated as:
(52) |
where, ||.|| stands for Euclidean norm
In this example, the computer code SIGMAW has been used as the FEM code. The reason for the selection of this software is that it has slip elements and the interface layer can be modeled easily, also the finite element modeling is the most similar to mesh-free simulation. By increasing the number of elements and a comparison of the results with the initial solution, η is calculated at every stage of the analysis. Figure 7 shows the variation of η with respect to the number of elements. According to the Fig. 7, the finite element model with 3420 elements gives the best results.
To investigate the convergence of proposed method, three mesh-free models and their Voronoi diagrams are assumed (Fig. 8). The first model contains 39 nodes which are 1 m apart from each other.
| Fig. 8: | Mesh-free models for reinforced column |
| Fig. 9: | Variation of ed with respect to the number of nodes |
In
the second and third models the distances between nodes decreased to 0.5
and 0.25 m, respectively and consequently 95 and 225 nodes are generated. The relative error of displacement between the mesh-free and the FEM is
determined by:
(53) |
where, ed is the relative error of displacement, UMFM is the displacement vector in mesh-free method, UFEM is the displacement vector in FEM. Figure 9 shows the variation of ed with respect to the distances between nodes. It is obvious that by decreasing the distances between nodes, the relative error is also decreased. This condition confirms the convergence of the method.
To investigate the effect of the size of support domain on the accuracy of results, the values of ed are plotted against different values of the radiuses of support in Fig. 10.
| Fig. 10: | Variation of ed with respect to the radiuses of supports |
| Fig. 11: | Effect of shape parameters on the value of ed |
Due to the insufficient number of nodes in the support domain, the 39 nodes model has no results for the radiuses value lower than 1.5. On the other hand, as it is obvious from Fig. 10, an increase in the number of nodes in the support domain (i.e., increasing the value of the radius of support) does not necessarily increase the accuracy of results.
The effect of shape parameters on the accuracy of results is investigated by computing ed for various values of c and q in 225 nodes model at the constant value of radius of support (0.7). The results of analysis are shown in Fig. 11.
| Fig. 12: | Irregular nodal arrangement for mesh-free model of reinforced column |
| Fig. 13: | Effect of nodal arrangement on the value of ed |
According to this figure, the analysis is more sensitive to the value of c than to the value of q and for the values of c lower than 0.5 or upper than 2, the error increases significantly.
In the last part of this section, the effect of nodal arrangement on the accuracy of results is investigated. Hence an irregular arrangement of nodes is considered for 225 nodes model. This model and the Voronoi cell of each node are shown in Fig. 12. The parameter ed is calculated for different radiuses of support and the results are shown in Fig. 13 for regular and irregular nodal arrangements. As is obvious from the Fig. 13, the disturbance in nodal arrangement deteriorates the accuracy of results, but the decrease in accuracy is not to the extent to destroy the stability of the proposed method.
CONCLUSION
General conclusions of this study can be expressed as follows:
| • | A truly mesh-free method is implemented for the elastic analysis of discontinues medium. In this method, the shape functions are constructed by the radial basis point interpolation method which is enriched by polynomial basis terms. The integration is also performed by the nodal method to eliminate the dependency of method to the background mesh. |
| • | The concept of a linkage element is used for interface modeling and the calculations are just based on the displacement of adjacent nodes to the interface layer. So, there is no need of enrichments such as Lagrangian coefficients or penalty method which increase the complexity of solution. |
| • | Comparison of results between two methods of integration showed that the nodal integration impose less errors to the norm of energy (which is an indicator for the numerical error in stress and strain) and hence the accuracy of proposed method is confirmed. |
| • | Suitable choice for the values of shape parameters and the radius of support guarantees the stability of method. This is also true for the irregular arrangement of nodes. |