This contribution is a continuation of our contribution denoted as PART I, where the discretized contact problem for elasto-perfectly plastic bodies was studied and suitable numerical methods were introduced. In particular, frictionless contact boundary conditions and Hencky’s material model with the von Mises criterion are considered. Here we describe some implementation details and present several numerical examples

Cermak, Martin

Haslinger, Jaroslav

Sysala, Stanislav

English

UPCommons

Numerical solution of perfect plastic problems with contact: part II - numerical realization

UPCommons. Portal del coneixement obert de la UPC

XII International Conference on Computational Plasticity. Fundamentals and ApplicationsCOMPLAS XIIE. On˜ate, D.R.J. Owen, D. Peric and B. Sua´rez (Eds)NUMERICAL SOLUTION OF PERFECT PLASTICPROBLEMS WITH CONTACT: PART II - NUMERICALREALIZATIONMARTIN CERMAK∗,‡, JAROSLAV HASLINGER†,‡, STANISLAV SYSALA‡∗VSˇB–Technical University of OstravaIT4Innovations17. listopadu 15, 708 33 Ostrava, Czech Republice-mail: martin.cermak@vsb.cz, web page: http://www.vsb.cz†Charles University in PragueDepartment of Numerical Mathematics, Faculty of Mathematics and PhysicsSokolovska´ 83, 186 75 Prague, Czech Republice-mail: hasling@karlin.mff.cuni.cz, web page: http://www.mff.cuni.cz‡Institute of Geonics AS CR, v.v.i.Department of IT4InnovationsStudentska´ 1768, 708 00 Ostrava, Czech Republice-mail: stanislav.sysala@ugn.cas.cz, web page: http://www.ugn.cas.czKey words: Perfect plasticity, Contact, Domain decompositionAbstract. This contribution is a continuation of our contribution denoted as PART I,where the discretized contact problem for elasto-perfectly plastic bodies was studied andsuitable numerical methods were introduced. In particular, frictionless contact boundaryconditions and Hencky’s material model with the von Mises criterion are considered. Herewe describe some implementation details and present several numerical examples.1 INTRODUCTIONThis contribution is a continuation of PART I ([2]). We briefly describe some imple-mentation details and present several numerical examples. Some of the results will appearin [3].Algorithms ALG1-ALG4 introduced in PART I ([2]) have a common feature: in eachiteration it is necessary to solve a quadratic programming problem. For numerical real-ization of this inner problem, we use a combination of the TFETI domain decompositionmethod [5] and the SMALSE-M method [6], similarly as in [1, 4]. The whole contactproblem for elastic-perfectly plastic bodies is implemented in MatLab within the MatSol1Numerical solution of perfect plastic problems with contact: part II - numerical realization  999Martin Cermak, Jaroslav Haslinger, Stanislav Sysalalibrary [8]. Another approach for numerical realization of contact problems for elasto-plastic bodies can be found in [10].This contribution is organized as follows. The algebraic formulation of the discretizedproblem is introduced in Section 2. In Section 3, some details concerning numericalrealization of the inner problem are mentioned. In Section 4, numerical results of sev-eral model examples for elasto-perfectly plastic problems (with and without contact) arepresented. Concluding remarks can be found in Section 5.2 ALGEBRAIC FORMULATION OF THE PROBLEMIn accordance with PART I ([2]), we consider two polyhedric domains Ω1,Ω2 ⊆ R3,whose boundaries are decomposed as follows: ∂Ωj = Γju ∪ Γjf ∪ Γjc, j = 1, 2, where Γju,Γjf , Γjc are open and mutually disjoint. On Γju = ∅, the structure is fixed, while surfacetractions are applied on Γjf , j = 1, 2. Further, Γc := Γ1c = Γ2c is a bounded contact zone,where the frictionless contact boundary conditions are prescribed.Since we apply the TFETI domain decomposition method [5], we decompose Ω1 andΩ2 into disjoint polyhedric subdomains Ωj,p ⊂ Ωj, p = 1, 2, . . . , sj, j = 1, 2, see Figure1. Continuity of a solution between two adjacent subdomains and the satisfaction ofthe Dirichlet condition will be taken into account by adding equality constraints into thedefinition of K, see below.Figure 1: Scheme of the geometry and domain decompositionLet T j,ph be a regular partition of Ωj,p, j = 1, 2, p = 1, 2, . . . , sj, into tetrahedronssuch that the nodes on the contact and common parts of boundaries of two adjacentsubdomains coincide. On any Ωj,p, we construct the spaces of continuous and piecewiselinear displacements over T j,ph . Thus any global displacement vector v ∈ Rn has thefollowing structurev =(vT1,1,vT1,2, . . . ,vT1,s1,vT2,1, . . . ,vT2,s2)T,21000Martin Cermak, Jaroslav Haslinger, Stanislav Sysalawhere vj,p denotes the displacement vector on Ωj,p, j = 1, 2, p = 1, 2, . . . , sj. Furthermore,the set of kinematically admissible displacement vectors is given byK := {v ∈ Rn | BEv = o, BIv ≤ o} .By means of the matrix BE ∈ RmE×n, we express the gluing conditions among neighboringsubdomains and the satisfaction of the homogeneous Dirichlet boundary condition, whilethe matrix BI ∈ RmI×n is used to express the non-penetration condition on Γc.The algebraic formulation of problem (P)h,ζ , ζ ≥ 0, introduced in PART I ([2]) readsas follows:find u ∈ K :{Jζ(u) ≤ Jζ(v) ∀v ∈ K,Jζ(v) := Ψ(v)− ζfTv, v ∈ Rn.Here f =(fT1,1, . . . , fT1,s1, fT2,1, . . . , fT2,s2)T ∈ Rn is the load vector corresponding to the appliesvolume and surface forces. The function Ψ is the algebraic form of the inner energyfunctional and has the following structure:Ψ(v) =(Ψ1,1(v1,1)T , . . . ,Ψ1,s1(v1,s1)T ,Ψ2,1(v2,1)T , . . . ,Ψ2,s2(v2,s2)T)T.Similarly, one can rewrite problems (P)∗h,ζ , (S)h,ζ,r, (S)αh , (P)αh introduced in PART I([2]) into the algebraic form.3 NOTES TO NUMERICAL REALIZATION OF INNER PROBLEMSAs we have mentioned, algorithms ALG1-ALG4 introduced in PART I ([2]) lead to thequadratic programming problems in each iteration. The scheme of these inner subprob-lems solved in the k-th iterative step is the following:find uk ∈ Kk : Jk(uk) ≤ Jk (v) ∀v ∈ Kk, (1)whereJk(v) :=12vTKρkv − fTk v, v ∈ Kk, (2)andKk :={v ∈ Rn ; B˜Ev = cE,k, BIv ≤ cI,k}.HereKρk ∈ Rn×n is a symmetric and positive semidefinite matrix representing the functionΣo,ρ(uk), ρ ∈ [0, 1], introduced in PART I ([2]). Let us recall that Σo,1(uk) = C, i.e.we obtain the elastic stiffness matrix for ρ = 1 (see ALG2). The computed vector ukrepresents either δuk in ALG1, ALG3, ALG4, or uk in ALG2. The form of the load vectorfk depends on the particular algorithm, but always it contains f . For ALG1, ALG2,ALG4, it holds that B˜E = BE and cE,k = o. For ALG3, the contraints L(v) = α isadditionally included into B˜Ev = cE,k.31001Martin Cermak, Jaroslav Haslinger, Stanislav SysalaThe TFETI domain decomposition method has been analyzed e.g. in [5, 6, 7]. Allconstraints appearing in the definition ofKk will be released by using Lagrange multipliers,see Figure 1. In particular, we use two types of the Lagrange multipliers, namely λI ∈RmI , λI ≥ o related to the inequality constraints, and λE ∈ Rm˜E related to the equalityconstraints. To simplify our notation, we denoteλ =[λEλI], B =[B˜EBI], ck =[cE,kcI,k],andΛ = {λ = (λTE,λTI )T ∈ Rm˜E+mI : λI ≥ o}.Then the Lagrangian associated with (1) reads asLk(v,λ) =12vTKρkv − fTk v + λT (Bv − ck), v ∈ Rn, λ ∈ Λ. (3)Using convexity of Lk and of the constraints, we can use the classical duality theory toreformulate problem (1) to getJk(uk) = minv∈KkJk(v) = minv∈Rnsupλ∈ΛLk(v,λ) = maxλ∈Λinfv∈RnLk(v,λ). (4)The max-inf problem in (4) is solved by algorithm SMALSE-M [6]. This algorithm isbased on the active set strategy and it combines three steps: CG with preconditioningbased on orthogonal projectors, expansion, and proportioning.4 NUMERICAL EXPERIMENTSIn this section, we illustrate the efficiency of the proposed algorithms on several modelexamples. First, we will apply ALG1, ALG2 and ALG3 to a model contact problemfor two bodies. Then we will compare the loading paths computed by ALG1, ALG3 andALG4 for a plane strain problem considering only one body without the contact boundaryconditions.4.1 3D contact problemLet Ω1 and Ω2 be two (3000×1000×1000) blocks with the configuration seen in Figure1. Observe that Ω1, Ω2 are fixed along the left, and right lateral face, respectively. Thecontact zone occupies 80% of the area of the faces. The load L is represented by theconstant surface traction f = (0, 0, 100) acting on the upper face of Ω2 and correspondsto ζ = 1. Ω1 and Ω2 are decomposed into 24 subdomains. The mesh Th contains 6 000nodes and 18 434 tetrahedrons. Both bodies Ωi are made of the same elasto-perfectlyplastic material which is characterized by the Young modulus E = 206 900, the Poissonratio ν = 0.29 and γ = 450√2/3 representing the initial yield stress.41002Martin Cermak, Jaroslav Haslinger, Stanislav SysalaIn the first set of numerical experiments, we investigate convergence of ALG1 andALG2 in dependence on ζ, ρ, and r, respectively. In particular, we initiate ALG1,2 usingthe zero displacements, strain and stress vectors and choose the stopping criterion givenby the relative displacement error (see PART I [2]) with the tolerance u = 1e− 4.For ALG1, the dependence of number of iterations on ζ and ρ is shown in Table 1.Table 1: Number of iterations in dependence on ζ and ρ for ALG1ζ \ ρ 0.00 0.05 0.10 0.15 0.20 0.30 0.50 1.001.00 5 5 6 6 6 7 8 101.30 7 11 13 16 18 24 30 481.50 11 41 65 87 107 143 205 331One can observe that the smallest number of iterations is for ρ = 0 and this numberincreases with increasing ζ and ρ. On the other hand, if we use ALG1 without damping,then convergence for ρ = 0 is problematic for higher values of ζ. But the influence ofdamping turns out to be minimal for ρ > 0.If ALG2 is used, convergence is faster for r between 0.10 and 0.30 independently of ζ.Moreover for r < 0.05 or r > 0.50 and higher values of ζ, convergence is too slow.The stress and displacement fields obtained by ALG1 and ALG2 are very similar.Furthermore, the results obtained for other stopping criteria mentioned in PART I ([2])are also similar.In the second set of experiments we compare the loading paths obtained by ALG1,ALG2 and ALG3 for the fixed values r = 0.2, ρ = 0 and ρ = 0.1.For ALG1,2, we use adaptively chosen values of ζ:0 = ζ0 < ζ1 < ζ2 < . . . < ζi < . . . < ζN .The incrementsζi := ζi−ζi−1 decrease for increasing i and are determined in dependenceon convergence of ALG1,2. For any i = 1, 2, . . ., the corresponding displacements ui andthe values αi = L(ui) are computed. ALG1 is initiated by using the zero displacementsfor ρ = 0 and ui−1 for ρ = 0.1. Let us note that the initial choice ui−1 for ρ = 0 leads tothe loss of convergence for larger i. ALG2 is initiated by using the zero strain and stressvectors for any i = 1, 2, . . .. Further, we use the stopping criterion given by the relativedisplacement error with the tolerance u = 1e − 4. The loading process terminates ifαi > 7e9 for some i. This value of α is sufficiently large to estimate ζlim, see Figure 2.For ALG3 we use ρ = 0. In the loading process, we keep the constant incrementα = αel/2 ≈ 0.27e9, where αel is the maximal value of α, for which the solution to(P)αh coincides with the solution to the corresponding contact problem for elastic bodies,see the circle in Figure 2. Then α increases up to the value 7e9. For a fixed value of αi,51003Martin Cermak, Jaroslav Haslinger, Stanislav Sysalai = 1, 2, . . ., we use the same stopping criterion and the tolerance as in ALG1,2. AlgorithmALG3 is iniciated using ui−1 computed in the previous step.The comparison of the loading paths for all ALG1-ALG3 is shown in Figure 2.Figure 2: Comparison of loading paths.One can see that the curves practically coincide. They are increasing and ζlim ≈ 1.55.Notice that ALG3 converges here even for ρ = 0. This does not hold for ALG1.Finally, the stress and displacements fields at the end of the loading process are depictedin Figure 3 and 4, respectively. From Figure 3 on the left, one can detect plastic parts ofthe bodies, where the maximal value γ = 450√2/3 is attained. From the right picture,one can observe possible shear zones in the vicinity of the fixed parts of the bodies.4.2 Plane strain problem without contactThe configuration of a 2D domain Ω is as Figure 5. On the left and bottom parts of itsboundary, we prescribe the homogeneous Dirichlet boundary conditions only in the normaldirection, i.e. the symmetry conditions are prescribed. On the top of the domain, thesurface traction in the normal direction is applied. The elasto-perfectly plastic materialis considered and characterized by the Young modulus E = 206 900, the Poisson ratio61004Martin Cermak, Jaroslav Haslinger, Stanislav SysalaFigure 3: Von Mises stress for σ (left) and Cε(u) (right) at the end of the loading process.Figure 4: Displacements at the end of the loading process.ν = 0.29 and γ = 450√2/3.First, we compare the loading paths obtained by ALG1 (ρ = 0), ALG3 (ρ = 0.1) andALG4 (ρ = 0) for fixed ρ similarly as in the previous example. To this end, we use themesh with 7600 elements and 3927 nodal points. Again, the obtained curves practicallycoincide, see Figure 6.Notice that we could not use ALG3 with ρ = 0 for determining the whole loading pathdue to loss of convergence. Moreover, for ρ = 0.1, ALG4 is a little bit faster than ALG3.Secondly, we will compare the loading paths in dependence on the discretization pa-rameter h. To this end, we use ALG3 with ρ = 0.1 for three different meshes with 3960,7600 and 20664 elements. The corresponding loading paths for α ∈ 〈0, 1008〉 are shownin Figure 7. The computed values of ζ corresponding to α = 1008 are 1.077, 1.066, and71005Martin Cermak, Jaroslav Haslinger, Stanislav SysalafΩ1 910Figure 5: Geometry of the problem.Figure 6: Comparison of loading paths.1.055, respectively. It means that the values of the limit load parameter decrease withdecreasing h, which is in accordance with the expected theoretical results.Finally, we illustrate the stress and displacement states at the end of the loadingprocess, see Figure 8 and 9.81006Martin Cermak, Jaroslav Haslinger, Stanislav SysalaFigure 7: Comparison of loading paths in dependence on the discretization parameter.Figure 8: Von Mises stress for σ (left) and Cε(u) (right) at the end of the loading process.5 CONCLUSIONSThis contribution deals with numerical realization of ALG1-ALG4, which were intro-duced in PART I ([2]) as possible tools for solving elasto-perfectly plastic contact problems.91007Martin Cermak, Jaroslav Haslinger, Stanislav SysalaFigure 9: Displacement at the end of the loading process.The efficiency of the algorithms was tested on several model examples.We see that convergence of ALG2 depends on the parameter r appearing in the aug-mented Lagrangian. However it is not clear how to find an optimal value of r. Further,we analyzed how the regularization parameter ρ influences convergence of ALG1, ALG3and ALG4. The smaller value of ρ is, faster convergence is. On the other hand, we foundsituations, when the choice ρ = 0 led to loss of convergence for all of these algorithms.We also compared two different strategies for controlling the loading process: or by theload parameter ζ or by the parameter α representing the work of external forces. Bothstrategies have some advantages and disadvantages. For example, numerical realizationof (P)αh can be more complicated than for (P)h,ζ . On the other hand, (P)αh has a solutionfor any α ≥ 0 and thus controlling the process by α is simpler than by ζ.Acknowledgments. This work was supported by the European Regional DevelopmentFund in the IT4Innovations Centre of Excellence project (CZ.1.05/1.1.00/02.0070) andby the project SPOMECH - Creating a multidisciplinary R&D team for reliable solutionof mechanical problems, reg. no. CZ.1.07/2.3.00/20.0070 within Operational Programme’Education for competitiveness’ funded by Structural Funds of the European Union andstate budget of the Czech Republic.REFERENCES[1] Cermak, M.: Scalable algorithms for solving elasto-plastic problems. Ph.D. thesis,VSB-TU Ostrava (2012).[2] Cermak, M., Haslinger, J., Sysala, S.: Numerical solutions of perfect plastic problemswith contact: PART I - theory. Submitted as a contribution for Complas XII.101008Martin Cermak, Jaroslav Haslinger, Stanislav Sysala[3] Cermak, M., Haslinger, J., Kozubek, T., Sysala, S.: Discretization and numerical re-alization of contact problems for elastic-perfectly plastic bodies. PART II – numericalrealization, limit analysis. To appear.[4] Cermak, M., Kozubek, T., Sysala, S., Valdman J.: A TFETI Domain DecompositionSolver for Elastoplastic Problems. arXiv:1205.1926 [math.NA].[5] Dosta´l, Z., Hora´k, D., Kucˇera, R.: Total feti - an easier implementable variant ofthe feti method for numerical solution of elliptic pde. Communications in NumericalMethods in Engineering. (2006) 22: 1155–1162.[6] Dosta´l, Z., Kozubek, T.: An optimal algorithm and superrelaxation for minimiza-tion of a quadratic function subject to separable convex constraints with applications,Mathematical Programming. (2012) 135, 1-2: 195–220.[7] Dosta´l, Z., Kozubek, T., Vondra´k, V., Brzobohaty´, T., Markopoulos, A.: Scalabletfeti algorithm for the solution of multibody contact problems of elasticity. Interna-tional Journal for Numerical Methods in Engineering. (2012) 82: 1384–1405.[8] Kozubek, T., Markopoulos, A., Brzobohaty´, T., Kucˇera, R., Vondra´k, V., Dosta´l,Z.: Matsol - matlab efficient solvers for problems in engineering. Available fromhttp://matsol.vsb.cz/ (2012).[9] Sysala, S., Haslinger, J., Hlava´cˇek, I., Cermak, M.: Discretization and numerical re-alization of contact problems for elastic-perfectly plastic bodies. PART I – discretiza-tion, limit analysis. Submitted.[10] Wohlmuth, B.: Variationally consistent discretization schemes and numerical algo-rithms for contact problems. Acta Numerica (2011) 20: 569–734.111009

https://upcommons.upc.edu/bitstream/2117/182990/1/COMPLAS2013-92_Numerical%20solution%20of%20perfect%20plastic.pdf

Numerical solution of perfect plastic problems with contact: part II - numerical realization

Abstract

Similar works

Full text

Available Versions

UPCommons

UPCommons. Portal del coneixement obert de la UPC