- Email: [email protected]
A general 3D contact smoothing method based on radial point interpolation
Journal of Computational and Applied Mathematics 257 (2014) 1–13
Contents lists available at ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
A general 3D contact smoothing method based on radial point interpolation Xiaoxiang Qian, Huina Yuan, Mozhen Zhou, Bingyin Zhang ∗ State Key Laboratory of Hydroscience and Engineering, Department of Hydraulic Engineering, Tsinghua University, Beijing 100084, China
highlights • • • • •
A general 3D contact smoothing method using radial point interpolation is proposed. The method constructs local contact surfaces with polynomial and radial bases. The constructed surface is smooth and passes exactly through the surface nodes. The method shows good approximation in both coordinates and normal vector. It can effectively reduce contact force oscillation and solve the chatter problem.
article
info
Article history: Received 9 May 2013 Received in revised form 4 August 2013 Keywords: Contact Finite element method Surface smoothing Radial point interpolation method
abstract This paper presents a general 3D contact smoothing method based on the meshfree radial point interpolation method to improve the numerical simulation of contact problems. In particular, a locally smooth contact surface is constructed from the scattered surface nodes by point interpolation using the combination of polynomial and radial bases. With such bases, this method reproduces smooth surfaces even for coarse meshes and the constructed surface passes exactly through the surface nodes. Results for contact problems involving deformable bodies are included to demonstrate its advantages. © 2013 Elsevier B.V. All rights reserved.
1. Introduction Significant progress has been made on the finite element (FE) analysis of contact problems. However, an accurate and efficient representation of the nonlinear geometry of contacting surfaces remains a challenge. When a standard master–slave approach is used, piecewise continuous FE discretization of the master contact surface may cause two types of problems. One type of problems is associated with the searching of contact segments. With the commonly used node-to-segment (NTS) algorithm, a slave node may have multiple or no candidates of contact segments due to C 0 continuity of the master surface [1,2]. To circumvent these problems, alternative contact search algorithms have been developed, such as the hierarchyterritory algorithm [3], the pinball algorithm [4,5], and the inside–outside algorithm [6]. Recently, some new techniques have also been proposed to deal with these special situations within the framework of the NTS algorithm [7]. The other type of problems, often met in large motion/deformation problems, is associated with the jump of normal vector when a slave node slides from one segment to an adjacent segment on the master contact surface. These jumps, induced by spatial discretization, are non-physical and can cause serious errors in the calculation of contact forces. Moreover, slave nodes situated near an edge may ‘chatter’ during equilibrium iterations and thus impede convergence. To overcome these problems, a general idea is to smooth the master surface. For 2D contact, spline-based smoothing techniques have typically been used, such as the cubic Hermite interpolation [8], the cubic Bézier spline [9], the Overhauser spline [10,11]
∗
Corresponding author. Tel.: +86 10 62787349; fax: +86 10 62785593. E-mail address: [email protected] (B. Zhang).
0377-0427/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cam.2013.08.014
2
X. Qian et al. / Journal of Computational and Applied Mathematics 257 (2014) 1–13
Fig. 1. Different types of surface meshes: (a) regular quadrilateral mesh, (b) triangular mesh, (c) irregular quadrilateral mesh, and (d) hybrid mesh.
Fig. 2. Notation for the frictional contact problem.
and the non-uniform rational B-splines [12]. Unfortunately, an extension of such methods to 3D arbitrary meshes proves to be much more difficult. Fig. 1 shows different types of surface meshes for 3D contact problems. 3D generalizations of 2D spline curves through tensor product only work for regular quadrilateral meshes (Fig. 1(a)) where each inside quadrilateral facet is surrounded by eight neighboring quadrilateral facets [13,14]. Smooth descriptions using quartic Bézier surfaces [15] were proposed for contact surfaces composed of triangular facets particularly (Fig. 1(b)). Another 3D approach [16] uses Gregory patches to interpolate irregular quadrilateral meshes (Fig. 1(c)) where a quadrilateral facet has an arbitrary number of neighboring quadrilateral facets. The Gregory patch augments the third-order Bézier patch with additional internal control points to enable tangent plane continuity in general. With special treatments [17] developed for patches involving inflection points, the Nagata patch interpolation algorithm can be applied to both triangular and quadrilateral faceted finite elements [18,17,19]. However, a consensual good solution for contact smoothing of complex surfaces composed of hybrid meshes (Fig. 1(d)), where each facet, surrounded by an arbitrary number of neighboring facets, can be either quadrilateral or triangular, remains a challenge. Due to the geometrical complexity of 3D bodies, hybrid surface meshes are often encountered when mixed 3D element types, e.g. tetrahedral and hexahedral elements, are used. Chamoret et al. [20] proposed a smoothing procedure based on a meshfree technique denoted as diffuse approximation. This approach creates a locally smooth contact surface from the scattered surface nodes using a least-squares approximation with a polynomial basis. Belytschko et al. [21] proposed a similar approach, where the smoothing is done implicitly by constructing smooth signed distance functions for the bodies. Without using element connectivity, these methods are applicable to arbitrary surface meshes, such as those shown in Fig. 1. However, the constructed contact surface does not pass through the surface nodes exactly, which introduces some inaccuracies with regard to the real geometry of contact bodies. The goal of this work is a general 3D contact smoothing method based on a meshfree technique denoted as the radial point interpolation method (RPIM) [22]. The basic concept is to construct a locally smooth contact surface from the scattered surface nodes using point interpolation with both polynomial and radial basis functions. Similar to the methods proposed in [20,21], this approach also applies to arbitrary surface meshes. Moreover, the smoothed contact surface passes exactly through all the surface nodes of the master body in the vicinity of any slave node where contact is considered. 2. Notation and preliminaries This section briefly reviews the fictional contact problem and the notion to be used in this paper [23,24]. Consider a frictional contact problem depicted in Fig. 2, where the number of potential contact bodies is limited to two for convenience. The initial configurations of the two bodies (at time t = 0) are represented by open sets 0 Ω 1 , 0 Ω 2 ⊂ R3 , where X 1 ∈ 0 Ω 1 , X 2 ∈ 0 Ω 2 are two material points. The motions of the two bodies are defined through mappings: x = t ϕ k (X k )
t k
k = 1, 2
(1)
where x , x are the current positions of X , X (at time t ≥ 0) and Ω , Ω are the current configurations of the two bodies. Although these two bodies are interchangeable with respect to their mechanics, one (t Ω 2 ) is designated as the master and the other (t Ω 1 ) as the slave to describe their contact. Correspondingly, the potential contacting surfaces of the two bodies, t γ 2 and t γ 1 , are designated as the master and slave contact surfaces respectively. t 1 t 2
1
2
t
1 t
2
X. Qian et al. / Journal of Computational and Applied Mathematics 257 (2014) 1–13
3
The bodies are governed by the conservation of mass, momentum and energy, the constitutive relations, and the initial and boundary conditions as usual. In addition, contact between the bodies is subjected to the normal and frictional contact constraints. First, the bodies cannot interpenetrate, called the impenetrability condition. This condition can be written as t
gN (t x1S ) = t n2P · (t x1S − t x2P )
= min t x1S − t x2 sign(t n2 · (t x1S − t x2 )) ≥ 0 ∀t x1S ∈ t γ 1
(2)
t x2 ∈t γ 2
where t x2P is the projection (closest) point of a slave point t x1S ∈ t γ 1 onto the master surface t γ 2 , t n2 is the unit outward normal to the master t Ω 2 , and t gN is the gap function. In the example shown in Fig. 2, the points of the contacting surfaces have interpenetrated (a situation often encountered during numerical simulation) and thus t gN < 0. Second, it is assumed that the normal traction across the contact interface cannot be tensile, i.e. t
FN1 = −t FN2 = t F 1 · t n2 = −t F 2 · t n2 ≥ 0
(3)
where t FN1 , t FN2 are the normal tractions on t γ 1 , t γ 2 respectively. Finally, the frictional behavior of the contact interface follows a specified friction model. Here the Coulomb friction model is considered:
¯ T = t υ1T − t υ2T = 0 ⇔ t FT1 < µ t FN1 Stick: t υ
(4)
υ¯ T = υ − υ ̸= 0 ⇔ t FT1 = µ t FN1 Slip: t t υ¯ T · t FT1 < 0, υ¯ T × t FT1 = 0
(5)
t
t
1 T
t
2 T
¯ T is the relative tangential velocity between a pair of contacting points, t FT1 = t F 1 − t FN1 t n2 is the tangential where t υ traction, and µ is the coefficient of friction. The last inequation and equation in (5) represent that the tangential velocity and tangential contact force are collinear. Within the framework of an incremental solution method, the weak form of the contact problem at time t + ∆t (assuming that the solution at time t has been obtained) can be written as 1 ,2
t +∆t
k δ Wint − t +∆t δ WLk − t +∆t δ WIk − t +∆t δ WCk = 0
(6)
k=
where δ Wint , δ WL , δ WI , and δ WC are the internal energy, the work of external forces, the work of inertial forces, and the work of contact forces respectively. In particular, the work of contact forces is given by t +∆ t
δ WC =
1,2
t +∆t
δ WCk =
k=
1 ,2 k=
t +∆t t +∆t γ k
Fik δ uki t +∆t ds =
t +∆t t +∆t γ 1
Fi1 δ ui t +∆t ds
(7)
where δ u is the variation of the relative displacement between a pair of contacting points. The impenetrability condition can also be written in an incremental form as t +∆ t
gN = t +∆t n2 · (t +∆t x1 − t +∆t x2 ) = t +∆t n2 · (u1 − u2 ) + t +∆t n2 · (t x1 − t x2 )
= u1N − u2N + tr gN ≥ 0
(8)
where u1N , u2N are the normal displacement increments of a pair of contacting points over ∆t time and tr gN is the residual gap from the last time step. Accordingly, the tangential velocity υ T in (4) and (5) can be replaced by the tangential displacement increment uT . The details of solving (6) by way of the finite element method can be found in [25–27]. The main concern of this paper is contact surface smoothing that plays a significant role in accurate location of the projection point and continuity of the normal vector. 3. Contact surface smoothing The contact surface smoothing procedure based on the meshfree radial point interpolation method (RPIM) is presented in this section. Resulting from FE discretization of the contact bodies, the contacting surfaces are originally represented by collections of quadrilateral and/or triangular facets. When a master–slave approach is used, only the master contact surface requires smoothing. The interpenetration of any slave node can be computed in terms of the gap function. 3.1. Surface smoothing based on RPIM Consider a slave node S (xS ∈ γ 1 ) for contact searching. For simplicity, the superscripts representing time and body are omitted hereinafter when there is no ambiguity. The basic idea is to create a smooth description of the contact surface using surface nodes of the master body Ω 2 in the vicinity of xS . Thus, it is first necessary to find the surface nodes on γ 2 in the neighborhood of xS . For two newly activated contact bodies, this can be done through a global search and a preliminary
4
X. Qian et al. / Journal of Computational and Applied Mathematics 257 (2014) 1–13
Fig. 3. A domain of influence for surface smoothing.
local search. Here the position code algorithm [28] and the inside–outside algorithm [6] are used to locate the temporary contact facet MF and the temporary projection point Q (xQ ). The vicinity of xS , also referred to as the domain of influence in the context of RPIM, can then be determined via a multilevel neighboring search for MF. A domain of influence for surface smoothing is shown in Fig. 3 as an example, where the temporary contact facet MF is shaded. Next, the meshfree radial point interpolation method [22] is applied to construct a locally smooth contact surface γS over the domain of influence ΩS of xS , using the information (coordinates) of all the surface nodes in ΩS (Fig. 3). Let (xi1 , xi2 , xi3 ) denote the coordinates of a master node xi ∈ ΩS in a reference frame Rd . The local approximation of the contact surface for the slave node xS is described by the following equation: x3 = fS (x1 , x2 ) =
n
Ri (x)ai +
i=1
m
Pj (x)bj = R T (x)a + P T (x)b
(9)
j=1
where n is the number of master nodes in ΩS , Ri , i = 1, . . . , n is a term of the radial basis, Pj , j = 1, . . . , m is a term of the polynomial basis, and ai , i = 1, . . . , n and bj , j = 1, . . . , m are corresponding coefficients. The radial basis function Ri takes the form: Ri (x) = Ri (x − xi ) = [ri2 + (Rdc )2 ]q
(10)
(x1 − xi1 )2 + (x2 − xi2 )2 is the distance between x and xi projected onto the x1 x2 plane, dc is the average where ri = projected node distance in ΩS , and R, q are selected shape parameters. As for the polynomial basis, commonly used linear and quadratic bases are given by P T (x) = [1
x2 ]T
x1
P (x) = [1 T
x1
x2
m=3 x1 x1
(11)
x1 x2
T
x2 x2 ]
m = 6.
(12)
The coefficient vectors a and b can be determined by making the surface pass through all the master nodes in ΩS , which leads to the following system of equations:
χ 0
=
C DT
where χ = [x13
x23
R1 (x1 ) R1 (x2 )
C =
.. . R1 (xn )
D 0
a a =G b b
···
(13)
xn3 ]T is a vector of the x3 components of all the master nodes xi ∈ ΩS , and
R2 (x1 ) R2 (x2 )
.. . R2 (xn )
··· ··· .. . ···
Rn (x1 ) Rn (x2 )
.. , . Rn (xn )
1 1
D = . .
. 1
x11 x21
x12 x22
xn1
xn2
.. .
.. .
··· ··· .. . ···
Pm (x1 ) Pm (x2 )
.. . Pm (xn )
(14)
are matrices formed by the radial and polynomial basis functions. With a and b solved from (13), the x3 component of any point x on the constructed local contact surface can be computed from x3 = fS (x1 , x2 ) = R T (x)
a
P T (x)
= R T (x) b
−1 P T (x) G
χ 0
= ΦT (x)χ
(15)
where Φ(x) is a vector of resulting RPIM shape functions. From the forms of the radial and polynomial basis functions, the constructed surface γS is C n continuous both inside and across patches whereby n represents an arbitrary level of continuity. In particular, the RPIM shape functions have the Kronecker delta property, i.e.
Φi (xj ) =
1, 0,
i=j i ̸= j
which ensures that the constructed surface passes exactly through all the master nodes in the domain of influence.
(16)
X. Qian et al. / Journal of Computational and Applied Mathematics 257 (2014) 1–13
5
Fig. 4. The 2D natural coordinate system for point interpolation.
3.2. Preconditioning Further examination of (13) shows that the interpolation will fail if projection of the local contact surface in the vicinity of xS onto the x1 x2 plane overlaps. To prevent this problem, a preconditioning method is proposed. Let point Q in the surface Q Q element MF be the control point of the domain of influence, as shown in Fig. 3. Its natural coordinates ξ Q (ξ1 , ξ2 ) in the element and global coordinates xQ (x1 , x2 , x3 ) in the reference frame Rd are related by Q
ne
xQ =
Q
Q
Nα (ξ Q )xα
(17)
α=1
where ne is the number of nodes in the surface element (either quadrilateral or triangular), xα , α = 1, . . . , ne is the α -th node, and Nα is the corresponding shape function, depending on the type of the element. A suitable 2D natural coordinate system for point interpolation can be obtained by rotating the plane parallel to the x1 x2 plane and passing through point Q by a certain angle. The procedure is described as follows. (1) Calculate the unit outward normal n at point Q : n=
τ1 × τ2 = ni ei ∥τ 1 × τ 2 ∥
(18)
where ei , i = 1, 2, 3, is the coordinate basis of Rd , and τ j = ∂ xQ /∂ξj , j = 1, 2, is the directional derivative of xQ with Q
respect to ξ (2) Determine the pivot r and rotation angle θ : Q j .
r =
e3 × n
∥e3 × n∥
,
θ = arccos(n · e3 ).
(19)
(3) Rotate the plane F0 (parallel to the x1 x2 plane and passing through point Q ) around r by angle θ . The resulting reference plane FR is shown in Fig. 4, with the coordinate basis given by
e′1
e′2 =
1
n21 + n22
n21 n3 + n22
n1 n2 (n3 − 1) −n1 (n21 + n22 )
n1 n2 (n3 − 1) n22 n3 + n21
−n2 (n21 + n22 )
n1 (n21 + n22 )
n2 (n21 + n22 ) [e1 , e2 ] .
(20)
n3 (n21 + n22 )
(4) Project the master nodes (interpolation points) onto the FR plane and calculate the natural coordinates x′i (x′1i , x′2i ): ′ x′i = e1
e′2
T
(xi − xQ ).
(21)
With such treatment, the surface smoothing can be done by taking x′ (x′1 , x′2 ) as the spatial variable and x(x1 , x2 , x3 ) as the field function for interpolation. On repeating the point interpolation process described above for each coordinate component xi , i = 1, 2, 3, the equation for a smoothed local contact surface can be written as x(x′ ) = ΦT (x′ )χ = ΦT1 (x′ )χ1
ΦT2 (x′ )χ2
ΦT3 (x′ )χ3
T
(22)
where χi = [ ] is a vector composed of the xi components of all the master nodes in the domain of influence. ··· Using this preconditioning approach, the overlap problem can be avoided in most situations. A few examples are shown in Fig. 5. For situations (a)–(c), the overlap problems are eliminated. However, there is still overlap for situation (d) and a possible remedy is using a smaller domain of influence. x1i
x2i
xni T
3.3. Surface smoothing examples Here two examples are included to illustrate the performance of the proposed RPIM surface smoothing method. Two local surface meshes resulting from FE discretization of a sphere and a saddle are shown Fig. 6. Their analytical expressions
6
X. Qian et al. / Journal of Computational and Applied Mathematics 257 (2014) 1–13
Fig. 5. Preconditioning examples.
Fig. 6. Local FE surface meshes: (a) spherical surface and (b) saddle surface.
Fig. 7. Projection of the local FE surface meshes onto the x1 x2 plane.
Fig. 8. Smoothed surfaces: (a) spherical surface and (b) saddle surface.
are given by (a) Spherical surface: (x1 − 2.2731)2 + (x2 − 1.0495)2 + x23 = 5.12 (b) Saddle surface: x3 = 15(x1 − 1.6843)2 /14 − (x2 − 1.1113)2 /2.
(23)
Their projections onto the x1 x2 plane are the same, as shown in Fig. 7, consisting of both triangular and quadrilateral facets. The constructed surfaces, as displayed in Fig. 8, are smooth both inside and across patches and pass through all the surface nodes. To evaluate the deviation of numerical descriptions from the analytic surface, the coordinate error and the error of normal vector are defined as
e0 =
x˜ − x ds
ds,
en =
˜ · n)ds arccos(n
ds
(24)
X. Qian et al. / Journal of Computational and Applied Mathematics 257 (2014) 1–13
7
Table 1 Interpolation errors of different methods. Surfaces
Coordinate error
Error of normal vector
FE
RPIM
FE
RPIM 0.01426 (0.817°) 0.05004 (2.867°)
Spherical surface
0.02937
0.00393
0.07957 (4.559°)
Saddle surface
0.03207
0.01464
0.18033 (10.332°)
where x(x1 , x2 , x3 ) and n denote the coordinates and unit outward normal calculated from the analytical expression, ˜ denote their numerical approximations through either FE interpolation or RPIM smoothing. The and x˜ (˜x1 , x˜ 2 , x˜ 3 ) and n interpolation errors of the two methods are summarized in Table 1, where RPIM shows much smaller interpolation errors, especially in the approximation of the normal vector that is particularly important for contact analysis. 4. Contact search and adjustment of the time step This section describes the contact search process following the surface smoothing and the adjustment of time step to reduce interpenetration. 4.1. Location of the projection point and gap calculation On the basis of the smoothed master contact surface γS in the vicinity of the slave node xS given by (22), the projection point can be obtained by solving the following system of equations: H (x′ ) =
H1 H2
=
(xs − x) · x,1 =0 (xs − x) · x,2
(25)
where x,1 , x,2 are the directional derivatives of x with respect to x′ . This system of equations can be solved iteratively using the Newton–Raphson method and x′0 = [0 0]T can be conveniently taken as the initial value since the temporary projection point xQ is the origin of the 2D natural coordinate system (Fig. 4), as indicated by (21). On denoting the solution of (25) by x′P , the projection point P and corresponding unit outward normal are given by xP = x(x′P ),
nP =
xP ,1 × xP ,2 ∥xP ,1 × xP ,2 ∥
(26)
where xP ,1 , xP ,2 are the directional derivatives of xP . With xP and nP in place, the gap function t gN can be calculated from (2), according to which, the contact status between the slave node xS and the corresponding local contact surface γS can then be determined. t gN > 0 means gap, t gN < 0 means interpenetration, and t gN = 0 means that the slave node lies on the local contact surface. In practical contact searching algorithms, threshold values are usually used to consider the precision of numerical computation. 4.2. Contact tracking and adjustment of the time step For contact bodies with a contact history, the contact searching procedure is different from that of newly activated contact bodies. In a time step t → t + ∆t, configurations of the contact bodies will change with deformation and/or motion. Thus the contact information needs to be updated. However, a global search is usually unnecessary since the contact information at time t is known. The temporary contact facet t +∆t MF and temporary projection point t +∆t xQ can be obtained through a local tracking process, as described in [6]. Then, a locally smooth contact surface t +∆t γS can be constructed using the RPIM smoothing method described in Section 3, from which the projection point t +∆t xP at time t + ∆t is located and the gap function t +∆t gN is calculated. Due to incremental calculation, the contacting points may interpenetrate at time t + ∆t. Theoretically, this situation can be eliminated by decreasing the time step. However, complete elimination of interpenetration may lead to a very small time step and thus a significant increase in computation time for complex contact problems. On the other hand, excessive interpenetration, if left untreated, will impair the accuracy of the numerical solution. Here an allowed interpenetration gp , specified by the user, is introduced to deal with this contradiction. The projection point t +∆t xP will be taken as the final contacting point if t +∆t gN ≥ −gp . Otherwise, the time step has to be reduced. Assuming a linear response in time t → t + ∆t, a reduction coefficient is suggested as
λ=
tg r N
t r gN t +∆ t g N
−
(27)
8
X. Qian et al. / Journal of Computational and Applied Mathematics 257 (2014) 1–13
Fig. 9. Reduction of the time step.
Fig. 10. Contact between the foundation and the subsoil.
where tr gN is the residual gap of the last time step, as shown in Fig. 9. In case that tr gN < 0, (27) is replaced by
λ=
t gp r gN tg t +∆t g N N r
+
−
.
(28)
When there are multiple interpenetrated contacting points, the minimum reduction coefficient should be taken. A number of reductions may be required due to the nonlinearity of contact problems. Once all the contacting points (without interpenetration or with interpenetration smaller than the allowed value gp ) are identified, an iterative procedure [25] is employed to solve the frictional contact problem of this time step, during which a constraint is applied to the normal displacement of each slave node in contact to impose the impenetrability condition [29]. A contacting slave node is forced to be on the contact surface and may slide on it, subject to the friction conditions. Therefore, the interpenetration is removed. This treatment might introduce undesired stress changes that can be dealt with during the iterative procedure. 5. Numerical examples A general finite element analysis program with implicit time integration is developed, where an improved iterative procedure (formulated in terms of the relative displacement freedoms in a local coordinate system) based on the work in [25] is employed to solve the nonlinear system equations. The doubly connected edge list (DCEL) data structure [30] is adopted to facilitate the contact search and a direct constraint method [29] is used to incorporate the contact constraints in the variational formulation. To demonstrate the advantages of the proposed RPIM surface smoothing method in contact analysis, two numerical examples of 3D contact problems are presented and compared with the results of the traditional finite element method in this section. 5.1. Contact between the foundation and the subsoil The contact between a square foundation and the subsoil is considered. The upper structure is simplified as a uniformly distributed vertical load applied to the foundation. The 3D finite element model is shown in Fig. 10. The dimensions of the foundation are 6 m × 6 m × 0.6 m and its top surface is subjected to a normal pressure of 1 MPa. The subsoil is simulated as a block with dimensions 10 m × 10 m × 4 m, fixed on the bottom and horizontally confined on the four sides. Here both the foundation and the subsoil are assumed to be elastic to concentrate on the interface behavior. Their Young’s modulus and Poisson’s ratio are taken as E = 100 MPa, ν = 0.2. Friction is modeled using Coulomb’s law with µ = 0.1. In the numerical calculation, the foundation is chosen as the slave and the subsoil as the master. The element size of the quadrilateral mesh on the slave contact surface (bottom of the foundation) is 0.2 m and that on the master contact surface (top surface of the subsoil) is 0.4 m (Fig. 12(a)). The master contact surface is described by the original FE surface mesh and the proposed RPIM surface smoothing method respectively. Fig. 11 shows the distribution of the normal contact stress along the intersection between the contact interface and the center plane, from which it can be seen that the two contact surface description methods give similar results. The normal contact stress is relatively uniform and slightly smaller than the applied normal pressure in the middle part of the contact
X. Qian et al. / Journal of Computational and Applied Mathematics 257 (2014) 1–13
9
Fig. 11. Distribution of the normal contact stress (quadrilateral mesh).
Fig. 12. Meshes for the master contact surface: (a) quadrilateral mesh, (b) triangular mesh, and (c) hybrid mesh.
Fig. 13. Distributions of the normal contact stress: (a) triangular mesh and (b) hybrid mesh.
interface, while stress concentration occurs at the boundary of the contact interface. A further comparison of the results shows that the normal contact stress obtained using an FE surface mesh fluctuates clearly, especially at the boundary of the contact interface, whereas the result with RPIM smoothing is much smoother. To demonstrate the applicability of the proposed RPIM surface smoothing method to arbitrary surface meshes, calculations have been conducted using triangular and hybrid meshes. The meshes for the master contact surface are shown in Fig. 12(b) and (c) and the results are shown in Fig. 13. As can be seen from the display, the RPIM smoothing method produces smoother results in both cases. 5.2. Cylindrical contactor sliding in a half-tube This example simulates the sliding of a cylindrical contactor in a half-tube [15]. The 3D FE model is shown in Fig. 14. The half-tube has the inner radius of 3 m and thickness of 1 m. It is 12 m long and fixed at the bottom. Its Young’s modulus and Poisson’s ratio are E2 = 30 MPa and ν2 = 0.2 respectively. The cylindrical contactor, with dimensions 4 m × 4 m × 1 m, has the curvature radius of 3 m in the contact surface. It is characterized by Young’s modulus E1 = 10 MPa and Poisson’s ratio ν1 = 0.2. Under a uniformly distributed normal pressure of 1 MPa applied to its top surface, the contactor is in tight contact with the half-tube. Displacements are applied to the rear surface of the contactor to make it slide in the half-tube at a speed of 0.1 m/time step. The frictional behavior between them is modeled using Coulomb’s law with µ = 0.1. Both the contactor and the half-tube are discretized by hexahedral elements, resulting in quadrilateral surface meshes. In the numerical analysis, the contactor is designated as the slave and the half-tube as the master. The element size on the contactor surface is about 1/3 m × 0.25 m and that on the surface of the half-tube is about 1 m × 0.75 m. The FE discretization of the master contact surface is relatively coarse and thus is very demanding on the smoothing procedure. 30 time steps of sliding are simulated. Under applied displacements, the contactor has slid by 3 m, i.e. 3 elements. The master
10
X. Qian et al. / Journal of Computational and Applied Mathematics 257 (2014) 1–13
Fig. 14. FE model for the cylindrical contactor — half-tube sliding simulation.
Fig. 15. Total friction for FE mesh and RPIM smoothing.
contact surface is described by the original FE surface mesh and the proposed RPIM surface smoothing method respectively for comparison. The total friction on the bottom of the contactor is calculated. It equals the total reaction on the rear surface of the contactor where displacements are applied. The results using the two contact surface description methods are compared in Fig. 15. The total friction — time (displacement) curve obtained using FE surface mesh displays clear oscillation. The period is about 10 time steps, i.e. 1 m (the length of an element on the half-tube). The result with RPIM smoothing is much smoother, which indicates that the proposed contact surface smoothing method is capable of eliminating the rough response that characterizes the solution using FE surface meshes. Another observation is that the friction force is smaller when the RPIM is adopted. The reason might be that the contact described by the FE meshes is actually between two polyhedrons. Due to the presence of ridges, it is possible that the relative sliding of two polyhedrons under compression will induce larger normal stress and friction force. This phenomenon is not remarkable for the contact between the foundation and the subsoil (Fig. 11) since their contact surface is quite flat. Distributions of the normal contact stress along the curves specified in Fig. 16 at the last time step are shown in Fig. 17. Comparable results have been obtained using the two contact surface description methods. Along the transverse curves F, TM, and R, the normal contact stress distribution is U-shaped (small at the center and large at both sides), reflecting the stress concentration phenomenon. Along the longitudinal curves L and LM, the normal stress on the contact interface is large in the front and small in the middle and rear, indicating that the front part of the contactor’s bottom surface is subjected to a larger friction since friction is proportional to normal stress. Moreover, comparison between the results shows that the RPIM smoothing method produces stress distributions much smoother than the FE surface mesh description. 5.3. Convergence analysis To compare the convergence of both methods, i.e., the FE surface mesh and RPIM smoothing, the iterative process for the initial step (application of the normal pressure) of the contactor sliding problem is analyzed. A personal computer with R Intel⃝ CoreTM i5 CPU 760 @ 2.80 GHz, 8.0 GB RAM and ST3500418AS 7200 RPM Hard Drive is used to perform the required numerical analysis. Fig. 18(a) shows the variation of the residual during the iterative process, where residual refers to the L2 norm of the residual nodal force vector. With the RPIM smoothing method, the number of iterations is significantly reduced. Fig. 18(b) plots the number of nodes with different contact status at each iteration, from which it can be seen that the contact status of nodes stabilize faster when the RPIM smoothing method is used. The CPU time with the FE surface mesh is 12.38 s while the RPIM smoothing method takes a much shorter time, 5.086 s. However, when the contact surface is quite flat and the deformation is relatively small, the number of iterations using both methods may be comparable and the RPIM smoothing method may take a longer time due to the additional time required for contact surface smoothing. For example,
X. Qian et al. / Journal of Computational and Applied Mathematics 257 (2014) 1–13
11
Fig. 16. Curves for stress distribution plotting.
Fig. 17. Distributions of the normal contact stress along (a) the front curve F, (b) the transverse middle curve TM, (c) the rear curve R, (d) the left curve L, and (e) the longitudinal middle curve LM.
the number of iterations using both methods for the foundation and the subsoil example is the same, 9. The CPU times are 18.454 and 18.888 min for FE surface mesh and RPIM smoothing respectively. On the basis of its formulation, the RPIM smoothing method may also have the potential of overcoming the ‘chatter’ problem, often encountered when solving contact problems with FE surface meshes. A different FE mesh for the contactor sliding problem is shown in Fig. 19, where the mesh line Y = 2.5 m on the bottom of the contactor coincides with a mesh line on the inner surface of the half-tube. The variations of the residual and the number of nodes with different contact status at different iterations using both methods are displayed in Fig. 20. Due to the coincidence of mesh lines and the symmetry of this problem, slave nodes on the bottom line Y = 2.5 m on the contactor ‘jump’ between neighboring facets on the inner
12
X. Qian et al. / Journal of Computational and Applied Mathematics 257 (2014) 1–13
Fig. 18. Convergence analysis: (a) residual and (b) contact status.
Fig. 19. A different FE mesh for the contactor sliding problem.
Fig. 20. Convergence analysis (a different FE mesh): (a) residual and (b) contact status.
surface of the half-tube during the iterative process and hence impede the convergence when the FE surface mesh method is used. In contrast, it is clear from Fig. 20 that the RPIM smoothing method does not suffer from this problem. Although special treatments, e.g., allowing a slave node to contact with a master node instead of a contact facet, can be employed to deal with the particular problem encountered here, the proposed RPIM smoothing method provides a general way to avoid the ‘chatter’ problem. 6. Conclusions A general method for 3D contact surface smoothing has been proposed. This method is based on creating a smooth description of the local contact surface using the meshfree radial point interpolation method (RPIM). In particular, radial and polynomial basis functions are used to interpolate the scattered surface nodes and thus construct a smooth contact surface in the vicinity of any slave node likely to come into contact. With such a definition, this method applies to arbitrary surface meshes and is capable of reproducing smooth surfaces passing exactly through the surface nodes even for coarse meshes. The surface smoothing algorithm has been tested on different 3D surfaces. It shows a significant improvement in representing the real geometry of 3D bodies, especially for normal vector approximation. To demonstrate its advantages
X. Qian et al. / Journal of Computational and Applied Mathematics 257 (2014) 1–13
13
in contact analysis, numerical simulations of contact between deformable bodies have been performed, where adjustment of time step is employed to control the interpenetration. The results indicate that the proposed contact surface smoothing method can effectively reduce the oscillation of contact forces and non-smoothness of stress distributions. Through convergence analysis, the proposed RPIM smoothing method demonstrates a great potential of improving the convergence and overcoming the ‘chatter’ problem for contact analyses involving curved contact surfaces. Acknowledgments The authors would like to thank the financial support of the National Basic Research Program of China #2010CB732103, the National Natural Science Foundation of China #51209118, and the State Key Laboratory of Hydroscience and Engineering #2012-KY-02. References [1] D.J. Benson, J.O. Hallquist, A single surface contact algorithm for the post-buckling analysis of shell structures, Comput. Methods Appl. Mech. Eng. 91 (1990) 141–163. [2] J.O. Hallquist, G.L. Goudreau, D.J. Benson, Sliding interfaces with contact-impact in large-scale Lagrangian computations, Comput. Methods Appl. Mech. Eng. 51 (1985) 107–137. [3] Z.H. Zhong, L. Nilsson, A unified contact algorithm based on the territory concept, Comput. Methods Appl. Mech. Eng. 130 (1996) 1–16. [4] T. Belytschko, J.I. Lin, A three-dimensional impact-penetration algorithm with erosion, Comput. Struct. 25 (1987) 95–104. [5] T. Belytschko, M.O. Neal, Contact-impact by the pinball algorithm with penalty projection and Lagrangian methods, Internat. J. Numer. Methods Engrg. 31 (1991) 547–572. [6] S.P. Wang, E. Nakamachi, The inside-outside search algorithm for finite element analysis, Internat. J. Numer. Methods Engrg. 40 (1997) 3665–3685. [7] G. Zavarise, L. De Lorenzis, The node-to-segment algorithm for 2D frictionless contact: classical formulation and special cases, Comput. Methods Appl. Mech. Eng. 198 (2009) 3428–3451. [8] V. Padmanabhan, T.A. Laursen, A framework for development of surface smoothing procedures in large deformation frictional contact analysis, Finite Elem. Anal. Des. 37 (2001) 173–198. [9] P. Wriggers, L. Krstulovic-Opara, J. Korelc, Smooth C 1 -interpolations for two-dimensional frictional contact problems, Internat. J. Numer. Methods Engrg. 51 (2001) 1469–1495. [10] S. Ulaga, M. Ulbin, J. Flasker, Contact problems of gears using overhauser splines, Int. J. Mech. Sci. 41 (1999) 385–395. [11] N. El-Abbasi, S.A. Meguid, A. Czekanski, On the modelling of smooth contact surfaces using cubic splines, Internat. J. Numer. Methods Engrg. 50 (2001) 953–996. [12] M. Stadler, G.A. Holzapfel, J. Korelc, C n continuous modelling of smooth contact surfaces using NURBS and application to 2D problems, Internat. J. Numer. Methods Engrg. 57 (2003) 2177–2203. [13] G. Pietrzak, A. Curnier, Large deformation frictional contact mechanics: continuum formulation and augmented Lagrangian treatment, Comput. Methods Appl. Mech. Eng. 177 (1999) 351–381. [14] F.J. Wang, J.G. Cheng, Z.H. Yao, FFS contact searching algorithm for dynamic finite element analysis, Internat. J. Numer. Methods Engrg. 52 (2001) 655–672. [15] L. Krstulovic-Opara, P. Wriggers, J. Korelc, A C 1 -continuous formulation for 3D finite deformation frictional contact, Comput. Mech. 29 (2002) 27–42. [16] M.A. Puso, T.A. Laursen, A 3D contact smoothing method using Gregory patches, Internat. J. Numer. Methods Engrg. 54 (2002) 1161–1194. [17] D.M. Neto, M.C. Oliveira, L.F. Menezes, J.L. Alves, Improving Nagata patch interpolation applied for tool surface description in sheet metal forming simulation, Comput. Aided Des. 45 (2013) 639–656. [18] T. Hama, T. Nagata, C. Teodosiu, A. Makinouchi, H. Takuda, Finite-element simulation of springback in sheet metal forming using local interpolation for tool surfaces, Int. J. Mech. Sci. 50 (2008) 175–192. [19] D.M. Neto, M.C. Oliveira, L.F. Menezes, J.L. Alves, Nagata patch interpolation using surface normal vectors evaluated from the IGES file, Finite Elem. Anal. Des. 72 (2013) 35–46. [20] D. Chamoret, P. Saillard, A. Rassineux, J.M. Bergheau, New smoothing procedures in contact mechanics, J. Comput. Appl. Math. 168 (2004) 107–116. [21] T. Belytschko, W.J.T. Daniel, G. Ventura, A monolithic smoothing-gap algorithm for contact-impact based on the signed distance function, Internat. J. Numer. Methods Engrg. 55 (2002) 101–125. [22] J.G. Wang, G.R. Liu, A point interpolation meshless method based on radial basis functions, Internat. J. Numer. Methods Engrg. 54 (2002) 1623–1648. [23] T.A. Laursen, J.C. Simo, A continuum-based finite element formulation for the implicit solution of multibody large deformation frictional contact problems, Internat. J. Numer. Methods Engrg. 36 (1993) 3451–3485. [24] M.C. Wang, Finite Element Method, Tsinghua University Press, Beijing, 2003. [25] M.U. Rahman, R.E. Rowlands, R.D. Cook, T.L. Wilkinson, An iterative procedure for finite-element stress analysis of frictional contact problems, Comput. Struct. 18 (1984) 947–954. [26] A.B. Chaudhary, K.J. Bathe, A solution method for static and dynamic analysis of three-dimensional contact problems with friction, Comput. Struct. 24 (1986) 855–873. [27] I. Doghri, A. Muller, R.L. Taylor, A general three-dimensional contact procedure for implicit finite element codes, Eng. Comput. 15 (1998) 233–259. [28] M. Oldenburg, L. Nilsson, The position code algorithm for contact searching, Internat. J. Numer. Methods Engrg. 37 (1994) 359–386. [29] MSC. Software Corporation, Marc Volume A: Theory and User Information, MSC. Software Corporation, 2010, pp. 554–556. [30] M. de Berg, O. Cheong, M. van Kreveld, M. Overmars, Computational Geometry: Algorithms and Applications, Springer, Berlin, 2008.
Contents lists available at ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
A general 3D contact smoothing method based on radial point interpolation Xiaoxiang Qian, Huina Yuan, Mozhen Zhou, Bingyin Zhang ∗ State Key Laboratory of Hydroscience and Engineering, Department of Hydraulic Engineering, Tsinghua University, Beijing 100084, China
highlights • • • • •
A general 3D contact smoothing method using radial point interpolation is proposed. The method constructs local contact surfaces with polynomial and radial bases. The constructed surface is smooth and passes exactly through the surface nodes. The method shows good approximation in both coordinates and normal vector. It can effectively reduce contact force oscillation and solve the chatter problem.
article
info
Article history: Received 9 May 2013 Received in revised form 4 August 2013 Keywords: Contact Finite element method Surface smoothing Radial point interpolation method
abstract This paper presents a general 3D contact smoothing method based on the meshfree radial point interpolation method to improve the numerical simulation of contact problems. In particular, a locally smooth contact surface is constructed from the scattered surface nodes by point interpolation using the combination of polynomial and radial bases. With such bases, this method reproduces smooth surfaces even for coarse meshes and the constructed surface passes exactly through the surface nodes. Results for contact problems involving deformable bodies are included to demonstrate its advantages. © 2013 Elsevier B.V. All rights reserved.
1. Introduction Significant progress has been made on the finite element (FE) analysis of contact problems. However, an accurate and efficient representation of the nonlinear geometry of contacting surfaces remains a challenge. When a standard master–slave approach is used, piecewise continuous FE discretization of the master contact surface may cause two types of problems. One type of problems is associated with the searching of contact segments. With the commonly used node-to-segment (NTS) algorithm, a slave node may have multiple or no candidates of contact segments due to C 0 continuity of the master surface [1,2]. To circumvent these problems, alternative contact search algorithms have been developed, such as the hierarchyterritory algorithm [3], the pinball algorithm [4,5], and the inside–outside algorithm [6]. Recently, some new techniques have also been proposed to deal with these special situations within the framework of the NTS algorithm [7]. The other type of problems, often met in large motion/deformation problems, is associated with the jump of normal vector when a slave node slides from one segment to an adjacent segment on the master contact surface. These jumps, induced by spatial discretization, are non-physical and can cause serious errors in the calculation of contact forces. Moreover, slave nodes situated near an edge may ‘chatter’ during equilibrium iterations and thus impede convergence. To overcome these problems, a general idea is to smooth the master surface. For 2D contact, spline-based smoothing techniques have typically been used, such as the cubic Hermite interpolation [8], the cubic Bézier spline [9], the Overhauser spline [10,11]
∗
Corresponding author. Tel.: +86 10 62787349; fax: +86 10 62785593. E-mail address: [email protected] (B. Zhang).
0377-0427/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cam.2013.08.014
2
X. Qian et al. / Journal of Computational and Applied Mathematics 257 (2014) 1–13
Fig. 1. Different types of surface meshes: (a) regular quadrilateral mesh, (b) triangular mesh, (c) irregular quadrilateral mesh, and (d) hybrid mesh.
Fig. 2. Notation for the frictional contact problem.
and the non-uniform rational B-splines [12]. Unfortunately, an extension of such methods to 3D arbitrary meshes proves to be much more difficult. Fig. 1 shows different types of surface meshes for 3D contact problems. 3D generalizations of 2D spline curves through tensor product only work for regular quadrilateral meshes (Fig. 1(a)) where each inside quadrilateral facet is surrounded by eight neighboring quadrilateral facets [13,14]. Smooth descriptions using quartic Bézier surfaces [15] were proposed for contact surfaces composed of triangular facets particularly (Fig. 1(b)). Another 3D approach [16] uses Gregory patches to interpolate irregular quadrilateral meshes (Fig. 1(c)) where a quadrilateral facet has an arbitrary number of neighboring quadrilateral facets. The Gregory patch augments the third-order Bézier patch with additional internal control points to enable tangent plane continuity in general. With special treatments [17] developed for patches involving inflection points, the Nagata patch interpolation algorithm can be applied to both triangular and quadrilateral faceted finite elements [18,17,19]. However, a consensual good solution for contact smoothing of complex surfaces composed of hybrid meshes (Fig. 1(d)), where each facet, surrounded by an arbitrary number of neighboring facets, can be either quadrilateral or triangular, remains a challenge. Due to the geometrical complexity of 3D bodies, hybrid surface meshes are often encountered when mixed 3D element types, e.g. tetrahedral and hexahedral elements, are used. Chamoret et al. [20] proposed a smoothing procedure based on a meshfree technique denoted as diffuse approximation. This approach creates a locally smooth contact surface from the scattered surface nodes using a least-squares approximation with a polynomial basis. Belytschko et al. [21] proposed a similar approach, where the smoothing is done implicitly by constructing smooth signed distance functions for the bodies. Without using element connectivity, these methods are applicable to arbitrary surface meshes, such as those shown in Fig. 1. However, the constructed contact surface does not pass through the surface nodes exactly, which introduces some inaccuracies with regard to the real geometry of contact bodies. The goal of this work is a general 3D contact smoothing method based on a meshfree technique denoted as the radial point interpolation method (RPIM) [22]. The basic concept is to construct a locally smooth contact surface from the scattered surface nodes using point interpolation with both polynomial and radial basis functions. Similar to the methods proposed in [20,21], this approach also applies to arbitrary surface meshes. Moreover, the smoothed contact surface passes exactly through all the surface nodes of the master body in the vicinity of any slave node where contact is considered. 2. Notation and preliminaries This section briefly reviews the fictional contact problem and the notion to be used in this paper [23,24]. Consider a frictional contact problem depicted in Fig. 2, where the number of potential contact bodies is limited to two for convenience. The initial configurations of the two bodies (at time t = 0) are represented by open sets 0 Ω 1 , 0 Ω 2 ⊂ R3 , where X 1 ∈ 0 Ω 1 , X 2 ∈ 0 Ω 2 are two material points. The motions of the two bodies are defined through mappings: x = t ϕ k (X k )
t k
k = 1, 2
(1)
where x , x are the current positions of X , X (at time t ≥ 0) and Ω , Ω are the current configurations of the two bodies. Although these two bodies are interchangeable with respect to their mechanics, one (t Ω 2 ) is designated as the master and the other (t Ω 1 ) as the slave to describe their contact. Correspondingly, the potential contacting surfaces of the two bodies, t γ 2 and t γ 1 , are designated as the master and slave contact surfaces respectively. t 1 t 2
1
2
t
1 t
2
X. Qian et al. / Journal of Computational and Applied Mathematics 257 (2014) 1–13
3
The bodies are governed by the conservation of mass, momentum and energy, the constitutive relations, and the initial and boundary conditions as usual. In addition, contact between the bodies is subjected to the normal and frictional contact constraints. First, the bodies cannot interpenetrate, called the impenetrability condition. This condition can be written as t
gN (t x1S ) = t n2P · (t x1S − t x2P )
= min t x1S − t x2 sign(t n2 · (t x1S − t x2 )) ≥ 0 ∀t x1S ∈ t γ 1
(2)
t x2 ∈t γ 2
where t x2P is the projection (closest) point of a slave point t x1S ∈ t γ 1 onto the master surface t γ 2 , t n2 is the unit outward normal to the master t Ω 2 , and t gN is the gap function. In the example shown in Fig. 2, the points of the contacting surfaces have interpenetrated (a situation often encountered during numerical simulation) and thus t gN < 0. Second, it is assumed that the normal traction across the contact interface cannot be tensile, i.e. t
FN1 = −t FN2 = t F 1 · t n2 = −t F 2 · t n2 ≥ 0
(3)
where t FN1 , t FN2 are the normal tractions on t γ 1 , t γ 2 respectively. Finally, the frictional behavior of the contact interface follows a specified friction model. Here the Coulomb friction model is considered:
¯ T = t υ1T − t υ2T = 0 ⇔ t FT1 < µ t FN1 Stick: t υ
(4)
υ¯ T = υ − υ ̸= 0 ⇔ t FT1 = µ t FN1 Slip: t t υ¯ T · t FT1 < 0, υ¯ T × t FT1 = 0
(5)
t
t
1 T
t
2 T
¯ T is the relative tangential velocity between a pair of contacting points, t FT1 = t F 1 − t FN1 t n2 is the tangential where t υ traction, and µ is the coefficient of friction. The last inequation and equation in (5) represent that the tangential velocity and tangential contact force are collinear. Within the framework of an incremental solution method, the weak form of the contact problem at time t + ∆t (assuming that the solution at time t has been obtained) can be written as 1 ,2
t +∆t
k δ Wint − t +∆t δ WLk − t +∆t δ WIk − t +∆t δ WCk = 0
(6)
k=
where δ Wint , δ WL , δ WI , and δ WC are the internal energy, the work of external forces, the work of inertial forces, and the work of contact forces respectively. In particular, the work of contact forces is given by t +∆ t
δ WC =
1,2
t +∆t
δ WCk =
k=
1 ,2 k=
t +∆t t +∆t γ k
Fik δ uki t +∆t ds =
t +∆t t +∆t γ 1
Fi1 δ ui t +∆t ds
(7)
where δ u is the variation of the relative displacement between a pair of contacting points. The impenetrability condition can also be written in an incremental form as t +∆ t
gN = t +∆t n2 · (t +∆t x1 − t +∆t x2 ) = t +∆t n2 · (u1 − u2 ) + t +∆t n2 · (t x1 − t x2 )
= u1N − u2N + tr gN ≥ 0
(8)
where u1N , u2N are the normal displacement increments of a pair of contacting points over ∆t time and tr gN is the residual gap from the last time step. Accordingly, the tangential velocity υ T in (4) and (5) can be replaced by the tangential displacement increment uT . The details of solving (6) by way of the finite element method can be found in [25–27]. The main concern of this paper is contact surface smoothing that plays a significant role in accurate location of the projection point and continuity of the normal vector. 3. Contact surface smoothing The contact surface smoothing procedure based on the meshfree radial point interpolation method (RPIM) is presented in this section. Resulting from FE discretization of the contact bodies, the contacting surfaces are originally represented by collections of quadrilateral and/or triangular facets. When a master–slave approach is used, only the master contact surface requires smoothing. The interpenetration of any slave node can be computed in terms of the gap function. 3.1. Surface smoothing based on RPIM Consider a slave node S (xS ∈ γ 1 ) for contact searching. For simplicity, the superscripts representing time and body are omitted hereinafter when there is no ambiguity. The basic idea is to create a smooth description of the contact surface using surface nodes of the master body Ω 2 in the vicinity of xS . Thus, it is first necessary to find the surface nodes on γ 2 in the neighborhood of xS . For two newly activated contact bodies, this can be done through a global search and a preliminary
4
X. Qian et al. / Journal of Computational and Applied Mathematics 257 (2014) 1–13
Fig. 3. A domain of influence for surface smoothing.
local search. Here the position code algorithm [28] and the inside–outside algorithm [6] are used to locate the temporary contact facet MF and the temporary projection point Q (xQ ). The vicinity of xS , also referred to as the domain of influence in the context of RPIM, can then be determined via a multilevel neighboring search for MF. A domain of influence for surface smoothing is shown in Fig. 3 as an example, where the temporary contact facet MF is shaded. Next, the meshfree radial point interpolation method [22] is applied to construct a locally smooth contact surface γS over the domain of influence ΩS of xS , using the information (coordinates) of all the surface nodes in ΩS (Fig. 3). Let (xi1 , xi2 , xi3 ) denote the coordinates of a master node xi ∈ ΩS in a reference frame Rd . The local approximation of the contact surface for the slave node xS is described by the following equation: x3 = fS (x1 , x2 ) =
n
Ri (x)ai +
i=1
m
Pj (x)bj = R T (x)a + P T (x)b
(9)
j=1
where n is the number of master nodes in ΩS , Ri , i = 1, . . . , n is a term of the radial basis, Pj , j = 1, . . . , m is a term of the polynomial basis, and ai , i = 1, . . . , n and bj , j = 1, . . . , m are corresponding coefficients. The radial basis function Ri takes the form: Ri (x) = Ri (x − xi ) = [ri2 + (Rdc )2 ]q
(10)
(x1 − xi1 )2 + (x2 − xi2 )2 is the distance between x and xi projected onto the x1 x2 plane, dc is the average where ri = projected node distance in ΩS , and R, q are selected shape parameters. As for the polynomial basis, commonly used linear and quadratic bases are given by P T (x) = [1
x2 ]T
x1
P (x) = [1 T
x1
x2
m=3 x1 x1
(11)
x1 x2
T
x2 x2 ]
m = 6.
(12)
The coefficient vectors a and b can be determined by making the surface pass through all the master nodes in ΩS , which leads to the following system of equations:
χ 0
=
C DT
where χ = [x13
x23
R1 (x1 ) R1 (x2 )
C =
.. . R1 (xn )
D 0
a a =G b b
···
(13)
xn3 ]T is a vector of the x3 components of all the master nodes xi ∈ ΩS , and
R2 (x1 ) R2 (x2 )
.. . R2 (xn )
··· ··· .. . ···
Rn (x1 ) Rn (x2 )
.. , . Rn (xn )
1 1
D = . .
. 1
x11 x21
x12 x22
xn1
xn2
.. .
.. .
··· ··· .. . ···
Pm (x1 ) Pm (x2 )
.. . Pm (xn )
(14)
are matrices formed by the radial and polynomial basis functions. With a and b solved from (13), the x3 component of any point x on the constructed local contact surface can be computed from x3 = fS (x1 , x2 ) = R T (x)
a
P T (x)
= R T (x) b
−1 P T (x) G
χ 0
= ΦT (x)χ
(15)
where Φ(x) is a vector of resulting RPIM shape functions. From the forms of the radial and polynomial basis functions, the constructed surface γS is C n continuous both inside and across patches whereby n represents an arbitrary level of continuity. In particular, the RPIM shape functions have the Kronecker delta property, i.e.
Φi (xj ) =
1, 0,
i=j i ̸= j
which ensures that the constructed surface passes exactly through all the master nodes in the domain of influence.
(16)
X. Qian et al. / Journal of Computational and Applied Mathematics 257 (2014) 1–13
5
Fig. 4. The 2D natural coordinate system for point interpolation.
3.2. Preconditioning Further examination of (13) shows that the interpolation will fail if projection of the local contact surface in the vicinity of xS onto the x1 x2 plane overlaps. To prevent this problem, a preconditioning method is proposed. Let point Q in the surface Q Q element MF be the control point of the domain of influence, as shown in Fig. 3. Its natural coordinates ξ Q (ξ1 , ξ2 ) in the element and global coordinates xQ (x1 , x2 , x3 ) in the reference frame Rd are related by Q
ne
xQ =
Q
Q
Nα (ξ Q )xα
(17)
α=1
where ne is the number of nodes in the surface element (either quadrilateral or triangular), xα , α = 1, . . . , ne is the α -th node, and Nα is the corresponding shape function, depending on the type of the element. A suitable 2D natural coordinate system for point interpolation can be obtained by rotating the plane parallel to the x1 x2 plane and passing through point Q by a certain angle. The procedure is described as follows. (1) Calculate the unit outward normal n at point Q : n=
τ1 × τ2 = ni ei ∥τ 1 × τ 2 ∥
(18)
where ei , i = 1, 2, 3, is the coordinate basis of Rd , and τ j = ∂ xQ /∂ξj , j = 1, 2, is the directional derivative of xQ with Q
respect to ξ (2) Determine the pivot r and rotation angle θ : Q j .
r =
e3 × n
∥e3 × n∥
,
θ = arccos(n · e3 ).
(19)
(3) Rotate the plane F0 (parallel to the x1 x2 plane and passing through point Q ) around r by angle θ . The resulting reference plane FR is shown in Fig. 4, with the coordinate basis given by
e′1
e′2 =
1
n21 + n22
n21 n3 + n22
n1 n2 (n3 − 1) −n1 (n21 + n22 )
n1 n2 (n3 − 1) n22 n3 + n21
−n2 (n21 + n22 )
n1 (n21 + n22 )
n2 (n21 + n22 ) [e1 , e2 ] .
(20)
n3 (n21 + n22 )
(4) Project the master nodes (interpolation points) onto the FR plane and calculate the natural coordinates x′i (x′1i , x′2i ): ′ x′i = e1
e′2
T
(xi − xQ ).
(21)
With such treatment, the surface smoothing can be done by taking x′ (x′1 , x′2 ) as the spatial variable and x(x1 , x2 , x3 ) as the field function for interpolation. On repeating the point interpolation process described above for each coordinate component xi , i = 1, 2, 3, the equation for a smoothed local contact surface can be written as x(x′ ) = ΦT (x′ )χ = ΦT1 (x′ )χ1
ΦT2 (x′ )χ2
ΦT3 (x′ )χ3
T
(22)
where χi = [ ] is a vector composed of the xi components of all the master nodes in the domain of influence. ··· Using this preconditioning approach, the overlap problem can be avoided in most situations. A few examples are shown in Fig. 5. For situations (a)–(c), the overlap problems are eliminated. However, there is still overlap for situation (d) and a possible remedy is using a smaller domain of influence. x1i
x2i
xni T
3.3. Surface smoothing examples Here two examples are included to illustrate the performance of the proposed RPIM surface smoothing method. Two local surface meshes resulting from FE discretization of a sphere and a saddle are shown Fig. 6. Their analytical expressions
6
X. Qian et al. / Journal of Computational and Applied Mathematics 257 (2014) 1–13
Fig. 5. Preconditioning examples.
Fig. 6. Local FE surface meshes: (a) spherical surface and (b) saddle surface.
Fig. 7. Projection of the local FE surface meshes onto the x1 x2 plane.
Fig. 8. Smoothed surfaces: (a) spherical surface and (b) saddle surface.
are given by (a) Spherical surface: (x1 − 2.2731)2 + (x2 − 1.0495)2 + x23 = 5.12 (b) Saddle surface: x3 = 15(x1 − 1.6843)2 /14 − (x2 − 1.1113)2 /2.
(23)
Their projections onto the x1 x2 plane are the same, as shown in Fig. 7, consisting of both triangular and quadrilateral facets. The constructed surfaces, as displayed in Fig. 8, are smooth both inside and across patches and pass through all the surface nodes. To evaluate the deviation of numerical descriptions from the analytic surface, the coordinate error and the error of normal vector are defined as
e0 =
x˜ − x ds
ds,
en =
˜ · n)ds arccos(n
ds
(24)
X. Qian et al. / Journal of Computational and Applied Mathematics 257 (2014) 1–13
7
Table 1 Interpolation errors of different methods. Surfaces
Coordinate error
Error of normal vector
FE
RPIM
FE
RPIM 0.01426 (0.817°) 0.05004 (2.867°)
Spherical surface
0.02937
0.00393
0.07957 (4.559°)
Saddle surface
0.03207
0.01464
0.18033 (10.332°)
where x(x1 , x2 , x3 ) and n denote the coordinates and unit outward normal calculated from the analytical expression, ˜ denote their numerical approximations through either FE interpolation or RPIM smoothing. The and x˜ (˜x1 , x˜ 2 , x˜ 3 ) and n interpolation errors of the two methods are summarized in Table 1, where RPIM shows much smaller interpolation errors, especially in the approximation of the normal vector that is particularly important for contact analysis. 4. Contact search and adjustment of the time step This section describes the contact search process following the surface smoothing and the adjustment of time step to reduce interpenetration. 4.1. Location of the projection point and gap calculation On the basis of the smoothed master contact surface γS in the vicinity of the slave node xS given by (22), the projection point can be obtained by solving the following system of equations: H (x′ ) =
H1 H2
=
(xs − x) · x,1 =0 (xs − x) · x,2
(25)
where x,1 , x,2 are the directional derivatives of x with respect to x′ . This system of equations can be solved iteratively using the Newton–Raphson method and x′0 = [0 0]T can be conveniently taken as the initial value since the temporary projection point xQ is the origin of the 2D natural coordinate system (Fig. 4), as indicated by (21). On denoting the solution of (25) by x′P , the projection point P and corresponding unit outward normal are given by xP = x(x′P ),
nP =
xP ,1 × xP ,2 ∥xP ,1 × xP ,2 ∥
(26)
where xP ,1 , xP ,2 are the directional derivatives of xP . With xP and nP in place, the gap function t gN can be calculated from (2), according to which, the contact status between the slave node xS and the corresponding local contact surface γS can then be determined. t gN > 0 means gap, t gN < 0 means interpenetration, and t gN = 0 means that the slave node lies on the local contact surface. In practical contact searching algorithms, threshold values are usually used to consider the precision of numerical computation. 4.2. Contact tracking and adjustment of the time step For contact bodies with a contact history, the contact searching procedure is different from that of newly activated contact bodies. In a time step t → t + ∆t, configurations of the contact bodies will change with deformation and/or motion. Thus the contact information needs to be updated. However, a global search is usually unnecessary since the contact information at time t is known. The temporary contact facet t +∆t MF and temporary projection point t +∆t xQ can be obtained through a local tracking process, as described in [6]. Then, a locally smooth contact surface t +∆t γS can be constructed using the RPIM smoothing method described in Section 3, from which the projection point t +∆t xP at time t + ∆t is located and the gap function t +∆t gN is calculated. Due to incremental calculation, the contacting points may interpenetrate at time t + ∆t. Theoretically, this situation can be eliminated by decreasing the time step. However, complete elimination of interpenetration may lead to a very small time step and thus a significant increase in computation time for complex contact problems. On the other hand, excessive interpenetration, if left untreated, will impair the accuracy of the numerical solution. Here an allowed interpenetration gp , specified by the user, is introduced to deal with this contradiction. The projection point t +∆t xP will be taken as the final contacting point if t +∆t gN ≥ −gp . Otherwise, the time step has to be reduced. Assuming a linear response in time t → t + ∆t, a reduction coefficient is suggested as
λ=
tg r N
t r gN t +∆ t g N
−
(27)
8
X. Qian et al. / Journal of Computational and Applied Mathematics 257 (2014) 1–13
Fig. 9. Reduction of the time step.
Fig. 10. Contact between the foundation and the subsoil.
where tr gN is the residual gap of the last time step, as shown in Fig. 9. In case that tr gN < 0, (27) is replaced by
λ=
t gp r gN tg t +∆t g N N r
+
−
.
(28)
When there are multiple interpenetrated contacting points, the minimum reduction coefficient should be taken. A number of reductions may be required due to the nonlinearity of contact problems. Once all the contacting points (without interpenetration or with interpenetration smaller than the allowed value gp ) are identified, an iterative procedure [25] is employed to solve the frictional contact problem of this time step, during which a constraint is applied to the normal displacement of each slave node in contact to impose the impenetrability condition [29]. A contacting slave node is forced to be on the contact surface and may slide on it, subject to the friction conditions. Therefore, the interpenetration is removed. This treatment might introduce undesired stress changes that can be dealt with during the iterative procedure. 5. Numerical examples A general finite element analysis program with implicit time integration is developed, where an improved iterative procedure (formulated in terms of the relative displacement freedoms in a local coordinate system) based on the work in [25] is employed to solve the nonlinear system equations. The doubly connected edge list (DCEL) data structure [30] is adopted to facilitate the contact search and a direct constraint method [29] is used to incorporate the contact constraints in the variational formulation. To demonstrate the advantages of the proposed RPIM surface smoothing method in contact analysis, two numerical examples of 3D contact problems are presented and compared with the results of the traditional finite element method in this section. 5.1. Contact between the foundation and the subsoil The contact between a square foundation and the subsoil is considered. The upper structure is simplified as a uniformly distributed vertical load applied to the foundation. The 3D finite element model is shown in Fig. 10. The dimensions of the foundation are 6 m × 6 m × 0.6 m and its top surface is subjected to a normal pressure of 1 MPa. The subsoil is simulated as a block with dimensions 10 m × 10 m × 4 m, fixed on the bottom and horizontally confined on the four sides. Here both the foundation and the subsoil are assumed to be elastic to concentrate on the interface behavior. Their Young’s modulus and Poisson’s ratio are taken as E = 100 MPa, ν = 0.2. Friction is modeled using Coulomb’s law with µ = 0.1. In the numerical calculation, the foundation is chosen as the slave and the subsoil as the master. The element size of the quadrilateral mesh on the slave contact surface (bottom of the foundation) is 0.2 m and that on the master contact surface (top surface of the subsoil) is 0.4 m (Fig. 12(a)). The master contact surface is described by the original FE surface mesh and the proposed RPIM surface smoothing method respectively. Fig. 11 shows the distribution of the normal contact stress along the intersection between the contact interface and the center plane, from which it can be seen that the two contact surface description methods give similar results. The normal contact stress is relatively uniform and slightly smaller than the applied normal pressure in the middle part of the contact
X. Qian et al. / Journal of Computational and Applied Mathematics 257 (2014) 1–13
9
Fig. 11. Distribution of the normal contact stress (quadrilateral mesh).
Fig. 12. Meshes for the master contact surface: (a) quadrilateral mesh, (b) triangular mesh, and (c) hybrid mesh.
Fig. 13. Distributions of the normal contact stress: (a) triangular mesh and (b) hybrid mesh.
interface, while stress concentration occurs at the boundary of the contact interface. A further comparison of the results shows that the normal contact stress obtained using an FE surface mesh fluctuates clearly, especially at the boundary of the contact interface, whereas the result with RPIM smoothing is much smoother. To demonstrate the applicability of the proposed RPIM surface smoothing method to arbitrary surface meshes, calculations have been conducted using triangular and hybrid meshes. The meshes for the master contact surface are shown in Fig. 12(b) and (c) and the results are shown in Fig. 13. As can be seen from the display, the RPIM smoothing method produces smoother results in both cases. 5.2. Cylindrical contactor sliding in a half-tube This example simulates the sliding of a cylindrical contactor in a half-tube [15]. The 3D FE model is shown in Fig. 14. The half-tube has the inner radius of 3 m and thickness of 1 m. It is 12 m long and fixed at the bottom. Its Young’s modulus and Poisson’s ratio are E2 = 30 MPa and ν2 = 0.2 respectively. The cylindrical contactor, with dimensions 4 m × 4 m × 1 m, has the curvature radius of 3 m in the contact surface. It is characterized by Young’s modulus E1 = 10 MPa and Poisson’s ratio ν1 = 0.2. Under a uniformly distributed normal pressure of 1 MPa applied to its top surface, the contactor is in tight contact with the half-tube. Displacements are applied to the rear surface of the contactor to make it slide in the half-tube at a speed of 0.1 m/time step. The frictional behavior between them is modeled using Coulomb’s law with µ = 0.1. Both the contactor and the half-tube are discretized by hexahedral elements, resulting in quadrilateral surface meshes. In the numerical analysis, the contactor is designated as the slave and the half-tube as the master. The element size on the contactor surface is about 1/3 m × 0.25 m and that on the surface of the half-tube is about 1 m × 0.75 m. The FE discretization of the master contact surface is relatively coarse and thus is very demanding on the smoothing procedure. 30 time steps of sliding are simulated. Under applied displacements, the contactor has slid by 3 m, i.e. 3 elements. The master
10
X. Qian et al. / Journal of Computational and Applied Mathematics 257 (2014) 1–13
Fig. 14. FE model for the cylindrical contactor — half-tube sliding simulation.
Fig. 15. Total friction for FE mesh and RPIM smoothing.
contact surface is described by the original FE surface mesh and the proposed RPIM surface smoothing method respectively for comparison. The total friction on the bottom of the contactor is calculated. It equals the total reaction on the rear surface of the contactor where displacements are applied. The results using the two contact surface description methods are compared in Fig. 15. The total friction — time (displacement) curve obtained using FE surface mesh displays clear oscillation. The period is about 10 time steps, i.e. 1 m (the length of an element on the half-tube). The result with RPIM smoothing is much smoother, which indicates that the proposed contact surface smoothing method is capable of eliminating the rough response that characterizes the solution using FE surface meshes. Another observation is that the friction force is smaller when the RPIM is adopted. The reason might be that the contact described by the FE meshes is actually between two polyhedrons. Due to the presence of ridges, it is possible that the relative sliding of two polyhedrons under compression will induce larger normal stress and friction force. This phenomenon is not remarkable for the contact between the foundation and the subsoil (Fig. 11) since their contact surface is quite flat. Distributions of the normal contact stress along the curves specified in Fig. 16 at the last time step are shown in Fig. 17. Comparable results have been obtained using the two contact surface description methods. Along the transverse curves F, TM, and R, the normal contact stress distribution is U-shaped (small at the center and large at both sides), reflecting the stress concentration phenomenon. Along the longitudinal curves L and LM, the normal stress on the contact interface is large in the front and small in the middle and rear, indicating that the front part of the contactor’s bottom surface is subjected to a larger friction since friction is proportional to normal stress. Moreover, comparison between the results shows that the RPIM smoothing method produces stress distributions much smoother than the FE surface mesh description. 5.3. Convergence analysis To compare the convergence of both methods, i.e., the FE surface mesh and RPIM smoothing, the iterative process for the initial step (application of the normal pressure) of the contactor sliding problem is analyzed. A personal computer with R Intel⃝ CoreTM i5 CPU 760 @ 2.80 GHz, 8.0 GB RAM and ST3500418AS 7200 RPM Hard Drive is used to perform the required numerical analysis. Fig. 18(a) shows the variation of the residual during the iterative process, where residual refers to the L2 norm of the residual nodal force vector. With the RPIM smoothing method, the number of iterations is significantly reduced. Fig. 18(b) plots the number of nodes with different contact status at each iteration, from which it can be seen that the contact status of nodes stabilize faster when the RPIM smoothing method is used. The CPU time with the FE surface mesh is 12.38 s while the RPIM smoothing method takes a much shorter time, 5.086 s. However, when the contact surface is quite flat and the deformation is relatively small, the number of iterations using both methods may be comparable and the RPIM smoothing method may take a longer time due to the additional time required for contact surface smoothing. For example,
X. Qian et al. / Journal of Computational and Applied Mathematics 257 (2014) 1–13
11
Fig. 16. Curves for stress distribution plotting.
Fig. 17. Distributions of the normal contact stress along (a) the front curve F, (b) the transverse middle curve TM, (c) the rear curve R, (d) the left curve L, and (e) the longitudinal middle curve LM.
the number of iterations using both methods for the foundation and the subsoil example is the same, 9. The CPU times are 18.454 and 18.888 min for FE surface mesh and RPIM smoothing respectively. On the basis of its formulation, the RPIM smoothing method may also have the potential of overcoming the ‘chatter’ problem, often encountered when solving contact problems with FE surface meshes. A different FE mesh for the contactor sliding problem is shown in Fig. 19, where the mesh line Y = 2.5 m on the bottom of the contactor coincides with a mesh line on the inner surface of the half-tube. The variations of the residual and the number of nodes with different contact status at different iterations using both methods are displayed in Fig. 20. Due to the coincidence of mesh lines and the symmetry of this problem, slave nodes on the bottom line Y = 2.5 m on the contactor ‘jump’ between neighboring facets on the inner
12
X. Qian et al. / Journal of Computational and Applied Mathematics 257 (2014) 1–13
Fig. 18. Convergence analysis: (a) residual and (b) contact status.
Fig. 19. A different FE mesh for the contactor sliding problem.
Fig. 20. Convergence analysis (a different FE mesh): (a) residual and (b) contact status.
surface of the half-tube during the iterative process and hence impede the convergence when the FE surface mesh method is used. In contrast, it is clear from Fig. 20 that the RPIM smoothing method does not suffer from this problem. Although special treatments, e.g., allowing a slave node to contact with a master node instead of a contact facet, can be employed to deal with the particular problem encountered here, the proposed RPIM smoothing method provides a general way to avoid the ‘chatter’ problem. 6. Conclusions A general method for 3D contact surface smoothing has been proposed. This method is based on creating a smooth description of the local contact surface using the meshfree radial point interpolation method (RPIM). In particular, radial and polynomial basis functions are used to interpolate the scattered surface nodes and thus construct a smooth contact surface in the vicinity of any slave node likely to come into contact. With such a definition, this method applies to arbitrary surface meshes and is capable of reproducing smooth surfaces passing exactly through the surface nodes even for coarse meshes. The surface smoothing algorithm has been tested on different 3D surfaces. It shows a significant improvement in representing the real geometry of 3D bodies, especially for normal vector approximation. To demonstrate its advantages
X. Qian et al. / Journal of Computational and Applied Mathematics 257 (2014) 1–13
13
in contact analysis, numerical simulations of contact between deformable bodies have been performed, where adjustment of time step is employed to control the interpenetration. The results indicate that the proposed contact surface smoothing method can effectively reduce the oscillation of contact forces and non-smoothness of stress distributions. Through convergence analysis, the proposed RPIM smoothing method demonstrates a great potential of improving the convergence and overcoming the ‘chatter’ problem for contact analyses involving curved contact surfaces. Acknowledgments The authors would like to thank the financial support of the National Basic Research Program of China #2010CB732103, the National Natural Science Foundation of China #51209118, and the State Key Laboratory of Hydroscience and Engineering #2012-KY-02. References [1] D.J. Benson, J.O. Hallquist, A single surface contact algorithm for the post-buckling analysis of shell structures, Comput. Methods Appl. Mech. Eng. 91 (1990) 141–163. [2] J.O. Hallquist, G.L. Goudreau, D.J. Benson, Sliding interfaces with contact-impact in large-scale Lagrangian computations, Comput. Methods Appl. Mech. Eng. 51 (1985) 107–137. [3] Z.H. Zhong, L. Nilsson, A unified contact algorithm based on the territory concept, Comput. Methods Appl. Mech. Eng. 130 (1996) 1–16. [4] T. Belytschko, J.I. Lin, A three-dimensional impact-penetration algorithm with erosion, Comput. Struct. 25 (1987) 95–104. [5] T. Belytschko, M.O. Neal, Contact-impact by the pinball algorithm with penalty projection and Lagrangian methods, Internat. J. Numer. Methods Engrg. 31 (1991) 547–572. [6] S.P. Wang, E. Nakamachi, The inside-outside search algorithm for finite element analysis, Internat. J. Numer. Methods Engrg. 40 (1997) 3665–3685. [7] G. Zavarise, L. De Lorenzis, The node-to-segment algorithm for 2D frictionless contact: classical formulation and special cases, Comput. Methods Appl. Mech. Eng. 198 (2009) 3428–3451. [8] V. Padmanabhan, T.A. Laursen, A framework for development of surface smoothing procedures in large deformation frictional contact analysis, Finite Elem. Anal. Des. 37 (2001) 173–198. [9] P. Wriggers, L. Krstulovic-Opara, J. Korelc, Smooth C 1 -interpolations for two-dimensional frictional contact problems, Internat. J. Numer. Methods Engrg. 51 (2001) 1469–1495. [10] S. Ulaga, M. Ulbin, J. Flasker, Contact problems of gears using overhauser splines, Int. J. Mech. Sci. 41 (1999) 385–395. [11] N. El-Abbasi, S.A. Meguid, A. Czekanski, On the modelling of smooth contact surfaces using cubic splines, Internat. J. Numer. Methods Engrg. 50 (2001) 953–996. [12] M. Stadler, G.A. Holzapfel, J. Korelc, C n continuous modelling of smooth contact surfaces using NURBS and application to 2D problems, Internat. J. Numer. Methods Engrg. 57 (2003) 2177–2203. [13] G. Pietrzak, A. Curnier, Large deformation frictional contact mechanics: continuum formulation and augmented Lagrangian treatment, Comput. Methods Appl. Mech. Eng. 177 (1999) 351–381. [14] F.J. Wang, J.G. Cheng, Z.H. Yao, FFS contact searching algorithm for dynamic finite element analysis, Internat. J. Numer. Methods Engrg. 52 (2001) 655–672. [15] L. Krstulovic-Opara, P. Wriggers, J. Korelc, A C 1 -continuous formulation for 3D finite deformation frictional contact, Comput. Mech. 29 (2002) 27–42. [16] M.A. Puso, T.A. Laursen, A 3D contact smoothing method using Gregory patches, Internat. J. Numer. Methods Engrg. 54 (2002) 1161–1194. [17] D.M. Neto, M.C. Oliveira, L.F. Menezes, J.L. Alves, Improving Nagata patch interpolation applied for tool surface description in sheet metal forming simulation, Comput. Aided Des. 45 (2013) 639–656. [18] T. Hama, T. Nagata, C. Teodosiu, A. Makinouchi, H. Takuda, Finite-element simulation of springback in sheet metal forming using local interpolation for tool surfaces, Int. J. Mech. Sci. 50 (2008) 175–192. [19] D.M. Neto, M.C. Oliveira, L.F. Menezes, J.L. Alves, Nagata patch interpolation using surface normal vectors evaluated from the IGES file, Finite Elem. Anal. Des. 72 (2013) 35–46. [20] D. Chamoret, P. Saillard, A. Rassineux, J.M. Bergheau, New smoothing procedures in contact mechanics, J. Comput. Appl. Math. 168 (2004) 107–116. [21] T. Belytschko, W.J.T. Daniel, G. Ventura, A monolithic smoothing-gap algorithm for contact-impact based on the signed distance function, Internat. J. Numer. Methods Engrg. 55 (2002) 101–125. [22] J.G. Wang, G.R. Liu, A point interpolation meshless method based on radial basis functions, Internat. J. Numer. Methods Engrg. 54 (2002) 1623–1648. [23] T.A. Laursen, J.C. Simo, A continuum-based finite element formulation for the implicit solution of multibody large deformation frictional contact problems, Internat. J. Numer. Methods Engrg. 36 (1993) 3451–3485. [24] M.C. Wang, Finite Element Method, Tsinghua University Press, Beijing, 2003. [25] M.U. Rahman, R.E. Rowlands, R.D. Cook, T.L. Wilkinson, An iterative procedure for finite-element stress analysis of frictional contact problems, Comput. Struct. 18 (1984) 947–954. [26] A.B. Chaudhary, K.J. Bathe, A solution method for static and dynamic analysis of three-dimensional contact problems with friction, Comput. Struct. 24 (1986) 855–873. [27] I. Doghri, A. Muller, R.L. Taylor, A general three-dimensional contact procedure for implicit finite element codes, Eng. Comput. 15 (1998) 233–259. [28] M. Oldenburg, L. Nilsson, The position code algorithm for contact searching, Internat. J. Numer. Methods Engrg. 37 (1994) 359–386. [29] MSC. Software Corporation, Marc Volume A: Theory and User Information, MSC. Software Corporation, 2010, pp. 554–556. [30] M. de Berg, O. Cheong, M. van Kreveld, M. Overmars, Computational Geometry: Algorithms and Applications, Springer, Berlin, 2008.