- Email: [email protected]
United elements method for general contact–impact problems
Comput. Methods Appl. Mech. Engrg. 191 (2001) 843±860
www.elsevier.com/locate/cma
United elements method for general contact±impact problems K. Farahani *, M. Mo®d, A. Vafai Department of Civil Engineering, Sharif University of Technology, P.O. Box 11365-9313, Tehran, Iran Received 2 May 2000; received in revised form 4 September 2000
Abstract This work presents the details of a ®nite element-based numerical method to analyse general contact/impact problems, called united elements method (UEM). In contrast to usual methods, that is Lagrange Multipliers and Penalty Function, this method is conceptual, suitable for computer programming, and free from inherent de®ciencies of those methods. The key points of this procedure are elimination of the normal degree of freedom of the contactor node from the system energy potential functional in terms of target surface motion, and proper calculation of target element unbalanced forces, to incorporate the contact constraints. After a review on classic formulations, the mathematical and matrix presentation of the united element stiness matrix are discussed. Solution of some examples adopting a Newton±Raphson-based solution, together with a normal return scheme, support the capability and accuracy of the procedure in static and dynamic problems. Ó 2001 Published by Elsevier Science B.V.
1. Introduction Contact/impact problems are one of the most frequent problems of science and industry for design purposes or failure forces of structural systems. These problems are inherently nonlinear since the contact area is obtained as a part of the solution. During the last three decades, dierent numerical methods have been developed for the solution of these problems. Using these methods, contact/impact problems accompanying highly nonlinear factors due to large deformations, friction, and constitutive equations have been successfully treated. The core of these methods is how the contact constraints are included in the solution process. Most of the methods use variational formulations to incorporate the contact constraints between contacted bodies. The two constraints are: (1) equal displacement over the contacted areas, which is the impenetration condition, and (2) equal and reverse reactions, which is in fact Newtons's third law. In the variational methods, contact constraints are incorporated in the formulation from the beginning. Variational principles of contact problems are usually formulated in classical incremental forms using Lagrange Multipliers [1±7] or Penalty function methods [8±19], Perturbed Lagrangian method which is the combination of these two [20,21], gap elements [22±25] which is based on penalty formulation, and augmented Lagrangian method [26±28] which has its roots in constraint optimization theories. In this work, focusing on a node-to-surface contact situation which is the concern of most contact problems, and regarding the variational principles of motion and system potential energy functional, the steps of an incremental procedure is presented for contact/impact problems in detail. This method is based on physical understanding and concepts of the phenomena using ®nite element method. The primary concern of the presented procedure called united elements method (UEM) [29] is to cover some de®ciencies of the Lagrange Multipliers and Penalty-based formulations.
*
Corresponding author. E-mail address: [email protected] (K. Farahani).
0045-7825/01/$ - see front matter Ó 2001 Published by Elsevier Science B.V. PII: S 0 0 4 5 - 7 8 2 5 ( 0 1 ) 0 0 2 9 4 - 8
844
K. Farahani et al. / Comput. Methods Appl. Mech. Engrg. 191 (2001) 843±860
The four main steps of this method are to: (1) combine or unite the contactor and target elements by eliminating normal or another dependent d.o.f., (2) solve the incremental equations of motion and impose the contact condition between contactor node and target surface in subsequent iterations as long as their reaction is compression, (3) consider and exert the ``contactor node internal force'' as the ``target surface external force'', and (4) continue the iterations until the convergence criteria is met and the nonlinear FE equations converge. An interpretation of this technique is to treat the contact forces between two elements or two bodies, as the internal force of a super element or a united body which is the sum of the two, hence the name UEM. In this work, without loss of generality, the contact is assumed to be frictionless, the contactor node is massless, and the impact occurs without energy loss. Therefore, the only energy consuming source is possible plastic deformation. Also, tensor quantities are presented in bold characters and index notation is adopted. Regarding these assumptions, after a review on the usual formulations, the four fundamental steps of UEM procedure are explained together with some static and dynamic numerical example supports. Comparing or combination of this numerical technique with other methods is open for future and further researches. 2. Mathematical presentation of classical formulations Consider the system of two contacted bodies shown in Fig. 1. As widely used in contact literature, they are optionally called Contactor and Target bodies, with current volumes V C and V T and boundaries CC and CT , respectively, and common contact zone Cc . Also, UnC and UnT are the displacements of each body along their common opposite normals ~ n and ~ n. Before the commencement of contact, for each body we can write CXu [ CXt CX ;
CXu \ CXt ;;
2:1
where X stands for C or T, and ; is the null set. Right subscripts u and t show boundaries with prescribed displacements and tractions, respectively. However, after contact and developing interacting tractions s between the bodies, as long as their normal reaction is compression, we have CC \ CT Cc ;
V C \ V T ;:
2:2
The boundary condition compatibility or impenetration equation is UnC x; y; z UnT x; y; z;
x; y; z 2 Cc ;
with x; y; z as a space point on contact zone Cc .
Fig. 1. Contactor and target bodies: (a) continuum model; (b) ®nite element model.
2:3
K. Farahani et al. / Comput. Methods Appl. Mech. Engrg. 191 (2001) 843±860
845
Therefore, the constrained continuum mechanics problem to be solved is: i rij;j qbi qU
in V X ;
f r h F in V
X
2:4
;
rij ni tj in Ct X ; Ui Ui in C X u ; g P 0; g UnC UnT
2:5 2:6 2:7
G
in Cc ;
2:8
where rij ; bi ; tj ; Ui and ni are Cauchy stress, body force, traction, displacement vector, and outward normal vector components, respectively, with barred parameters as prescribed quantities. Also, F is deformation gradient tensor, and f( ) and h( ) are tensor time or space differential operators or functionals. Thus, (2.4) and (2.5) are equilibrium and constitutive equations, and (2.6) and (2.7) are traction and displacement boundary conditions, respectively. Also, (2.8) presents the contact constraint in which G and g are initial and current gap dimension functions, respectively. According to Fig. 1(a), applying the principle of virtual work for the bodies and ignoring the damping forces, the weak form of (2.4)±(2.8) may be written as Z Z Z Z Z C T Sij dEij dV qUi dUi dV sn d Un Un G dC qbi dUi dV ti dUi dC; 2:9 V
V
Cc
V
Ct
where S and E are two energetically conjugate stress and strain pairs with reference to the current con i is second time derivative of displacement vector Ui . Also, d means any in®nitesimal ®guration, and U variation which is consistent with displacement boundary conditions. In (2.9) it is assumed that Ct CCt [ CTt ;
V V C [ V T:
2:10
Eq. (2.9) with sn as unknown multiplier, present the Lagrange Multipliers approach for the contact problem solution. It explicitly satis®es impenetration condition (2.3) or (2.8), but results in additional unknowns sn and nonpositive-de®nite system of equations after ®nite element discretization. As a second approach, if the contact force or normal traction sn is assumed to be proportional to bodies penetration, using a penalty proportion parameter a, we have sn ag or sn a UnC
UnT G:
2:11
Hence, the no-overlap condition (2.8) is intentionally violated. Also numerical errors may arise because of a large parameter a and limited number of computer decimal digits. In this case the virtual work equation takes the form Z Z Z Z Z i dUi dV Sij dEij dV qU agdg dC qbi dUi dV ti dUi dC: 2:12 V
V
Cc
V
Ct
In both (2.9) and (2.12), the third term on left-hand side presents the virtual work of contact forces. It is also noted that dG 0. As a third approach called UEM, let us consider the two contacted bodies as a united body. Since the departure is not anticipated, the contact interactions forces are now as internal forces of a whole body, and the work of traction forces or the third term of left-hand side in (2.9) or (2.12) is unnecessary. Therefore, the virtual work equation takes the reduced form Z Z Z Z i dUi dV Sij dEij dV qU qbi dUi dV ti dUi dC; 2:13 V
V
V
Ct
which is the regular virtual work equation. After ®nite element discretization shown in Fig. 1(b), and using approximate displacement functions, (2.8) must be included in numerical calculation of the ®rst integral in left-hand side of (2.13). That is, depending on the applied elements in the contact zone of the bodies and assuming a node-to-surface contact, dUnC may be determined in terms of target surface
846
K. Farahani et al. / Comput. Methods Appl. Mech. Engrg. 191 (2001) 843±860
variational motion. Hence dUnC is omitted from the virtual work equation (2.13) or even the system potential energy functional.
3. Matrix presentation of classic formulations Assuming a ®nite element spatial discretization of the bodies, the static system potential functional before the contact takes place, may be considered as P 12U t KU
U tf ;
3:1
where K
KC 0
0 ; KT
U
UC ; UT
f
fC ; fT
where K, U, and f are stiffness matrix, displacement vector, and force vector of the whole system, respectively. Right superscript shows the body to which the quantity belongs, and t means transpose. Assuming in®nitesimal deformations, minimization of (3.1) relative to independent displacement variables will result in the system equilibrium equations. The stiffness matrix zero off-diagonals show that the two bodies have not yet contacted. After the contact occurs, the new constrained system potential functional in Lagrange Multipliers form may be written as PL 12U t KU
U t f gt k;
3:2
where g is nodal gap function vector, and k is independent Lagrange Multipliers vector. The ®rst-order Taylor's series expansion of g in terms of U may be given as g g0
og U: oU
3:3
Supposing that the constant square matrix g0 og=oU , (3.2) may be developed as PL 12U t KU
U t f g0 g0 U t k:
3:4
Hence, we can write oPL dU t KU
t
t
dU t f g0 dU k g0 g0 U dk:
Since for the equilibrium dPL 0, we arrive at t U f K g0 : g0 k g0 0
3:5
3:6
In this formulation, additional terms in the unknowns vector and diagonal zero terms in the enlarged coecient matrix are observed. In the Penalty form, system potential may be established as PP 12U t KU
U t f 12gt KP g;
3:7
where KP is the constant penalty matrix. Therefore oPP dU t KU
dU t f dgt KP g:
Again, making use of (3.3) and setting dPP 0 results in i t t K g0 KP g0 U f g0 KP g0 ;
3:8
3:9
K. Farahani et al. / Comput. Methods Appl. Mech. Engrg. 191 (2001) 843±860
847
which shows additional penalty terms in the stiness matrix and its positive de®niteness without enlargement.
4. UEM as an alternative algorithm In this section, the four steps of the contact alternative solution algorithm UEM is presented. The aim of the algorithm is to unite the two bodies so that their interacting forces turn into internal forces. For simplicity, an element-based presentation will be given and discussed ®rst. As mentioned before, these four steps are: (1) uniting the contactor and target elements or bodies by eliminating all of the dependent degrees, (2) solve the resulted equations and return the contactor node of the target surface, (3) impose the ``contactor node internal force'' as the ``target surface external force'', and (4) continue iterations until the convergence is reached. After this four steps are thoroughly outlined, the mathematical and matrix presentation of the method will be covered, respectively. 4.1. Step 1: uniting elements Since this procedure is element-based, it is easier to explain the procedure with a simple C0 continuous element. However, the method is general and can be extended to any other kind of elements, e.g. solids, beams, shells or any combination of them. As shown in Fig. 2, consider the Contactor and Target element with current or tangent stiness matrices K C and K T , respectively. The term current or tangent implies possible incorporation of any sort of nonlinearities in the problem. Ignoring the friction eect, it is assumed that the contact occurs between the contactor node 4 and target surface 2±3. As long as the two elements are in touch and their reaction is compression, it is simple to show that node 4 has only one independent d.o.f. as shown in Fig. 2(c). That is, if the target surface and tangent d.o.f. are restrained, the normal d.o.f. cannot move. Therefore, the normal degree can be eliminated from the set of dependent variables. A simple transformation is needed to rotate the contactor node d.o.f.'s of element C normal and parallel to the target surface as 0
K C TK C T t ;
4:1
where right superscript t means transpose and T is a rotation matrix.
Fig. 2. Contactor and target elements: (a) two contacted elements; (b) rotated d.o.f.'s; (c) dependency of normal d.o.f.; (d) multicontactor elements.
848
K. Farahani et al. / Comput. Methods Appl. Mech. Engrg. 191 (2001) 843±860
After the transformation (4.1), if we impose unit in®nitesimal deformation to u3 ; . . . ; u6 one by one, with other d.o.f.'s being restrained, according to the superposition principle and noting that in this case line 2±3 remains line, for in®nitesimal deformations we can write: u8 b sin au3 b cos au4 1
b sin au5 1
b cos au6 ;
4:2
where 0 6 b 6 1 and a are shown in Fig. 2(a). Eq. (4.2) may be called uniting or impenetration condition. 0 Therefore, having the stiness matrices K C and K T of contactor and target elements, we can write the incremental potential energy functional of the two elements as P 12U t KU n
U tP 0
T t
12 fU g fU C g
t
o K T 0
( fU T g ) n 0o 0 UC KC 0
( ) n o fP T g 0 t t n 0o ; fU T g fU C g PC
4:3:a
or P 12f u1
u2
u11
u12 g
KT 0
9 8 u1 > > > > > > > > > = < u2 > 0 .. 0 KC > > > . > > > > > > u11 > ; : u12
f u1
u2
u11
9 8 p1 > > > > > > > > > = < p2 > . .. : u12 g > > > > > > > > > p11 > ; : p12
4:3:b
Putting u8 from contact constraint (4.2) into (4.3.a) and (4.3.b), and forming the system of equations oP 0 oui
i 6 8
4:4
and rearranging the resulted equations, the united stiness matrix is obtained as
4:5 where C1 b sin a;
C2 b cos a;
C3 1
b sin a;
C4 1
b cos a
4:6
K. Farahani et al. / Comput. Methods Appl. Mech. Engrg. 191 (2001) 843±860
849
and Unite( ) is considered a matrix function with arguments K T ; K C , slope parameter a, and position parameter b. Since the stiness of normal freedom is accounted for in other degrees via additional terms in (4.5), it is eliminated from K U by omitting its corresponding row and column. However, this may disturb the sky line map, omitting of the row and column could be ignored by simple matrix manipulations [29]. The same united stiness matrix K U could be established from the de®nition of stiness matrix. That is, KijU is the reaction in d.o.f. i due to unit displacement in d.o.f. j, with other degrees being restrained [29]. In many cases, there may be several contactor elements Ci with stiness matrices K Ci versus one target element with K T , as shown in Fig. 2(d). In this case, the same procedure can be adopted to form a larger united element. It is noted that it is possible to eliminate other d.o.f.'s from Eq. (4.2) instead of u8 . Therefore, target surface degrees could be eliminated to unite the elements. Also, if the target surface degrees were rotated parallel to (4.1) shown in Fig. (3) as 0
K T T1 K T T1t
4:7
the uniting equation and consequently UE stiness matrix would get even simpler forms as u8 bu4 1
bu6
4:8
and
4:9
respectively. 4.2. Step 2: normal return After the K U is assembled in global stiness matrix and the equations are solved, since u8 is omitted from the equations, it must be computed in another way. For example u8 may consistently be computed from (4.2). However, since (4.2) is true only for in®nitesimal deformations, and in general ®nite movements occur, use of (4.2) causes the contactor node and target surface to either depart or penetrate. Therefore, an answer must be adopted for u8 which returns the contactor node to the target surface along the normal d.o.f. This step may be called normal return, similar to radial return in computational plasticity. When the node is returned back to the surface, contact boundary conditions are again met. However, since this forced
850
K. Farahani et al. / Comput. Methods Appl. Mech. Engrg. 191 (2001) 843±860
Fig. 3. Rotation of target surface degrees of freedom.
normal return will cause additional unbalanced forces, simultaneous satisfaction of boundary conditions and equilibrium is obtained via equilibrium iterations. 4.3. Step 3: contact force, from contactor node to target surface After obtaining incremental displacements and rotating them back into the global coordinates, the internal force vector is computed. As a usual nonlinear FEM scheme, the dierence between external and internal force vector of all the elements originates subsequent incremental deformations. However, the unbalanced force of target element must be computed from a dierent scheme. That is, the succeeding step is to externally load the target surface by normal internal force of contactor node (for a frictionless contact), then compute its unbalanced forces, in each equilibrium iteration. Hence, in addition to geometric and material nonlinearities, change and movement of contact force are other sources of target element unbalanced forces. 4.4. Step 4: convergence control The iterations must continue until the convergence criteria is reached. During the iterations, the contactor node may move from one target element to another. This is important since the contact forces must be imposed on proper target surface even if the UE matrix is not updated. It is noted that the time of contact and departure of elements must be determined with enough accuracy. Also, since the contact iterations may cause repeated spurious material loading and unloading in elasto-plastic analyses due to the normal return algorithm, special attention must be paid to the stress update scheme. For example, the reference stress state may be taken the one at the beginning of the time step. The ¯owchart given in Appendix A manifests a full Newton±Raphson scheme of the UEM algorithm. 4.5. Geometry approximation It is seen that Eqs. (4.1) and (4.2) of step 1 have most importance in constitution of K U and they completely depend on the geometry and formulation of the applied elements. For higher-order elements or target surfaces, these relations will take more complicated forms. However, for such cases some geometric approximations may be adopted to simplify (4.1) and (4.2) so as to speed up the calculations with any desired accuracy. These approximations may be compensated in other steps and equilibrium iterations. As an example, Fig. 4 shows an eight node isoparametric element with a secondorder target surface/curve which is approximated with two lines. This will highly reduce the form of impenetration equation and the number of additional terms in the united element stiness matrix. Hence, instead of writing K U Unite K T ; K C ; a; b, depending on the contact position we may write
K. Farahani et al. / Comput. Methods Appl. Mech. Engrg. 191 (2001) 843±860
851
Fig. 4. Geometry approximation: (a) exact geometry; (b) approximate geometry.
Fig. 5. Higher-order target surface: (a) contact position; (b) degrees of freedom.
K U Unite K T ; K C ; ai ; bi where i depends on the contact zone shown in Fig. 4(b). This approximation which offers ¯exibility to the procedure to utilize higher-order elements is used only in step 1, while other steps may continue with exact element geometry. To show the effect of higher-order target surface on the complexity uniting equation, consider the parabolic target surface shown in Fig. (5). It is easy to show that uniting equation for the rotated degrees can be written as 6 X oY oY oY oY oY oY u12 cos a u1 u2 u3 u4 u5 u6 M i ui ; 4:10 oX1 oY1 oX2 oY2 oX3 oY3 i1 X X6
where Xi and Yi are coordinates of node number i. Also, Y is the equation of target curve in XY plain as Y aX 2 bX c where 2 3 2 2 a X1 4b5 4X2 2 c X32
X1 X2 X3
3 12 3 Y1 1 1 5 4 Y2 5 : Y3 1
4:11
As an example, M1 in Eq. (4.10) is
3
M1 A X6 where A
X1
X2 X6 X2 X2
X3 X2 X3 X3
Y2 X3 X1
Y1 Y3 X2 X1
Y1 Y2 X3 X2
X1 and X6 X3 b X1
X3 :
Y3 ; X3
4:12
852
K. Farahani et al. / Comput. Methods Appl. Mech. Engrg. 191 (2001) 843±860
Therefore, in such cases, a large amount of computation time must be spent to extract the uniting equation coecients.
5. Mathematical presentation of UEM Following the procedure of Section 4, the mathematical form of uniting equation in a three-dimensional space may be derived to unite any arbitrary contacted bodies with any number of contactor nodes. Assume that the potential energy function of the bodies shown in Fig. 1(b) before contact is P de®ned as P P u1 ; . . . ; un ;
5:1
where n is the total number of freedom degrees of all the bodies. Therefore, the stiffness matrix of the system can be de®ned as Kij
o2 P : oui ouj
5:2
After contact, a series of uniting equations or contact constraints such as uk F u1 ; . . . ; ui ; . . . ; un ;
i 6 k;
5:3
may hold, where uk 's are the normal d.o.f. of contactor nodes obtained after rotation. Also, all the arguments of contact constraint function F( ) correspond to d.o.f.'s of the target surfaces. Substituting (5.3) for the system as into (5.1) yields as new potential P P u1 ; . . . ; un P ; 5:4:a uk F u1 ;...;ui ;...;un or P u 1 ; . . . ; ui ; . . . ; un ; P
i 6 k:
5:4:b
Therefore, the new tangent united stiness matrix of the structure, K U , is de®ned as KijU
o2 P ; oui ouj
i; j 6 k:
5:5
Taking the partial dierentiation relative to i and j and one by one yields o oP o oP oP oF U Kij ; ouj oui ouj oui oF oui uk F u1 ;...;ui ;...;un
5:6
or KijU
o2 P o2 P oF oui ouj oui oF ouj
o2 P o2 P oF 2 ouj oF oF ouj
oF oP o2 F : oui oF oui ouj
5:7
Therefore, linearizing the uniting function F will both eliminate the last term of right-hand side of (5.5) and simplify its differentiation. After rearrangement and regarding (5.2) and (5.7) may be rewritten as 2 2 o P oF o2 P oF o P oF oF oP o2 F KijU Kij ; 5:8 oui oF ouj ouj oF oui oF 2 oui ouj oF oui ouj which thoroughly exhibits the quality and origin of additional terms in the united stiness matrix. It is noted that equations such as (2.3), (2.8), and (5.3) may alternatively be called the no-overlap, impenetration, contact, or uniting constraint.
K. Farahani et al. / Comput. Methods Appl. Mech. Engrg. 191 (2001) 843±860
853
6. Matrix presentation of UEM Again, after ®nite element discretization of the bodies, the static system potential functional before the contact is as given in (3.1). However, after the contact, the contactor nodes normal displacement vector UnC , can be described in terms of target surfaces deformation vector U T as UnC F U T ;
6:1
which again indicates the uniting equation. In general, (6.1) can be written as UnC AU T ;
6:2
where A is the connection matrix. The higher is the order of the target surface, the more complicated are the components of A, as shown in Section 4.5. The potential function after contact may be presented as P 12U t KU
U tf ;
6:3
where
6:4
Eq. (6.4) is simply a partitioning of stiness matrix and force and displacement vectors, into contactor normal, target surface, and the rest of freedom degrees, corresponding to C, T, and R right superscripts, respectively. Zero o-diagonals show that there is no physical relation between them before contact. After substitution of (6.1) into (6.3), we can write:
6:5 Again, noting the matrix oF =oU T A; considering dP 0 and rearranging the equations result in t cc T T U A K A K TT At K CR K TR f At fnC : 6:6 K CR A K TR K RR UR fR The stiness matrix in (6.6) is in fact what we called UE matrix or K U in the level of structure. It is noted that here, the number of rows and columns of connection matrix A is equal to the number of contactor nodes and the number of target surface d.o.f.'s, respectively. For example, the matrices A and K CC , corresponding to contact equation (4.2) are A hb sin a b cos a 1 h 0i C K CC K22 :
b sin a 1
b cos ai;
6:7 6:8
Although most of the discussions above concerned the stiness part of the contacting structures, they can be generalized to the problems including friction, damping, and inertial forces. After the solution of linear system of (6.6), UnC can be obtained from (6.2) for in®nitesimal deformations. However, for ®nite deformations, the four steps of Section 4 must be used to arrive in simultaneous equilibrium and boundary conditions satisfaction. Therefore, for in®nitesimal deformations, Eqs. (6.6) and (6.2) present a positive de®nite coecient matrix and exact satisfaction of boundary conditions, respectively, with no stiness illconditioning or addition in number of unknowns.
854
K. Farahani et al. / Comput. Methods Appl. Mech. Engrg. 191 (2001) 843±860
7. Numerical examples Based on the proposed procedure, a computer program is developed and static and dynamic numerical examples are given below. In the following examples, modi®ed Newton±Raphson iterative solution together with normal return scheme are used to solve the nonlinear equilibrium equations. Material is assumed to be elastic or elasto-plastic and rate-independent and the eect of friction is ignored. 7.1. Example 1 Fig. 6 shows two elastic beams each constituted of eight plain stress four node elements. The mechanical and geometrical properties are given in the ®gure. The total Lagrangian incremental formulation is chosen. The adopted material constitutive equation for ®nite elastic strains is Sij 2lEij kdij Ekk ;
7:1
where S and E are second Piola±Kirchhoff stress and Green±Lagrange strain tensors, respectively. k and l are material constants, d is the Kronecker delta, and index notation is adopted [28]. This contact problem shows a snap-through behavior. It is noted that if the strains become ®nite, the eect of using (7.1) on material response in ®nite elastic strains is important and is beyond the scope of this work. However. the use of dierent constitutive rate models which are based on objective rates of Cauchy stress (real stress) is well established. Fig. 7 shows three important phases of deformation, that is, contact, snap through, and departure. For the snap through part, a displacement control scheme is utilized to capture the history of motion. The contact is assumed only for the contactor node A shown in Fig. 6. Curves 1 and 2 of Fig. 8 show load±de¯ection curve of the point under load and contactor nodes, respectively. Curve 3 of the ®gure is the history of contact force. Points A, B, and C of curve 2 are relevant to three phases of Fig. 7. 7.2. Example 2 Fig. 9 shows two bodies again each constituted of 8 plain stress four node elements, with given geometric, mechanical, and initial conditions. Again the constitutive equation (7.1) is used to model the elastic
Fig. 6. Snap-through contact.
Fig. 7. Deformation phases.
K. Farahani et al. / Comput. Methods Appl. Mech. Engrg. 191 (2001) 843±860
855
Fig. 8. Load±de¯ection curves at dierent nodes.
Fig. 9. Dynamic contact.
part of material behavior. The lower body is massless and the upper one has concentrated masses Mx My 1 at each node except for the lower nodes. Gravity eect is not included. The chosen dynamic time integration scheme is Newmark with parameters d 0:5; a 0:25, and time step Dt 0:01 s. The upper body has an elastic±perfect plastic behavior with von Misez yield stress Fy 800, while the lower body is elastic. Fig. 10 shows the impact force history between the bodies, and Fig. 11 shows the deformation history of the impact. Plastic deformation in contactor element is apparent. 7.3. Example 3 Fig. 12 shows two bodies each constituted of plain stress four node elements, with given geometric and mechanical conditions. The hatched elements are much stier than the other two elements and the upper body is moving toward the lower one with initial velocity V0 10 and gravity acceleration G 10. Again the constitutive equation (7.1) is used for elastic material behavior. The upper body is elastic±perfect plastic with yield stress Fy 7:5. Newmark time integration with parameters d 0:5; a 0:25 and time step Dt 0:01 s is used. The height parameter h shown in the ®gure must be large enough to reduce the con®nement effect on the two soft elements, and in this example h 60. Due to relative stiffness of the elements, the problem is equivalent to the one shown in Fig. 13. Exact solution of the elastic±perfect plastic impact problem shown in Fig. 13 can be presented in six time intervals as
856
K. Farahani et al. / Comput. Methods Appl. Mech. Engrg. 191 (2001) 843±860
Fig. 10. Elastic and elasto-plastic impact force history.
Fig. 11. Deformation history.
Fig. 12. Elasto-plastic impact.
K. Farahani et al. / Comput. Methods Appl. Mech. Engrg. 191 (2001) 843±860
Fig. 13. Equivalent problem.
1 2 3
u 10; u 0 0; u_ 0:31 P t P 0 ) u 5t2
10; 10t:
u1 u 3:5804; t1 t 0:31; p p sin 200t1 0:05 cos 200t1 0:0643 P t1 P 0 ) u1 p13:1 200 u2 u 4:3305; t2 t 0:3743; 0:0615 P t2 P 0 ) u2 70t22 5:609t2 :
Fig. 14. Comparison of displacement history.
Fig. 15. Comparison of impact force history.
0:05:
857
858
K. Farahani et al. / Comput. Methods Appl. Mech. Engrg. 191 (2001) 843±860
Fig. A.1. Full Newton±Raphson iterative scheme with UEM and normal return.
4
u3 u 3:8452; t3 t 0:4358; p 0:1161 P t3 P 0 ) u3 0:7 cos 200t3
0:05:
K. Farahani et al. / Comput. Methods Appl. Mech. Engrg. 191 (2001) 843±860
5 6
859
U4 u 3:8452; t4 t 0:5519; 1:9748 P t4 P 0 ) u4 5t42 9:874t4 : U5 u 3:8452; t5 t 2:5267; p p sin 200t5 0:05 cos 200t5 0:2323 P t5 P 0 ) u5 p9:874 200
0:05:
Figs. 14 and 15 are the comparison between exact and numerical deformation and impact force history of the problem after two subsequent contact and departure, respectively. In both of the above ®gures, good agreement is observed between the exact solution and the numerical results. 8. Conclusions Regarding the incremental variational principles of contact problem, or based on the incremental form of the potential energy functional, the detailed mathematical and numerical procedure of UEM is presented in detail. The method focuses on node-to-surface contact/impacts which is the concern of most contact/ impact problems, and node-to-surface contact/impacts which is the concern of most contact/impact problems, and imposes the constraint of impenetration accurately without the drawbacks of Lagrange Multipliers and Penalty Function methods. That is, the number of unknowns and solution scheme do not change, and the boundary conditions are satis®ed exactly without loss of accuracy. The main idea is to omit normal d.o.f. of the contactors node by eliminating it from the system potential functional via uniting equations, and to use a normal return scheme to return the contactor node back to the target surface. The procedure can also be interpreted as treatment of the contact force between two bodies as the internal force of a whole body and can also be extended to frictional problems. The elemental, mathematical, and matrix presentation of the uniting step are also presented. Illustrative examples are given to demonstrate the capability and accuracy of the method. Enhancing or comparing this technique with other methods may be the matter of future works. Appendix A For a general nonlinear FEM-based contact/impact computer code, the four steps of UEM combined with a full Newton±Raphson iterative scheme may be illustrated in the ¯owchart of Fig. A.1. References [1] T.J.R. Hughes, R.L. Taylor, J.L. Sackman, A. Curnier, W. Kanoknukulchai, A ®nite element method for a class of contact± impact problems, Computer Methods in Applied Mechanics and Engineering 8 (1976) 249±276. [2] K.J. Bathe, A.B. Chaudhary, A solution method for static and dynamic analysis of three-dimensional contact problems with friction, Computers and Structures 24 (1986) 855±873. [3] N. Asano, An approximate hybrid type of virtual work principle for two elastoimpact contact bodies, Bulletin of JSME 26 (1983) 1849±1856. [4] N. Asano, A hybrid type of virtual work principle for impact contact problems of two bodies, Bulletin of JSME 29 (1986) 1679± 1684. [5] N.J. Carpenter, R.L. Taylor, M.G. Katona, Lagrange constraints for transient ®nite element surface contact, International Journal of Numerical Methods in Engineering 32 (1991) 103±128. [6] Z.-H. Zhong, L.G. Nilsson, Lagrange multiplier approach for evaluation of friction in explicit ®nite-element analysis, International Journal of Numerical Methods in Engineering 10 (1994) 249±255. [7] J.C. Simo, P. Wriggers, K.H. Schweizerhof, R.L. Taylor, Finite deformation post-buckling analysis involving inelasticity and contact constraints, International Journal of Numerical Methods in Engineering 23 (1986) 779±800. [8] N. Asano, A virtual work principle using penalty function method for impact contact problems of two bodies, Bulletin of JSME 29 (1986) 731±736. [9] N. Asano, A penalty function type of virtual work principle for impact contact problems of two bodies, Bulletin of JSME 29 (1986) 3701±3709. [10] Y. Kanto, G. Yagawa, A dynamic contact buckling analysis by the penalty ®nite element method, International Journal of Numerical Methods in Engineering 29 (1990) 755±774.
860
K. Farahani et al. / Comput. Methods Appl. Mech. Engrg. 191 (2001) 843±860
[11] R.F. Kulak, Adaptive contact elements for three-dimensional explicit transient analysis, Computer Methods in Applied Mechanics and Engineering 72 (1989) 125±151. [12] A.B. Pifko, R. Winter, Theory and application of ®nite element analysis to structural crash simulation, Computers and Structures 13 (1981) 277±285. [13] T. Belytschko, M.O. Neal, Contact±impact by pinball algorithm with penalty and Lagrangian methods, International Journal of Numerical Methods in Engineering 13 (1991) 547±572. [14] H.M. De La Fuente, C.A. Fellipa, Ephemeral penalty functions for contact±impact dynamics, Finite Elements in Analysis and Design 9 (1991) 177±191. [15] I. Hunek, On a penalty formulation for contact±impact problems, Computers and Structures 11 (1993) 193±203. [16] N. Kikuchi, A smoothing technique for reduce integration penalty methods in contact problems, International Journal Numerical methods Engineering 18 (1982) 343±350. [17] T. Campos, J.T. Oden, N. Kikuchi, A numerical analysis of a class of contact problems with friction in elastostatics, Computer Methods in Applied Mechanics and Engineering 34 (1982) 821±845. [18] B. Nour-Omid, P. Wriggers, A two level iteration method for solution of contact problems, Computer Methods in Applied Mechanics and engineering 54 (1986) 131±144. [19] J.C. Simo, P. Wriggers, R.L. Taylor, A perturbed Lagrangian formulation for the ®nite element solution of contact problems, Computer Methods in Applied Mechanics and Engineering 50 (1985) 163±180. [20] P. Wriggers, W. Wagner, E. Stein, Algorithms for nonlinear contact constraints with application to stability problems of rods and shells, Computational Mechanics 2 (1987) 215±230. [21] S.H. Lee, Rudimentary considerations for adaptive gap/friction element based on the penalty method, Computers and Structures 47 (1993) 1043±1056. [22] J. Ghaboussi, E.L. Wilson, J. Isenberg, Finite element of rock joints and interfaces, Journal of Soil Mechanics and Foundations Division, ASCE 99 (1973) 833±848. [23] S. Desai, M.M. Zaman, J.G. Lightner, H.J. Siriwardane, Thin-layer element for interfaces and joints, International Journal of Numerical and Analytical Methods in Geomechanics 8 (1984) 19±43. [24] P. Kowalczyk, Finite±deformation interface formulation for frictionless contact problems, Communications in Numerical Methods in Engineering 10 (1994) 879±893. [25] J.C. Simo, T.A. Laursen, An augmented Lagrangian treatment of contact problems involving friction, Computers and Structures 42 (1992) 97±116. [26] J.-H. Heegaard, A. Curnier, An augmented Lagrangian method for discrete large-slip contact problems involving friction, International Journal of Numerical Methods in Engineering 36 (1993) 569±593. [27] T.A. Laursen, S. Govindjee, A note on the treatment of frictionless contact between non-smooth surfaces in fully-nonlinear problems, Communications in Numerical Methods in Engineering 10 (1994) 869±878. [28] K.J. Bathe, Finite Element Procedures, Prentice-Hall, Englewood Clis, NJ, 1996. [29] K. Farahani, M. Mo®d, A. Vafai, A solution method for general contact-impact problems, Computer Methods in Applied Mechanics and Engineering 187 (2000) 69±77.
www.elsevier.com/locate/cma
United elements method for general contact±impact problems K. Farahani *, M. Mo®d, A. Vafai Department of Civil Engineering, Sharif University of Technology, P.O. Box 11365-9313, Tehran, Iran Received 2 May 2000; received in revised form 4 September 2000
Abstract This work presents the details of a ®nite element-based numerical method to analyse general contact/impact problems, called united elements method (UEM). In contrast to usual methods, that is Lagrange Multipliers and Penalty Function, this method is conceptual, suitable for computer programming, and free from inherent de®ciencies of those methods. The key points of this procedure are elimination of the normal degree of freedom of the contactor node from the system energy potential functional in terms of target surface motion, and proper calculation of target element unbalanced forces, to incorporate the contact constraints. After a review on classic formulations, the mathematical and matrix presentation of the united element stiness matrix are discussed. Solution of some examples adopting a Newton±Raphson-based solution, together with a normal return scheme, support the capability and accuracy of the procedure in static and dynamic problems. Ó 2001 Published by Elsevier Science B.V.
1. Introduction Contact/impact problems are one of the most frequent problems of science and industry for design purposes or failure forces of structural systems. These problems are inherently nonlinear since the contact area is obtained as a part of the solution. During the last three decades, dierent numerical methods have been developed for the solution of these problems. Using these methods, contact/impact problems accompanying highly nonlinear factors due to large deformations, friction, and constitutive equations have been successfully treated. The core of these methods is how the contact constraints are included in the solution process. Most of the methods use variational formulations to incorporate the contact constraints between contacted bodies. The two constraints are: (1) equal displacement over the contacted areas, which is the impenetration condition, and (2) equal and reverse reactions, which is in fact Newtons's third law. In the variational methods, contact constraints are incorporated in the formulation from the beginning. Variational principles of contact problems are usually formulated in classical incremental forms using Lagrange Multipliers [1±7] or Penalty function methods [8±19], Perturbed Lagrangian method which is the combination of these two [20,21], gap elements [22±25] which is based on penalty formulation, and augmented Lagrangian method [26±28] which has its roots in constraint optimization theories. In this work, focusing on a node-to-surface contact situation which is the concern of most contact problems, and regarding the variational principles of motion and system potential energy functional, the steps of an incremental procedure is presented for contact/impact problems in detail. This method is based on physical understanding and concepts of the phenomena using ®nite element method. The primary concern of the presented procedure called united elements method (UEM) [29] is to cover some de®ciencies of the Lagrange Multipliers and Penalty-based formulations.
*
Corresponding author. E-mail address: [email protected] (K. Farahani).
0045-7825/01/$ - see front matter Ó 2001 Published by Elsevier Science B.V. PII: S 0 0 4 5 - 7 8 2 5 ( 0 1 ) 0 0 2 9 4 - 8
844
K. Farahani et al. / Comput. Methods Appl. Mech. Engrg. 191 (2001) 843±860
The four main steps of this method are to: (1) combine or unite the contactor and target elements by eliminating normal or another dependent d.o.f., (2) solve the incremental equations of motion and impose the contact condition between contactor node and target surface in subsequent iterations as long as their reaction is compression, (3) consider and exert the ``contactor node internal force'' as the ``target surface external force'', and (4) continue the iterations until the convergence criteria is met and the nonlinear FE equations converge. An interpretation of this technique is to treat the contact forces between two elements or two bodies, as the internal force of a super element or a united body which is the sum of the two, hence the name UEM. In this work, without loss of generality, the contact is assumed to be frictionless, the contactor node is massless, and the impact occurs without energy loss. Therefore, the only energy consuming source is possible plastic deformation. Also, tensor quantities are presented in bold characters and index notation is adopted. Regarding these assumptions, after a review on the usual formulations, the four fundamental steps of UEM procedure are explained together with some static and dynamic numerical example supports. Comparing or combination of this numerical technique with other methods is open for future and further researches. 2. Mathematical presentation of classical formulations Consider the system of two contacted bodies shown in Fig. 1. As widely used in contact literature, they are optionally called Contactor and Target bodies, with current volumes V C and V T and boundaries CC and CT , respectively, and common contact zone Cc . Also, UnC and UnT are the displacements of each body along their common opposite normals ~ n and ~ n. Before the commencement of contact, for each body we can write CXu [ CXt CX ;
CXu \ CXt ;;
2:1
where X stands for C or T, and ; is the null set. Right subscripts u and t show boundaries with prescribed displacements and tractions, respectively. However, after contact and developing interacting tractions s between the bodies, as long as their normal reaction is compression, we have CC \ CT Cc ;
V C \ V T ;:
2:2
The boundary condition compatibility or impenetration equation is UnC x; y; z UnT x; y; z;
x; y; z 2 Cc ;
with x; y; z as a space point on contact zone Cc .
Fig. 1. Contactor and target bodies: (a) continuum model; (b) ®nite element model.
2:3
K. Farahani et al. / Comput. Methods Appl. Mech. Engrg. 191 (2001) 843±860
845
Therefore, the constrained continuum mechanics problem to be solved is: i rij;j qbi qU
in V X ;
f r h F in V
X
2:4
;
rij ni tj in Ct X ; Ui Ui in C X u ; g P 0; g UnC UnT
2:5 2:6 2:7
G
in Cc ;
2:8
where rij ; bi ; tj ; Ui and ni are Cauchy stress, body force, traction, displacement vector, and outward normal vector components, respectively, with barred parameters as prescribed quantities. Also, F is deformation gradient tensor, and f( ) and h( ) are tensor time or space differential operators or functionals. Thus, (2.4) and (2.5) are equilibrium and constitutive equations, and (2.6) and (2.7) are traction and displacement boundary conditions, respectively. Also, (2.8) presents the contact constraint in which G and g are initial and current gap dimension functions, respectively. According to Fig. 1(a), applying the principle of virtual work for the bodies and ignoring the damping forces, the weak form of (2.4)±(2.8) may be written as Z Z Z Z Z C T Sij dEij dV qUi dUi dV sn d Un Un G dC qbi dUi dV ti dUi dC; 2:9 V
V
Cc
V
Ct
where S and E are two energetically conjugate stress and strain pairs with reference to the current con i is second time derivative of displacement vector Ui . Also, d means any in®nitesimal ®guration, and U variation which is consistent with displacement boundary conditions. In (2.9) it is assumed that Ct CCt [ CTt ;
V V C [ V T:
2:10
Eq. (2.9) with sn as unknown multiplier, present the Lagrange Multipliers approach for the contact problem solution. It explicitly satis®es impenetration condition (2.3) or (2.8), but results in additional unknowns sn and nonpositive-de®nite system of equations after ®nite element discretization. As a second approach, if the contact force or normal traction sn is assumed to be proportional to bodies penetration, using a penalty proportion parameter a, we have sn ag or sn a UnC
UnT G:
2:11
Hence, the no-overlap condition (2.8) is intentionally violated. Also numerical errors may arise because of a large parameter a and limited number of computer decimal digits. In this case the virtual work equation takes the form Z Z Z Z Z i dUi dV Sij dEij dV qU agdg dC qbi dUi dV ti dUi dC: 2:12 V
V
Cc
V
Ct
In both (2.9) and (2.12), the third term on left-hand side presents the virtual work of contact forces. It is also noted that dG 0. As a third approach called UEM, let us consider the two contacted bodies as a united body. Since the departure is not anticipated, the contact interactions forces are now as internal forces of a whole body, and the work of traction forces or the third term of left-hand side in (2.9) or (2.12) is unnecessary. Therefore, the virtual work equation takes the reduced form Z Z Z Z i dUi dV Sij dEij dV qU qbi dUi dV ti dUi dC; 2:13 V
V
V
Ct
which is the regular virtual work equation. After ®nite element discretization shown in Fig. 1(b), and using approximate displacement functions, (2.8) must be included in numerical calculation of the ®rst integral in left-hand side of (2.13). That is, depending on the applied elements in the contact zone of the bodies and assuming a node-to-surface contact, dUnC may be determined in terms of target surface
846
K. Farahani et al. / Comput. Methods Appl. Mech. Engrg. 191 (2001) 843±860
variational motion. Hence dUnC is omitted from the virtual work equation (2.13) or even the system potential energy functional.
3. Matrix presentation of classic formulations Assuming a ®nite element spatial discretization of the bodies, the static system potential functional before the contact takes place, may be considered as P 12U t KU
U tf ;
3:1
where K
KC 0
0 ; KT
U
UC ; UT
f
fC ; fT
where K, U, and f are stiffness matrix, displacement vector, and force vector of the whole system, respectively. Right superscript shows the body to which the quantity belongs, and t means transpose. Assuming in®nitesimal deformations, minimization of (3.1) relative to independent displacement variables will result in the system equilibrium equations. The stiffness matrix zero off-diagonals show that the two bodies have not yet contacted. After the contact occurs, the new constrained system potential functional in Lagrange Multipliers form may be written as PL 12U t KU
U t f gt k;
3:2
where g is nodal gap function vector, and k is independent Lagrange Multipliers vector. The ®rst-order Taylor's series expansion of g in terms of U may be given as g g0
og U: oU
3:3
Supposing that the constant square matrix g0 og=oU , (3.2) may be developed as PL 12U t KU
U t f g0 g0 U t k:
3:4
Hence, we can write oPL dU t KU
t
t
dU t f g0 dU k g0 g0 U dk:
Since for the equilibrium dPL 0, we arrive at t U f K g0 : g0 k g0 0
3:5
3:6
In this formulation, additional terms in the unknowns vector and diagonal zero terms in the enlarged coecient matrix are observed. In the Penalty form, system potential may be established as PP 12U t KU
U t f 12gt KP g;
3:7
where KP is the constant penalty matrix. Therefore oPP dU t KU
dU t f dgt KP g:
Again, making use of (3.3) and setting dPP 0 results in i t t K g0 KP g0 U f g0 KP g0 ;
3:8
3:9
K. Farahani et al. / Comput. Methods Appl. Mech. Engrg. 191 (2001) 843±860
847
which shows additional penalty terms in the stiness matrix and its positive de®niteness without enlargement.
4. UEM as an alternative algorithm In this section, the four steps of the contact alternative solution algorithm UEM is presented. The aim of the algorithm is to unite the two bodies so that their interacting forces turn into internal forces. For simplicity, an element-based presentation will be given and discussed ®rst. As mentioned before, these four steps are: (1) uniting the contactor and target elements or bodies by eliminating all of the dependent degrees, (2) solve the resulted equations and return the contactor node of the target surface, (3) impose the ``contactor node internal force'' as the ``target surface external force'', and (4) continue iterations until the convergence is reached. After this four steps are thoroughly outlined, the mathematical and matrix presentation of the method will be covered, respectively. 4.1. Step 1: uniting elements Since this procedure is element-based, it is easier to explain the procedure with a simple C0 continuous element. However, the method is general and can be extended to any other kind of elements, e.g. solids, beams, shells or any combination of them. As shown in Fig. 2, consider the Contactor and Target element with current or tangent stiness matrices K C and K T , respectively. The term current or tangent implies possible incorporation of any sort of nonlinearities in the problem. Ignoring the friction eect, it is assumed that the contact occurs between the contactor node 4 and target surface 2±3. As long as the two elements are in touch and their reaction is compression, it is simple to show that node 4 has only one independent d.o.f. as shown in Fig. 2(c). That is, if the target surface and tangent d.o.f. are restrained, the normal d.o.f. cannot move. Therefore, the normal degree can be eliminated from the set of dependent variables. A simple transformation is needed to rotate the contactor node d.o.f.'s of element C normal and parallel to the target surface as 0
K C TK C T t ;
4:1
where right superscript t means transpose and T is a rotation matrix.
Fig. 2. Contactor and target elements: (a) two contacted elements; (b) rotated d.o.f.'s; (c) dependency of normal d.o.f.; (d) multicontactor elements.
848
K. Farahani et al. / Comput. Methods Appl. Mech. Engrg. 191 (2001) 843±860
After the transformation (4.1), if we impose unit in®nitesimal deformation to u3 ; . . . ; u6 one by one, with other d.o.f.'s being restrained, according to the superposition principle and noting that in this case line 2±3 remains line, for in®nitesimal deformations we can write: u8 b sin au3 b cos au4 1
b sin au5 1
b cos au6 ;
4:2
where 0 6 b 6 1 and a are shown in Fig. 2(a). Eq. (4.2) may be called uniting or impenetration condition. 0 Therefore, having the stiness matrices K C and K T of contactor and target elements, we can write the incremental potential energy functional of the two elements as P 12U t KU n
U tP 0
T t
12 fU g fU C g
t
o K T 0
( fU T g ) n 0o 0 UC KC 0
( ) n o fP T g 0 t t n 0o ; fU T g fU C g PC
4:3:a
or P 12f u1
u2
u11
u12 g
KT 0
9 8 u1 > > > > > > > > > = < u2 > 0 .. 0 KC > > > . > > > > > > u11 > ; : u12
f u1
u2
u11
9 8 p1 > > > > > > > > > = < p2 > . .. : u12 g > > > > > > > > > p11 > ; : p12
4:3:b
Putting u8 from contact constraint (4.2) into (4.3.a) and (4.3.b), and forming the system of equations oP 0 oui
i 6 8
4:4
and rearranging the resulted equations, the united stiness matrix is obtained as
4:5 where C1 b sin a;
C2 b cos a;
C3 1
b sin a;
C4 1
b cos a
4:6
K. Farahani et al. / Comput. Methods Appl. Mech. Engrg. 191 (2001) 843±860
849
and Unite( ) is considered a matrix function with arguments K T ; K C , slope parameter a, and position parameter b. Since the stiness of normal freedom is accounted for in other degrees via additional terms in (4.5), it is eliminated from K U by omitting its corresponding row and column. However, this may disturb the sky line map, omitting of the row and column could be ignored by simple matrix manipulations [29]. The same united stiness matrix K U could be established from the de®nition of stiness matrix. That is, KijU is the reaction in d.o.f. i due to unit displacement in d.o.f. j, with other degrees being restrained [29]. In many cases, there may be several contactor elements Ci with stiness matrices K Ci versus one target element with K T , as shown in Fig. 2(d). In this case, the same procedure can be adopted to form a larger united element. It is noted that it is possible to eliminate other d.o.f.'s from Eq. (4.2) instead of u8 . Therefore, target surface degrees could be eliminated to unite the elements. Also, if the target surface degrees were rotated parallel to (4.1) shown in Fig. (3) as 0
K T T1 K T T1t
4:7
the uniting equation and consequently UE stiness matrix would get even simpler forms as u8 bu4 1
bu6
4:8
and
4:9
respectively. 4.2. Step 2: normal return After the K U is assembled in global stiness matrix and the equations are solved, since u8 is omitted from the equations, it must be computed in another way. For example u8 may consistently be computed from (4.2). However, since (4.2) is true only for in®nitesimal deformations, and in general ®nite movements occur, use of (4.2) causes the contactor node and target surface to either depart or penetrate. Therefore, an answer must be adopted for u8 which returns the contactor node to the target surface along the normal d.o.f. This step may be called normal return, similar to radial return in computational plasticity. When the node is returned back to the surface, contact boundary conditions are again met. However, since this forced
850
K. Farahani et al. / Comput. Methods Appl. Mech. Engrg. 191 (2001) 843±860
Fig. 3. Rotation of target surface degrees of freedom.
normal return will cause additional unbalanced forces, simultaneous satisfaction of boundary conditions and equilibrium is obtained via equilibrium iterations. 4.3. Step 3: contact force, from contactor node to target surface After obtaining incremental displacements and rotating them back into the global coordinates, the internal force vector is computed. As a usual nonlinear FEM scheme, the dierence between external and internal force vector of all the elements originates subsequent incremental deformations. However, the unbalanced force of target element must be computed from a dierent scheme. That is, the succeeding step is to externally load the target surface by normal internal force of contactor node (for a frictionless contact), then compute its unbalanced forces, in each equilibrium iteration. Hence, in addition to geometric and material nonlinearities, change and movement of contact force are other sources of target element unbalanced forces. 4.4. Step 4: convergence control The iterations must continue until the convergence criteria is reached. During the iterations, the contactor node may move from one target element to another. This is important since the contact forces must be imposed on proper target surface even if the UE matrix is not updated. It is noted that the time of contact and departure of elements must be determined with enough accuracy. Also, since the contact iterations may cause repeated spurious material loading and unloading in elasto-plastic analyses due to the normal return algorithm, special attention must be paid to the stress update scheme. For example, the reference stress state may be taken the one at the beginning of the time step. The ¯owchart given in Appendix A manifests a full Newton±Raphson scheme of the UEM algorithm. 4.5. Geometry approximation It is seen that Eqs. (4.1) and (4.2) of step 1 have most importance in constitution of K U and they completely depend on the geometry and formulation of the applied elements. For higher-order elements or target surfaces, these relations will take more complicated forms. However, for such cases some geometric approximations may be adopted to simplify (4.1) and (4.2) so as to speed up the calculations with any desired accuracy. These approximations may be compensated in other steps and equilibrium iterations. As an example, Fig. 4 shows an eight node isoparametric element with a secondorder target surface/curve which is approximated with two lines. This will highly reduce the form of impenetration equation and the number of additional terms in the united element stiness matrix. Hence, instead of writing K U Unite K T ; K C ; a; b, depending on the contact position we may write
K. Farahani et al. / Comput. Methods Appl. Mech. Engrg. 191 (2001) 843±860
851
Fig. 4. Geometry approximation: (a) exact geometry; (b) approximate geometry.
Fig. 5. Higher-order target surface: (a) contact position; (b) degrees of freedom.
K U Unite K T ; K C ; ai ; bi where i depends on the contact zone shown in Fig. 4(b). This approximation which offers ¯exibility to the procedure to utilize higher-order elements is used only in step 1, while other steps may continue with exact element geometry. To show the effect of higher-order target surface on the complexity uniting equation, consider the parabolic target surface shown in Fig. (5). It is easy to show that uniting equation for the rotated degrees can be written as 6 X oY oY oY oY oY oY u12 cos a u1 u2 u3 u4 u5 u6 M i ui ; 4:10 oX1 oY1 oX2 oY2 oX3 oY3 i1 X X6
where Xi and Yi are coordinates of node number i. Also, Y is the equation of target curve in XY plain as Y aX 2 bX c where 2 3 2 2 a X1 4b5 4X2 2 c X32
X1 X2 X3
3 12 3 Y1 1 1 5 4 Y2 5 : Y3 1
4:11
As an example, M1 in Eq. (4.10) is
3
M1 A X6 where A
X1
X2 X6 X2 X2
X3 X2 X3 X3
Y2 X3 X1
Y1 Y3 X2 X1
Y1 Y2 X3 X2
X1 and X6 X3 b X1
X3 :
Y3 ; X3
4:12
852
K. Farahani et al. / Comput. Methods Appl. Mech. Engrg. 191 (2001) 843±860
Therefore, in such cases, a large amount of computation time must be spent to extract the uniting equation coecients.
5. Mathematical presentation of UEM Following the procedure of Section 4, the mathematical form of uniting equation in a three-dimensional space may be derived to unite any arbitrary contacted bodies with any number of contactor nodes. Assume that the potential energy function of the bodies shown in Fig. 1(b) before contact is P de®ned as P P u1 ; . . . ; un ;
5:1
where n is the total number of freedom degrees of all the bodies. Therefore, the stiffness matrix of the system can be de®ned as Kij
o2 P : oui ouj
5:2
After contact, a series of uniting equations or contact constraints such as uk F u1 ; . . . ; ui ; . . . ; un ;
i 6 k;
5:3
may hold, where uk 's are the normal d.o.f. of contactor nodes obtained after rotation. Also, all the arguments of contact constraint function F( ) correspond to d.o.f.'s of the target surfaces. Substituting (5.3) for the system as into (5.1) yields as new potential P P u1 ; . . . ; un P ; 5:4:a uk F u1 ;...;ui ;...;un or P u 1 ; . . . ; ui ; . . . ; un ; P
i 6 k:
5:4:b
Therefore, the new tangent united stiness matrix of the structure, K U , is de®ned as KijU
o2 P ; oui ouj
i; j 6 k:
5:5
Taking the partial dierentiation relative to i and j and one by one yields o oP o oP oP oF U Kij ; ouj oui ouj oui oF oui uk F u1 ;...;ui ;...;un
5:6
or KijU
o2 P o2 P oF oui ouj oui oF ouj
o2 P o2 P oF 2 ouj oF oF ouj
oF oP o2 F : oui oF oui ouj
5:7
Therefore, linearizing the uniting function F will both eliminate the last term of right-hand side of (5.5) and simplify its differentiation. After rearrangement and regarding (5.2) and (5.7) may be rewritten as 2 2 o P oF o2 P oF o P oF oF oP o2 F KijU Kij ; 5:8 oui oF ouj ouj oF oui oF 2 oui ouj oF oui ouj which thoroughly exhibits the quality and origin of additional terms in the united stiness matrix. It is noted that equations such as (2.3), (2.8), and (5.3) may alternatively be called the no-overlap, impenetration, contact, or uniting constraint.
K. Farahani et al. / Comput. Methods Appl. Mech. Engrg. 191 (2001) 843±860
853
6. Matrix presentation of UEM Again, after ®nite element discretization of the bodies, the static system potential functional before the contact is as given in (3.1). However, after the contact, the contactor nodes normal displacement vector UnC , can be described in terms of target surfaces deformation vector U T as UnC F U T ;
6:1
which again indicates the uniting equation. In general, (6.1) can be written as UnC AU T ;
6:2
where A is the connection matrix. The higher is the order of the target surface, the more complicated are the components of A, as shown in Section 4.5. The potential function after contact may be presented as P 12U t KU
U tf ;
6:3
where
6:4
Eq. (6.4) is simply a partitioning of stiness matrix and force and displacement vectors, into contactor normal, target surface, and the rest of freedom degrees, corresponding to C, T, and R right superscripts, respectively. Zero o-diagonals show that there is no physical relation between them before contact. After substitution of (6.1) into (6.3), we can write:
6:5 Again, noting the matrix oF =oU T A; considering dP 0 and rearranging the equations result in t cc T T U A K A K TT At K CR K TR f At fnC : 6:6 K CR A K TR K RR UR fR The stiness matrix in (6.6) is in fact what we called UE matrix or K U in the level of structure. It is noted that here, the number of rows and columns of connection matrix A is equal to the number of contactor nodes and the number of target surface d.o.f.'s, respectively. For example, the matrices A and K CC , corresponding to contact equation (4.2) are A hb sin a b cos a 1 h 0i C K CC K22 :
b sin a 1
b cos ai;
6:7 6:8
Although most of the discussions above concerned the stiness part of the contacting structures, they can be generalized to the problems including friction, damping, and inertial forces. After the solution of linear system of (6.6), UnC can be obtained from (6.2) for in®nitesimal deformations. However, for ®nite deformations, the four steps of Section 4 must be used to arrive in simultaneous equilibrium and boundary conditions satisfaction. Therefore, for in®nitesimal deformations, Eqs. (6.6) and (6.2) present a positive de®nite coecient matrix and exact satisfaction of boundary conditions, respectively, with no stiness illconditioning or addition in number of unknowns.
854
K. Farahani et al. / Comput. Methods Appl. Mech. Engrg. 191 (2001) 843±860
7. Numerical examples Based on the proposed procedure, a computer program is developed and static and dynamic numerical examples are given below. In the following examples, modi®ed Newton±Raphson iterative solution together with normal return scheme are used to solve the nonlinear equilibrium equations. Material is assumed to be elastic or elasto-plastic and rate-independent and the eect of friction is ignored. 7.1. Example 1 Fig. 6 shows two elastic beams each constituted of eight plain stress four node elements. The mechanical and geometrical properties are given in the ®gure. The total Lagrangian incremental formulation is chosen. The adopted material constitutive equation for ®nite elastic strains is Sij 2lEij kdij Ekk ;
7:1
where S and E are second Piola±Kirchhoff stress and Green±Lagrange strain tensors, respectively. k and l are material constants, d is the Kronecker delta, and index notation is adopted [28]. This contact problem shows a snap-through behavior. It is noted that if the strains become ®nite, the eect of using (7.1) on material response in ®nite elastic strains is important and is beyond the scope of this work. However. the use of dierent constitutive rate models which are based on objective rates of Cauchy stress (real stress) is well established. Fig. 7 shows three important phases of deformation, that is, contact, snap through, and departure. For the snap through part, a displacement control scheme is utilized to capture the history of motion. The contact is assumed only for the contactor node A shown in Fig. 6. Curves 1 and 2 of Fig. 8 show load±de¯ection curve of the point under load and contactor nodes, respectively. Curve 3 of the ®gure is the history of contact force. Points A, B, and C of curve 2 are relevant to three phases of Fig. 7. 7.2. Example 2 Fig. 9 shows two bodies again each constituted of 8 plain stress four node elements, with given geometric, mechanical, and initial conditions. Again the constitutive equation (7.1) is used to model the elastic
Fig. 6. Snap-through contact.
Fig. 7. Deformation phases.
K. Farahani et al. / Comput. Methods Appl. Mech. Engrg. 191 (2001) 843±860
855
Fig. 8. Load±de¯ection curves at dierent nodes.
Fig. 9. Dynamic contact.
part of material behavior. The lower body is massless and the upper one has concentrated masses Mx My 1 at each node except for the lower nodes. Gravity eect is not included. The chosen dynamic time integration scheme is Newmark with parameters d 0:5; a 0:25, and time step Dt 0:01 s. The upper body has an elastic±perfect plastic behavior with von Misez yield stress Fy 800, while the lower body is elastic. Fig. 10 shows the impact force history between the bodies, and Fig. 11 shows the deformation history of the impact. Plastic deformation in contactor element is apparent. 7.3. Example 3 Fig. 12 shows two bodies each constituted of plain stress four node elements, with given geometric and mechanical conditions. The hatched elements are much stier than the other two elements and the upper body is moving toward the lower one with initial velocity V0 10 and gravity acceleration G 10. Again the constitutive equation (7.1) is used for elastic material behavior. The upper body is elastic±perfect plastic with yield stress Fy 7:5. Newmark time integration with parameters d 0:5; a 0:25 and time step Dt 0:01 s is used. The height parameter h shown in the ®gure must be large enough to reduce the con®nement effect on the two soft elements, and in this example h 60. Due to relative stiffness of the elements, the problem is equivalent to the one shown in Fig. 13. Exact solution of the elastic±perfect plastic impact problem shown in Fig. 13 can be presented in six time intervals as
856
K. Farahani et al. / Comput. Methods Appl. Mech. Engrg. 191 (2001) 843±860
Fig. 10. Elastic and elasto-plastic impact force history.
Fig. 11. Deformation history.
Fig. 12. Elasto-plastic impact.
K. Farahani et al. / Comput. Methods Appl. Mech. Engrg. 191 (2001) 843±860
Fig. 13. Equivalent problem.
1 2 3
u 10; u 0 0; u_ 0:31 P t P 0 ) u 5t2
10; 10t:
u1 u 3:5804; t1 t 0:31; p p sin 200t1 0:05 cos 200t1 0:0643 P t1 P 0 ) u1 p13:1 200 u2 u 4:3305; t2 t 0:3743; 0:0615 P t2 P 0 ) u2 70t22 5:609t2 :
Fig. 14. Comparison of displacement history.
Fig. 15. Comparison of impact force history.
0:05:
857
858
K. Farahani et al. / Comput. Methods Appl. Mech. Engrg. 191 (2001) 843±860
Fig. A.1. Full Newton±Raphson iterative scheme with UEM and normal return.
4
u3 u 3:8452; t3 t 0:4358; p 0:1161 P t3 P 0 ) u3 0:7 cos 200t3
0:05:
K. Farahani et al. / Comput. Methods Appl. Mech. Engrg. 191 (2001) 843±860
5 6
859
U4 u 3:8452; t4 t 0:5519; 1:9748 P t4 P 0 ) u4 5t42 9:874t4 : U5 u 3:8452; t5 t 2:5267; p p sin 200t5 0:05 cos 200t5 0:2323 P t5 P 0 ) u5 p9:874 200
0:05:
Figs. 14 and 15 are the comparison between exact and numerical deformation and impact force history of the problem after two subsequent contact and departure, respectively. In both of the above ®gures, good agreement is observed between the exact solution and the numerical results. 8. Conclusions Regarding the incremental variational principles of contact problem, or based on the incremental form of the potential energy functional, the detailed mathematical and numerical procedure of UEM is presented in detail. The method focuses on node-to-surface contact/impacts which is the concern of most contact/ impact problems, and node-to-surface contact/impacts which is the concern of most contact/impact problems, and imposes the constraint of impenetration accurately without the drawbacks of Lagrange Multipliers and Penalty Function methods. That is, the number of unknowns and solution scheme do not change, and the boundary conditions are satis®ed exactly without loss of accuracy. The main idea is to omit normal d.o.f. of the contactors node by eliminating it from the system potential functional via uniting equations, and to use a normal return scheme to return the contactor node back to the target surface. The procedure can also be interpreted as treatment of the contact force between two bodies as the internal force of a whole body and can also be extended to frictional problems. The elemental, mathematical, and matrix presentation of the uniting step are also presented. Illustrative examples are given to demonstrate the capability and accuracy of the method. Enhancing or comparing this technique with other methods may be the matter of future works. Appendix A For a general nonlinear FEM-based contact/impact computer code, the four steps of UEM combined with a full Newton±Raphson iterative scheme may be illustrated in the ¯owchart of Fig. A.1. References [1] T.J.R. Hughes, R.L. Taylor, J.L. Sackman, A. Curnier, W. Kanoknukulchai, A ®nite element method for a class of contact± impact problems, Computer Methods in Applied Mechanics and Engineering 8 (1976) 249±276. [2] K.J. Bathe, A.B. Chaudhary, A solution method for static and dynamic analysis of three-dimensional contact problems with friction, Computers and Structures 24 (1986) 855±873. [3] N. Asano, An approximate hybrid type of virtual work principle for two elastoimpact contact bodies, Bulletin of JSME 26 (1983) 1849±1856. [4] N. Asano, A hybrid type of virtual work principle for impact contact problems of two bodies, Bulletin of JSME 29 (1986) 1679± 1684. [5] N.J. Carpenter, R.L. Taylor, M.G. Katona, Lagrange constraints for transient ®nite element surface contact, International Journal of Numerical Methods in Engineering 32 (1991) 103±128. [6] Z.-H. Zhong, L.G. Nilsson, Lagrange multiplier approach for evaluation of friction in explicit ®nite-element analysis, International Journal of Numerical Methods in Engineering 10 (1994) 249±255. [7] J.C. Simo, P. Wriggers, K.H. Schweizerhof, R.L. Taylor, Finite deformation post-buckling analysis involving inelasticity and contact constraints, International Journal of Numerical Methods in Engineering 23 (1986) 779±800. [8] N. Asano, A virtual work principle using penalty function method for impact contact problems of two bodies, Bulletin of JSME 29 (1986) 731±736. [9] N. Asano, A penalty function type of virtual work principle for impact contact problems of two bodies, Bulletin of JSME 29 (1986) 3701±3709. [10] Y. Kanto, G. Yagawa, A dynamic contact buckling analysis by the penalty ®nite element method, International Journal of Numerical Methods in Engineering 29 (1990) 755±774.
860
K. Farahani et al. / Comput. Methods Appl. Mech. Engrg. 191 (2001) 843±860
[11] R.F. Kulak, Adaptive contact elements for three-dimensional explicit transient analysis, Computer Methods in Applied Mechanics and Engineering 72 (1989) 125±151. [12] A.B. Pifko, R. Winter, Theory and application of ®nite element analysis to structural crash simulation, Computers and Structures 13 (1981) 277±285. [13] T. Belytschko, M.O. Neal, Contact±impact by pinball algorithm with penalty and Lagrangian methods, International Journal of Numerical Methods in Engineering 13 (1991) 547±572. [14] H.M. De La Fuente, C.A. Fellipa, Ephemeral penalty functions for contact±impact dynamics, Finite Elements in Analysis and Design 9 (1991) 177±191. [15] I. Hunek, On a penalty formulation for contact±impact problems, Computers and Structures 11 (1993) 193±203. [16] N. Kikuchi, A smoothing technique for reduce integration penalty methods in contact problems, International Journal Numerical methods Engineering 18 (1982) 343±350. [17] T. Campos, J.T. Oden, N. Kikuchi, A numerical analysis of a class of contact problems with friction in elastostatics, Computer Methods in Applied Mechanics and Engineering 34 (1982) 821±845. [18] B. Nour-Omid, P. Wriggers, A two level iteration method for solution of contact problems, Computer Methods in Applied Mechanics and engineering 54 (1986) 131±144. [19] J.C. Simo, P. Wriggers, R.L. Taylor, A perturbed Lagrangian formulation for the ®nite element solution of contact problems, Computer Methods in Applied Mechanics and Engineering 50 (1985) 163±180. [20] P. Wriggers, W. Wagner, E. Stein, Algorithms for nonlinear contact constraints with application to stability problems of rods and shells, Computational Mechanics 2 (1987) 215±230. [21] S.H. Lee, Rudimentary considerations for adaptive gap/friction element based on the penalty method, Computers and Structures 47 (1993) 1043±1056. [22] J. Ghaboussi, E.L. Wilson, J. Isenberg, Finite element of rock joints and interfaces, Journal of Soil Mechanics and Foundations Division, ASCE 99 (1973) 833±848. [23] S. Desai, M.M. Zaman, J.G. Lightner, H.J. Siriwardane, Thin-layer element for interfaces and joints, International Journal of Numerical and Analytical Methods in Geomechanics 8 (1984) 19±43. [24] P. Kowalczyk, Finite±deformation interface formulation for frictionless contact problems, Communications in Numerical Methods in Engineering 10 (1994) 879±893. [25] J.C. Simo, T.A. Laursen, An augmented Lagrangian treatment of contact problems involving friction, Computers and Structures 42 (1992) 97±116. [26] J.-H. Heegaard, A. Curnier, An augmented Lagrangian method for discrete large-slip contact problems involving friction, International Journal of Numerical Methods in Engineering 36 (1993) 569±593. [27] T.A. Laursen, S. Govindjee, A note on the treatment of frictionless contact between non-smooth surfaces in fully-nonlinear problems, Communications in Numerical Methods in Engineering 10 (1994) 869±878. [28] K.J. Bathe, Finite Element Procedures, Prentice-Hall, Englewood Clis, NJ, 1996. [29] K. Farahani, M. Mo®d, A. Vafai, A solution method for general contact-impact problems, Computer Methods in Applied Mechanics and Engineering 187 (2000) 69±77.