The Symmetric Boundary Element Method based on the Galerkin hypotheses has found
application in the nonlinear analysis of plasticity and contact-detachment problems, but dealt with separately.
In this paper we wants to treat these complex phenomena together.
This method works in structures by introducing a subdivision into sub-structures, distinguished into macroelements,
where elastic behaviour is assumed, and bem-elements, where it is possible for plastic strains to
occur. In all the sub-structures, elasticity equations are written and regularity conditions in weighted (weak)
form and/or in nodal (strong) form between boundaries have to be introduced, to attain the solving equation
system

PANZECA, Teotista

TERRAVECCHIA, Silvio Salvatore

ZITO, Liborio

Archivio istituzionale della ricerca - Università di Palermo

MULTIDOMAIN SBEM ANALYSIS OF TWO DIMENSIONAL ELASTOPLASTIC-CONTACT PROBLEMS Panzeca T.1,a, Terravecchia S. 1,b and Zito L.1,c 1Dip.to di Ingegneria Civile, Ambientale e Aerospaziale, University of Palermo, Italy. atpanzeca@tiscali.it, bsilvioterravecchia@gmail.com, czitoliborio@tiscali.it. Keywords: symmetric BEM; elastoplasticity; contact-detachment; substructuring. Abstract: The Symmetric Boundary Element Method based on the Galerkin hypotheses has found application in the nonlinear analysis of plasticity and contact-detachment problems, but dealt with separately. In this paper we wants to treat these complex phenomena together. This method works in structures by introducing a subdivision into sub-structures, distinguished into macro-elements, where elastic behaviour is assumed, and bem-elements, where it is possible for plastic strains to occur. In all the sub-structures, elasticity equations are written and regularity conditions in weighted (weak) form and/or in nodal (strong) form between boundaries have to be introduced, to attain the solving equation system. Introduction The present paper shows a strategy to perform the elastoplastic-contact-detachment analysis by using the Symmetric Boundary Element Method based on the Galerkin (SGBEM) hypotheses in a simultaneous analysis obtained by solving the non-linear problems of elastoplasticity and contact-detachment using Linear Complementarity Problem (LCP) in an incremental approach. - Plasticity problems. In the plastic analysis carried out using the symmetric BEM, it is necessary to distinguish a computing phase for the elastic response to all the actions, including the volumetric (body forces and plastic strains) ones, and a subsequent one for plastic strain evaluation, stored during the loading process [1]. In a first phase the Somigliana Identities (SIs) of the displacements and the tractions, both evaluated on the boundary, are employed through a weighting process. In a subsequent phase the stresses have to be evaluated in each bem-element and a predictor-corrector process has to be performed in order to evaluate the plastic strains stored in the bem-elements where the stress violates the elastic yield domain. In both phases, strong singular integrals are involved in the domain integrals when stresses and tractions caused by volumetric actions have to be evaluated [1,2]. - Contact-detachment problems. On the basis of the boundary integral method, in its symmetric formulation, the frictionless unilateral contact between two elastic bodies was studied according to the Signorini formulation [3]. A boundary discretization by boundary elements of the two bodies in contact leads to an algebraic formulation in the form of linear complementarity problem. - Elastoplasticity and Contact-detachment. The analysis of two bodies in contact having elastoplastic behaviour can be performed simultaneously, using an LCP analysis by alternating the contact-detachment phenomenon with the plasticity. This proves to be advantageous when this analysis is carried out through the symmetric BEM, mainly for two reasons: The contact-detachment process proves to have immediate execution because it is carried out through comparison between generalized quantities evaluated along the boundary elements and the reference values in weighted form. At every step characterizing the previous phase, an elastoplastic analysis is made in accordance with the predictor-corrector strategy. 1. The equation system governing the elastoplastic-contact problems In this section the integral equations governing the elastoplastic-contact problems are shown.  Let us consider a bi-dimensional body having domain Ω  and boundary Γ , subjected to actions acting in its plane: - forces 2f  at the portion 2Γ  of free boundary, - displacements 1u  imposed at the portion 1Γ  of constrained boundary, - body forces b  and plastic strains ε p  in Ω . The external actions 2 1, ,f u b  may increase separately or simultaneously through the multiplier β . In the hypothesis that the physical and geometrical characteristics of the body are zone-wise variables, an appropriate subdivision of the domain into bem-elements is introduced. This subdivision involves the introduction of an interface boundary 0Γ  between contiguous bem-elements and, as a consequence, two new unknown quantities arising in the analysis problem, i.e. the displacements 0u  and the tractions 0t  vectors, both referring to interface boundaries.  Let us start by imposing for each bem-e the classical SIs: ( )( )( )Γ Γ Ω ΩΓ Γ Ω ΩΓ Γ Ω Ω= Γ + − Γ + Ω+ Ω= Γ + − Γ + Ω+ Ω= Γ + − Γ + Ω+ Ω∫ ∫ ∫ ∫∫ ∫ ∫ ∫∫ ∫ ∫ ∫u G f G u G ε G bt G f G u G ε G bσ G f G u G ε G bpuu ut u uuptu tt t tupu t ud d d dd d d dd d d dσσσ σ σσ σ (1a-c)  These provide the displacements, tractions and stresses in the unbounded domain caused by layered mechanical jumps f  and double-layered kinematical ones (−u ) as well as by volumetric actions ε p  (volumetric distortions) and b  (body forces) both in the Ω  domain. The operators G pq  are the Fundamental Solution matrices, whose symbolism was introduced by Maier and Polizzotto [4]; the sub-indices p=u,t,σ  and q=u,t,σ  indicate the effect and the dual quantity in an energetic sense associated with the cause, respectively. For each substructure we can write the following integral eqautions 0 00 0 00 0 0 00 0ˆˆˆ= + + += + + += + + +0 AX A X A p LZ A X A X A p Lσ a X a X a p lσσσ σ σσ σβββ  (2a-c) These are obtained by weighting all the coefficients of eqs.(1a,b) computed on the boundaries according to the Galerkin approach [3-5]. The vector X  collects the sub-vectors 1F , ( )2U− , whereas the vector 0X  collects the sub-vectors 0F , ( )0−U  along the interfaces. The vector 0Z  collects the generalized (or weighted) diplacement and traction vectors defined in the interface boundary elements, obtained as a weighted response to all the known actions, amplified by β , and unknown actions, regarding boundary and domain quantities. The vector σ  represents the stress, evaluated at the Gauss points, due to all the known actions, amplified by β , and unknown actions. The vector p  represents the plastic strains defining the plastic distribution =ε Ψ pp p  inside each bem-e. By performing a variable condensation through the replacement of the X  vector extracted from eq.(2a) into eqs.(2b,c), one obtains: 10 00 00 0 0 00 0ˆ[ ]ˆˆ−= − + += + += + +X A A X A p LZ D X D p Zσ d X d p σσσσ σσβββ  (3a-c) The latter are the characteristic equations of each bem-e. Because the body is subdivided into m bem-elements, for each of these the eqs.(3a-c) can be written. Thus we obtain three global relations connecting all the generalized quantities and the stresses related to the bem-elements considered, formally equal to eqs.(3a,c), but regarding the characteristic equations of the system. Let us introduce the strong and weak coupling conditions between adjacent bem-elements: 0 00   =       =X EξEZ 0strong coupling condtionsweak coupling condtions  (4a,b) 0ξ  being the nodal interface vector which collects the mechanical 0F  and kinematical ( )0−U  unknowns of the assembled system. Using eqs.(4a,b), eqs.(3a-c) become: 10 000 0 0 00 0ˆ[ ]ˆˆ−= − + ++ + == + +X A A Eξ A p LK ξ K p f 0σ k ξ k p σσσσ σσβββ   (5a-c) Let us perform a new variable condensation through the replacement of 0ξ  vector extracted from eq.(5b) into eqs.(5a,c), thus obtaining: 10 010 00 0 0ˆ[ ]ˆ[ ]ˆ−−= − + += − += +X A A Eξ A p Lξ K K p fσ Kp σeσσβββ  (6a-c) Eqs.(6a,b) provide the elastic solution in terms of nodal forces 1F , 0F  and in terms of nodal displacements 2−U , 0−U . In detail, the quantities 0F  on 0Γ  and 2−U  on 2Γ  govern the contact-detachment phenomenum. Eq.(6c) provides the stress at the strain points of each bem-e as a function of the volumetric plastic strain p  and of the external actions σˆe , the latter amplified by β . The matrix K , defined as the self-stress influence matrix of the assembled system, is a square matrix having 3mx3m dimensions, with m bem-elements. It is fully-populated, non-symmetric and semi-definite negative. The evaluation of this matrix only involves knowledge of the material elastic characteristics and of the structure geometry.  This equation is used to evaluate the trial stress in the predictor phase, whereas the first term is utilized to perform the corrector phase, in order to obtain the collapse load factor. The reader can refer to Zito et al. [5] for a more detailed discussion of the characteristics of this equation introduced for a multidomain SGBEM problem. 2. The incremental contact-detachment algorithm Inside the topic of the SGBEM, in order to reach the analytical solution to this frictionless contact-detachment problem, an iterative LCP procedure can be employed once the incremental elastic analysis has been performed using eqs.(6a,b). For this purpose we remember that the unknown vectors ( )0 1F n+ , ( )2 1UAn+  and ( )2 1UBn+ , obtained by eqs.(6a,b) refer to the nodes of the in-contact boundary and to the nodes of the detached one at the generic load increment 1+n . Indeed, the vector 0 0 0= = −F F FA B  represents the nodal forces of the body A, computed in the contact boundary 0ΓA  and the vectors 2UA  and 2UB  represent the nodal displacements of the free boundaries of  2ΓA  and 2ΓB  respectively. With reference to the system of the two in-contact bodies, whose boundaries are discretized into boundary elements, the contact-detachment phenomenon can be computed by rewriting in discrete form the classical Signorini equations via SGBEM [3] ( ) ( )( ) ( )( )2 1 2 1 1− − ≤N U U H 0 A A B2 n+ n+ n+  gap condition  (7a) ( )0 0 1 − ≤N F C 0 A n+  contact condition  (7b) ( ) ( )( ) ( )( ) ( )2 1 2 1 1 2 1 0⎡ ⎤ ⎡ ⎤− − =⎢ ⎥ ⎣ ⎦⎣ ⎦N U U H N F A A B A2 2n+ n+ n+ n+               complementarity condition on 2Γ  (7c) ( ) ( )( ) ( )0 00 1 0 1 0 1 0⎡ ⎤ ⎡ ⎤− − =⎣ ⎦⎣ ⎦N U U N F C A A B An+ n+ n+                         complementarity condition on 0Γ  (7d) where 0 0(   )=N n" "A Adiag  and 2 2(   )=N n" "A Adiag  are global matrices collecting the normal vector associated with the boundary 0ΓA , 2ΓA  of the body A. The vector ( )1H n+  collects all the nodal gaps between the corresponding nodes of the boundaries 2ΓA  and 2ΓB  at the load increment 1+n , in the zone of potential contact, whereas the vector C  collects the cohesion between the nodes which are in contact, in the zone of potential detachment 0Γ . In detail, through the eq.(7a), all the nodes of the free boundary 2Γ , where the condition ( ) ( )( ) ( )( )2 2− − =N U U H 0 A A B2 n+1 n+1 n+1  occurs, change into the contact boundary 0Γ , thus defining a new contact boundary. Vice versa, through eq.(7b), all the nodes of the contact boundary 0Γ , where the condition ( )0 0 − =N F C 0 A n+1  occurs, change into the free boundary 2Γ , thus defining a new detachment boundary. 3. The incremental elastoplastic analysis for active macro-zones A brief description of the strategy utilized for incremental elastoplastic analysis via Multidomain SGBEM, called elastoplastic active macro-zone analysis, is provided in this section. The complete version can be found in [5]. For each loading step and at each bem-e this analysis uses eq.(6c) both to evaluate the trial stresses in the predictor phase and to compute the plastic strains in the corrector phase. Let us start by evaluating the trial stresses, i.e. the purely elastic response at the load increment 1+n  in each m bem-element of the discretized body. For this purpose eq.(6c) provides all the predictors ( )*1+σ n  as a function of the plastic strain vector ( )p n , stored up at step n , of all the increments ∆p  inside step 1+n , and of the external load σˆe , amplified by ( )1+nβ : ( ) ( ) ( ) ( ) ( )*1 1 1 1ˆ          + + + += + = + ∆σ K p σ p p pen n n n nβ with   (8) where ( ) ( )1+ = + ∆n nβ β β  is the load factor and K  matrix is fully-populated and regards all the bem-elements, the place of nonlinear phenomenon. A check on the plastic consistency condition of the stresses, computed at appropriately chosen points inside each bem-e, is performed using the yield condition expressed in this context through the von Mises law, i.e.: ( ) ( ) ( )* 2121 1 1[ ] σ 0         1...+ + += − ≤ =σ σ σT yi n i n i nF with i mM  (9) where M  is a matrix of constants and σ y  the uni-axial yield stress. In the a bem-elements (with ≤a m ) where this inequality is violated, a return mapping phase is required in order to evaluate the plastic strain increments. This phase, called the corrector phase, uses the first term of eq.(6c) to obtain the elastoplastic solution at every bem-e where the plastic consistency condition is violated. In this phase the vector σ , representing the end step stress, as well as the increment of the volumetric plastic strain vector ∆p , are unknown quantities. The latter are the plastic strains to impose on every plastically active bem-element in order to have the stress on the yield boundary of the elastic domain, through which the direction of the plastic flow can be defined. Obviously, inside each loading step the corrector phase has to be repeated until all the predictors satisfy the plastic consistency conditions. In detail, the elastoplastic algorithm allows one to write, for all the active h bem-elements ( h = 1,...,a ), a nonlocal system at the 1+n  load step simultaneously in all the plastically active macro-zones identified in the previous predictor phase, i.e.: *( 1) ( 1) ( 1)+ + += + ∆σ σ K pa n a n aa a n   (10) ( 1) ( 1) ( 1) ( 1)[ ] , , [ ] =+ + + +≤ ∆ ≥ ∆σ 0 Λ 0 Λ σ 0a n a n a n a nF F   (10a-c) where eqs.(10a-c) are the plastic admissibility conditions for the a bem-elements. The index a characterizes the vectors and matrices connecting the mechanical and kinematical quantities relating to all the active a bem-elements. The Kaa  matrix coefficients are derived from the K  matrix present in eq.(6c), by extracting the blocks relating to the a plastically active bem-elements. In the following equations the subscript 1+n  has been omitted for simplicity. In the hypothesis that, for each h-th bem-e, the shape function defined in =ε Ψ pph p h  is the same as the shape function relating to the plastic multiplier, i.e. ψ Λ∆ = ∆h p hλ  with ψ 0≥p , the plastic strain increment for the a-th active bem-elements is expressed as:  [ ]∆ = ∆ ∂ = ∆σp Λ σ Λ σaa a a a a aF M   (11) The solving nonlinear system for all the active a bem-elements is the following: *[ ]⎧ − − ∆ =⎪⎨ =⎪⎩σ σ Λ K σ 0σ 0a a a aa a aaMF  (12) or, in explicit form, using the von Mises yield law and the plastic flow rule given by eq.(11) **212212σ 0σ 0⎧ − − ∆Λ − ∆Λ =⎪⎪⎪ − − ∆Λ − ∆Λ =⎪⎨ − =⎪⎪⎪⎪ − =⎩σ σ K σ K σ 0σ σ K σ K σ 0σ σσ σ…#…#1 1 1 11 1 a 1a aa a 1 a1 1 a aa aT1 1 yTa a yM MM MMM   (13) where σa  is the stress solution located on the yield surface of the elastic domain of all the active bem-elements, *σa  the elastic predictor, and ∆Λ K σa aa a aM  the corrective components (containing local and nonlocal contributions). The approximate solution of this nonlinear problem involving all the plastically active bem-elements, in terms of σa  and ∆Λa , can be obtained by applying the Newton-Raphson procedure. 4. The elastoplastic-contact-detachment procedure In this section the sequence of steps concerning the proposed procedure for obtaining the numerical results is described. Step 0: Imput data. Step 1: Load update ( ) ( )1+ = + ∆n nβ β β . Step 2: Nodal unknown update:  ( ) ( ) ( )( ) ( ) ( ) ( )( )( ) ( ) ( ) ( )( ) ( )( )10 00 0 01 1 0 1 0 1 0 110 01 1 1 2 1ˆˆ              −+ + + + +−+ + + +⎧ = − + → = →⎪⎨ = − − − →⎪⎩ξ K K p f X Eξ FX A A Eξ A p L Un n n n n nn n n n nββσσ. Step 3: Check on contact-detachment ( ) ( )( ) ( )( ) ( )02 1 2 1 1 0 1,+ + + +− − ≤ − ≤N U U H 0 N F C 0 A A B A2 n n n n . If true then go to Step 4 If false then modify topologically some boundaries 0 2 2 0  or/and  Γ →Γ Γ →Γ  and go to Step 2 Step 4: Computing the elastic predictor ( ) ( ) ( )*1 1 1 ˆ+ + += +σ K p σsn n nβ . Step 5: Check for yielding If  ( )*1[ ]    with 1...+ ≤ =σi nF Tol i m  then go to Step 1, If  ( )*1[ ]    with 1...+ > =σi nF Tol i m  then go to Step 6. Step 6: Identification of the active bem-elements (a being the active macro-zone) in the corrector phase, *[ ]⎧ − − ∆ =⎪⎨ =⎪⎩σ σ Λ K σ 0σ 0a a a aa a aaMF Step 7: Plastic strain vector update:  ( ) ( ) ( ) ( )1 1 1+ + += + ∆p p Λ σaa n a n a n a nM ,  go to Step 1   5. Numerical results In order to show the efficiency of the proposed method, a numerical test was performed in the hypothesis of elastic-contact-detachment only. Let us consider the detachment problem regarding a beam A supported by two elastic blocks B, without friction or sliding, symmetrically loaded. The analysis is performed on half the structure, as shown in Fig.1a. 2010 AB10 10 20qa) b) Fig.1. Beam supported on elastic blocks: a) geometric description,  b) strained shape obtained by iterative LCP analysis.  The geometrical and mechanical characteristics of the structure are the same as those utilized by Vodicka [6]. Thus the beam A, having unitary thickness, is characterized by the Young modulus 4aE 30.6x10 Mpa=  and Poisson ratio 0.3υ =  and is subjected to vertical force distribution q 1020 daN / m= . The body B is characterized by the Young modulus 6bE 30.6x10 Mpa=  and by same thickness and Poisson coefficient. In order to discretize the free and constrained boundaries of the solids A and B, a step p 2 cm=  and p 0.1cm=  were introduced, respectively on 2Γ  and 1Γ  boundaries. Fig.1b shows the final strained shape obtained by iterative LCP.  In Table 1 the detachment length is shown in comparison with the solution obtained by Vodicka [6]. In all the cases shown in Table 1 the detachment length proves to be very similar.  Method Detachment length [cm]SGBEM Iterative LCP 8.4  BEM (R. Vodicka [6]) 8.3  Table 1: Comparison of detachment lengths. References. [1] Panzeca T., Terravecchia S., Zito L., (2010). “Computational aspects in 2D SBEM analysis with domain inelastic actions”. Int. J. Numer. Meth.Engng., 82, 184-2004. [2] Gao X.W., (2002). “The radial integration method for evaluation of domain integrals with boundary-only discretization”. Engng. Anal. Bound. Elem., 26, 905-916. [3] Panzeca T., Salerno M., Terravecchia S., Zito L., (2008). “The symmetric boundary element method for unilateral contact problems”. Comput. Methods Appl. Mech. Engrg., 197, 2667–2679. [4] Maier G., Polizzotto C., (1987). “A Galerkin approach to boundary element elastoplastic analysis”. Comput. Methods Appl. Mech. Engrg., 60, 175-194. [5] Zito L., Cucco F., Parlavecchio E., Panzeca T., (2012). “Incremental elastoplastic analysis for active macro-zones”. Int. J. Numer. Meth.Engng., In press. [6] Vodicka R., (2000). “The first and the second-kind boundary integral equation systems for solution of frictionless contact problems”. Eng. Anal. Bound. Elem., 24, 407–426. 

Multidomain SBEM analysis for two dimensionalelastoplastic-contact problems

https://iris.unipa.it/retrieve/handle/10447/64414/45226/BETEQ%202012.pdf

Multidomain SBEM analysis for two dimensionalelastoplastic-contact problems

Abstract

Similar works

Full text

Available Versions

Archivio istituzionale della ricerca - Università di Palermo