In the present paper we compute upper and lower bounds for limit analysis in two and three dimensions. From the solution of the discretised upper and lower bound problems, and from the optimum displacement rate and stress fields, we compute an error estimate defined at the body elements and at their boundaries, which are applied in an adaptive remeshing strategy. In order to reduce the computational cost in 3D limit analysis, the tightness of the upper bound is relaxed and its computation avoided. Instead, the results of the lower bound are used to estimate elemental and edge errors. The theory has been implemented for Von Mises materials, and applied to two- and three-dimensions
examples.Peer Reviewe

Bonet Carbonell, Javier

Huerta, Antonio

Muñoz, José J.

Peraire Guitart, Jaume

In the present paper we compute upper and lower bounds for limit analysis in two and three dimensions. From the solution of the discretised upper and lower bound problems, and from the optimum displacement rate and stress fields, we compute an error estimate defined at the body elements and at their boundaries, which are applied in an adaptive remeshing strategy. In order to reduce the computational cost in 3D limit analysis, the tightness of the upper bound is relaxed and its computation avoided. Instead, the results of the lower bound are used to estimate elemental and edge errors. The theory has been implemented for Von Mises materials, and applied to two- and three-dimensions

examples.Peer Reviewe

UPCommons

Bounds and adaptivity for 3D limit analysis

UPCommons. Portal del coneixement obert de la UPC

7th Workshop on Numerical Methods in AppliedScience and Engineering (NMASE 08)Vall de Nu´ria, 9 a 11 de enero de 2008c©LaCa`N, www.lacan-upc.esBOUNDS AND ADAPTIVITY FOR 3D LIMIT ANALYSISJ. Mun˜oz1∗, J. Bonet2, A. Huerta1 and J. Peraire31: Dep. Applied Mathematics III, LaCa`NUniv. Poite`cnica de Catalunya (UPC)e-mail: j.munoz@upc.edu, web: http://www-lacan.upc.es2: Civil and Computational Engineering Centre, School of EngineeringUniversity of Wales, Swansea, UK3: Department of Aeronautics and AstronauticsMassachusetts Institute of Technology, USAKey words: adaptivity, SOCP, limit analysis, bounds, error estimatesAbstract. In the present paper we compute upper and lower bounds for limit analysisin two and three dimensions. From the solution of the discretised upper and lower boundproblems, and from the optimum displacement rate and stress fields, we compute an er-ror estimate defined at the body elements and at their boundaries, which are applied inan adaptive remeshing strategy. In order to reduce the computational cost in 3D limitanalysis, the tightness of the upper bound is relaxed and its computation avoided. Instead,the results of the lower bound are used to estimate elemental and edge errors. The theoryhas been implemented for Von Mises materials, and applied to two- and three-dimensionsexamples.1 INTRODUCTIONDuring the last decade, numerical methods for limit analysis have undergone greatprogress. This is mainly due to the recent development of fast optimisation techniques[Stu99, TTT03], and the increase of the computational capacity, which allow to afford thesolution of realistic problems.The models described in the literature may be mainly characterised by three aspects:(i) the use of strict bounds or estimates of the load factor, (ii) the optimisation techniquesused to solve the discretised max-min problem, and (iii) the use of an adaptive remeshingstrategy. We will here show two discretisations of the saddle-point problem of limitanalysis, one for the bound problem and another for the upper bound problem, whichcan be solved as a Second Order Cone Program (SOCP). Each one of them yield strictbounds, and yield optimum values of the displacement rate and stress field that are usedto design elemental and edge contributions to the bound gap. These values are in turnemployed in an adaptive remeshing strategy.1J. Mun˜oz1∗, J. Bonet2, A. Huerta1 and J. Peraire3Similar developments have been presented in reference [MBHP] for two-dimensionalproblems in Von Mises and Mohr-Coulomb plasticity. We extend here the essence of themethod to three-dimensional limit analysis with Von Mises plasticity. Although somethree-dimensional problems can be found in the literature [LZC04, VaG07, Va07], noneof them makes use of adaptive remeshing strategies..2 STATIC AND KINEMATIC PRINCIPLEWe are interested in finding the collapse load factor λ of a body Ω ∈ Rnsd (nsd = 2, 3is the number of space dimensions) subjected to homogeneous Dirichlet and Neumannboundary conditions on ∂Ω , i.e. u = 0 at Γu, and σn = λg at Γg, where λu and σ arethe displacement rate and stress fields, respectively, and Γu ∩ Γg = ∂Ω. In addition, theload per unit volume λf may be also considered in Ω. We further introduce the linearand bilinear forms:a(σ,u) =∫Ωσ : ε(u) dV,`(u) =∫Ωf · u dV +∫Γgu · g dΓ.where ε = 12(∇u+(∇u)T ) is the strain rate. The body is assumed to be a rigid-plasticmaterial, and therfore the stress field is subjected to belong to the set B of admissiblestresses, which in Von Mises plasticity, is given byB = {σ|devσ : devσ −23σ2y < 0}. (1)With this notation and definition at hand, the optimal load factor λ∗ is the solution ofthe following saddle-point problem:λ∗ = supλ,σ∈Binf`(u) = 1a(σ,u) = inf`(u) = 1supσ∈Ba(σ,u). (2)The first identity gives rise to the static principle of limit analysis, which states thatthe collapse load factor λ∗ corresponds to the maximum external potential that the bodycan sustain, while being in equilibrium. Indeed, the first identity may be also written as,λ∗ = supλ,σ∈Binfu(a(σ,u)− λ(`(u)− 1)) = supσ∈Binfa(σ,u) = λ`(u)λ= supλs.t.{a(σ,u) = λ`(u) ∀u ∈ Vσ ∈ B(3)The second identity in (2) gives rise to the kinematic principle, which defining theinternal rate of dissipation as D(u) = maxσ∈B∫Ωσ : ε(u), it states that λ∗ is the value of2J. Mun˜oz1∗, J. Bonet2, A. Huerta1 and J. Peraire3the minimum D(u) subjected to `(u) = 1. In summary, the optimum values of (λ,σ,u)satisfy the relation λ∗ = a(σ∗,u∗), and due to the saddle-point of structure of the op-timum, upper and lower bounds of the load factor, λLB and λUB, respectively, may beobtained obtained by satisfying just some of the supremum or infimum in (2), i.e.λLB = a(σ,u∗) ≤ λ∗ = a(σ∗,u∗) ≤ a(σ∗,u) = λUB. (4)3 DISCRETISATION AND BOUNDSWe will use discrete spaces ΣLB × VLB 3 (σLB,uLB) and ΣUB × VUB 3 (σUB,uUB)that preserve the first and second inequality in (4), respectively. The choices for suchspaces have been specified in Table 1, which are the same used in [Cir04, CPB]. By usingthem, the optimisation problem in (3) gives rise to the following two discrete problems(see [MBHP, CPB] for their derivation):λLB = sup λ λUB = supλs.t.{ALBσLB = fLBσLB ∈ BUBs.t.[AUB B]{ σUBtUB}= fUBσUB ∈ BtUB ∈ Bt(5)Bound Domain Space ConstrainUB Elements uLB ∈ P0 = VLB divu = 0σLB ∈ P1 = ΣLB σ ∈ BEdges uLB ∈ P1 = VLB Imposes edge equilibriumLB Elements uUB ∈ P1 = VUB divu = 0σUB ∈ P0 = ΣUB σ ∈ BEdges tUB ∈ P1 t ∈ BtTable 1: Spaces used in the upper and lower bound problem for the stress and displacement rate field.No continuity is enforced between elements or edges.The space Bt is defined as follows: Bt = {t| |tτ | ≤ σy/√3}, where tτ denotes thetangent component of the edge tension. In this way, we ensure that if the stresses at theboundary satisfy σ ∈ B, then t = σn ∈ Bt.We remark that the choices in Table 1 preserve the inequalities in (4) [Cir04, CPB,MBHP]. By using (ΣLB,VLB), we are satisfying the equilibrium equations of the con-tinuum, and therefore the minimisation is exactly achieved (over the restraint set ΣLB).On the other hand, resorting to spaces (ΣUB,VUB) we are able to obtain the maximum3J. Mun˜oz1∗, J. Bonet2, A. Huerta1 and J. Peraire3exactly (as a function of the restraint set VUB), or due to the use of a larger projectedspace Bt, obtain a higher value, and thus an upper bound to λ∗ [Cir04, MBHP].The expressions for matricesALB,AUB andBUB, and for vectors fLB and fUB have beenexplicitly given in [MBHP] for two-dimensional cases, but are easily extensible for threedimensions. We will here just comment and give details of the membership constraint ofthe stresses in (5), in order to turn the optimisation problem in (5) into the form of aSCOP.For 3D problems, and using the notationσT16 = (σ1 σ2 σ3 σ4 σ5 σ6) =( σxx σyy σzz σxy σxz σyz),the set B in (1) is expressed as,B = {σ|(σ1 − σ2)2 + (σ1 − σ3)2 + (σ2 − σ3)2 + 6(σ24 + σ25 + σ6)− 2σ2y < 0}.Furthermore, we will use the variables xT17 =(x1 x2 x3 x4 x5 x6 x7) instead of σ16. Theyare both related according to the following transformation:x27 = T−1σ16 , σ16 = Tx27 (6)withT−1 =0 −√6/2√6/2 0 0 0√(2) −√2/2√2/2 0 0 00 0 0√6 0 00 0 0 0√6 00 0 0 0 0√60 1 0 0 0 0,T =1/√6√2/2 0 0 0 10 0 0 0 0 12/√6 0 0 0 0 10 0 1/√6 0 0 00 0 0 1/√6 0 00 0 0 0 1/√6 0After using these expressions, it can be verified that the set B is equivalent to,B = {x17|x16 ∈ L6, x7 free, x1 =√2σy},where we have introduced the 6th dimensional Lorentz cone L6 = {x ∈ R6|x1 ≥√∑6i=2 x2i }. Similarly, in problem (5)b, the set Bt is modified by using the variableszT = (z1 z2 z3), and setting (z2 z3) = tTτ , which allows to define the admissible set for theboundary tensions as Bt = {z| z ∈ L3, z1 = σy/√3}.4 ERROR ESTIMATEThe successive refinement of the domain will yield tighter bound gaps ∆λ = λUB −λLB. However, the localisation of the displacement rate and tension field jumps thatcharacterises limit analysis, demands the design of effective adaptive remeshing strategies.In addition, it should be noted that for three-dimensional analysis, the necessarynumber of elements to obtain tight bounds is considerably higher than in two dimen-sions. More precisely, for a square domain (cube in 3D) with length L in each di-mension, and element size h, the number of elements is given by nelem2D = 2(L/h)2,4J. Mun˜oz1∗, J. Bonet2, A. Huerta1 and J. Peraire3while in 3D this is given by nelem3D = 6(L/h)3. In addition, the number of pri-mal variables (stresses) in 2D for the lower bound problem is given approximately bynvar2D = 3 × 3 × nelem2D = 18(L/h)2, while in 3D the number of primal nodal vari-ables is nvar3D = 4 × 6 × nelem3D = 144(L/h)3 = 8(L/h)nvar2D. Therfore, if uniformremeshing is used, and after 4 remeshing cycles, with similar final ratios L/h, the numberof primal variables in 3D analysis is at least 24 times higher than in 2D, which implies aconsiderable increase of the computational cost. For this reason, adaptive remeshing in3D analysis is a must. However, while adaptive strategies for 2D problems can be foundin [Cir04, CPB, LSKH05, MBHP], adaptive remeshing in 3D is lacking in the literature.It has been shown in [MBHP] that the total bound gap ∆λ may be expressed as thesum of elemental contributions ∆eλ and contributions of the edges ∆λξλ, i.e.,∆λ =nele∑e=1∆eλ +NI∑ξ=1∆ξλ. (7)The elemental and edge contributions are,respectively, given by,∆eλ = De(uUB)−(∫Ωe(−∇ · σLB) · uUB dV +∫∂Ωe(ne · σLB) · uUBdΓ)︸ ︷︷ ︸le(uUB)(8a)∆ξλ = Dξ(uUB) +∫ξe−e′σLBn · JuUBK dΓ (8b)where s = −JuUBK/|JuUBK|, and the expression of the elemental plastic dissipationrates De(uUB) and Dξ(uUB) are given, for Von Mises plasticity, by,De(uUB) =∫ΩeσY εeq dVDξ(uUB) =∫ξe−e′σY√3|JuUBK| dΓ.(9)We remark that due to the use linear edge tension fields, and the numerical integration,the dissipated energy at the edges Dξ(uUB) computed with the optimal values, differsfrom the analytical expression in (9). The reason for this discrepancy is illustrated inFigure 1 for a one dimensional edge and an hypothetical dissipation energy given byD(u) = max|t|≤ty∫ξtudΓ, with t and u scalar tension and displacement rate fields. Itcan be deduced that in 2D, and using 2 Gauss point quadrature (as it is done in ourtwo-dimensional case), the optimiser will choose the distribution of tensions in Figure bwhenever|ug1 + ug2| <ug2 − ug1√3.5J. Mun˜oz1∗, J. Bonet2, A. Huerta1 and J. Peraire3Indeed, the tension field is linear, and functions such as those given in Figure 1b are not.The optimiser will choose a linear function that maximises the integral in the dissipatedenergy, and thus will provide the optimal tension in Figure 1d. A similar reasoning appliesto the two-dimensional edges in 3D problems. Therefore, in order to compute correctly theedge contributions, we have considered the actual values of |tUBτ | and not the limit valueσY /√3. We note that in the latter case, a larger error contribution would be obtained,and therefore the equal sign in (7) would become a “≤”.−tytg2tg2(b)(d)(c)(a)tyug1ug2tu1u2−ty−tyu2u1tug2ug1−tytttyug1ug2u1u2u2u1ug2ug1tFigure 1: Comparison between analytical values (a,b) and numerical optimal values (c,d) obtained of thetension field t and displacement rates u for a one-dimensional case. Figures (a,c) correspond to a linearfield u with no sign variation, and Figures (b,d) a linear field u with a sign variation.5 UPPER BOUND ESTIMATEDue to the computational cost of the optimisation problems in (5), we present here amanner to obtain an upper bound without having to actually compute it. Instead of usinguUB in expressions (8), we will use the following averaged nodal displacements rates:u¯UBn =∑nee3n uUBe,nne(10)where ne is the number of elements connected to node n, and uUBe,n is the nodal dis-placement rate of node n in element e. The upper bound may be then computed asλ¯ = λLB + ∆λ(u¯UB), with ∆λ given in (7). Since we are using a non-optimal displace-ment rate, the minimisation of the min-max problem is not exactly satisfied, and thusa higher upper bound should be obtained. While this is true for plane stress analysis,where no constrain is solely imposed in the displacement rate field, in the other casesthe displacement rate u¯UB will in general not satisfy the constraint given in Table (1),namely divu = 0. Consequently, u¯UB is non-optimal and a non-feasible solution of the6J. Mun˜oz1∗, J. Bonet2, A. Huerta1 and J. Peraire3dual problem, and hence it cannot be used to compute the error contributions. Thenumerical experiments in the next section illustrate this fact.We also note that due to the reasoning described in the previous section, computingthe error with the analytical expressions in (9) yields elemental contributions that maybe higher than those computed as a SOCP, but in any case will not violate the strictnessof the upper the bound.6 RESULTS6.1 Vertical cut problem in 2DWe analyise the failure mechanism of a soil subjected to the gravity load and with avertical cut. The dimensions and the boundary conditions are indicated in Figure 2a. Thisproblem has been already analysed in [LS02a, LS02b, LSKH05, KHS05], and in particularin [MBHP] using the same formulation employed here and the same initial mesh shownin Figure 2. Nonetheless, we present it here to show the effects of averaging the lowerbound displacement rates to obtain a non-optimal upper bound solution.σn = στ = 0γu = 0u = 0u = 0(a) (b)Figure 2: Vertical cut problem geometry and boundary conditions (a), and initial mesh employed (b).The evolution of the bound gap for plane strain and plane stress is shown in Figure3, respectively. As mention in Section 5, in plane stress problems, when we average ateach node the displacement rates of the lower bound uLB, we are obtaining a non-optimalsolution, and thus a higher bound than the upper bound computed as a SOCP, λUB.However, in plane strain, the displacement rate given by u¯UB yield an estimated boundλ¯UB which may be higher or lower than λUB. The bounds in Figure 3 are in agreementwith this reasoning. However, it is worth pointing out that in both cases, the resultingerror estimate furnishes a lower bound which converges in a similar trend as the lowerbound computed with SOCP.7J. Mun˜oz1∗, J. Bonet2, A. Huerta1 and J. Peraire3 4.5 5 5.5 6 6.5 7 7.5 8 8.5 9 3  4  5  6  7  8  9λLog(nele)UB estimatedUB SOCP(a) 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3 3.1 3.2 3.3 3  3.5  4  4.5  5  5.5  6  6.5  7  7.5  8  8.5λLog(nele)UB estimatedUB SOCP(b)Figure 3: Evolution of the load factor when using computed upper bound with SOCP, and estimatedfrom the lower bound problem for the two-dimensional vertical cut problem in plane strain (a) and planestress (b).6.2 Vertical cut in problem in 3DThis a three-dimensional version of the previous problem, and therefore, similar sliplines should be expected. The problem has been run using an initial mesh with 55 ele-ments. Figure 5 shows the deformed mesh after 4 remeshing process, and the values of |tτ |at the internal edges. While the deformation pattern is very similar to the two-dimensionalcase, the evolution of the bounds in Figure 4 reveal a certain lack of convergence of thebounds. However, when using the estimated upper bound according to Section 5, anddespite the non-strictness of the upper bound, the lower bound converges. The source ofthis discrepancy in the two remeshing approaches is currently under investigation.7 CONCLUSIONSWe have extended the computation of bounds in two-dimensional limit analysis andVon Mises plasticity to three dimensions. We have also extended the error estimates andadaptive remeshing strategy.While the numerical model is theoretically able to provide such bounds, the numericaltests run so far have revealed the need to reduce the computational cost. More specifi-cally, with the SOCP solvers we have employed [Stu99, TTT03], the bottle neck of theoptimisation problem is not as much the time required, but the memory needed to solvepractical problems.In addition, alternative remeshing strategies are being studied to improve the conver-gence of the bound gap when the upper and lower problems are being solved. researcheffort is being applied to the reduction paid8J. Mun˜oz1∗, J. Bonet2, A. Huerta1 and J. Peraire3 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6 4  5  6  7  8  9  10λLog(nele)LB (UB estimated)UB (UB estimated)LB (SOCP)UB (SOCP)Figure 4: Evolution of the load factor when using computed upper bound with SOCP, and estimatedfrom the lower bound problem for the three-dimensional vertical cut problem.(a) (b)Figure 5: Effective stresses on the deformed mesh of the upper bound problem (a) and modulus of thetangent edge tensions (b) when using 4316 elements.REFERENCES[Cir04] H Ciria. Computation of Upper and Lower Bounds in Limit State Analysisusing Second-Order Cone Programming and Mesh Adaptivity. PhD thesis,Dep. Aeron. and Astron, MIT, USA, 2004.[CPB] H Ciria, J Peraire, and J Bonet. Mesh adaptive computation of upper and9J. Mun˜oz1∗, J. Bonet2, A. Huerta1 and J. Peraire3lower bounds in limit analysis. Int. J. Num. Meth. Engng.. Accepted.[KHS05] K Krabbenhoft, A V Lyamin M Hjiaj, and S W Sloan. A new discontinuousupper bound limit analysis formulation. Int. J. Num. Meth. Engng., 63:1069–1088, 2005.[LS02a] A V Lyamin and S W Sloan. Lower bound limit analysis using non-linearprogramming. Int. J. Num. Meth. Engng., 55:576–611, 2002.[LS02b] A V Lyamin and S W Sloan. Upper bound limit analysis using linear finiteelements and non-linear programming. Int. J. Num. Anal. Meth. Geomech.,26:181–216, 2002.[LSKH05] A V Lyamin, S W Sloan, C Krabbenhoft, and M Hjiaj. Lower bound limitanalysis with adaptive remeshing. Int. J. Num. Meth. Engng., 63:1961–1974,2005.[LZC04] Y Liu, X Zhang, and Z Cen. Numerical determination of limit loads forthree-dimensional structures using boundary element method. Eur. J. Mech.A/Solids, 23:129–138, 2004.[MBHP] J Mun˜oz, J Bonet, A Huerta, and J Peraire. Upper and lower bounds in limitanalysis: adaptive meshing strategies and discontinuous loading. Submitted.[Stu99] J F Sturm. Using SeDuMi 1.02, a MATLAB toolbox for optimization oversymmetric cones over symmetric cones. Optim. Meth. Soft., 11-12:625–653,1999. Avail. http://sedumi.mcmaster.ca.[TTT03] RH Tu¨tu¨ncu¨, KC Toh, and MJ Todd. Solving semidefinite-quadratic-linearprograms using SDPT3. Mathem. Progr. Ser. B, 95:189–217, 2003. Avail.http://www.math.nus.edu.sg/ mattohkc/sdpt3.html.[Va07] M Vicente da Silva and A Ant ao. A non-linear programming method approachfor upper bound limit analysis. Int. J. Num. Meth. Engng., 72:1192–1218, 2007.[VaG07] M Vicente da Silva, A Ant ao, and N Guerra. Numerical implementation ofthe kinematical theorem: Application to the determination of 3d passive earthpressures of cohesionless. In COMPLAS2007, Fundamentals and applications,Barcelona, Spain, 5-7 September 2007. ECCOMAS, Owen et al (eds).Con el apoyo de Universitat Polite`cnica de Catalunya y E.T.S. d’Enginyers de Camins, Canals i Portsde Barcelona10

http://creativecommons.org/licenses/by-nc-nd/3.0/es/

Bounds and adaptivity for 3D limit analysis

Abstract

Similar works

Full text

Available Versions

UPCommons

UPCommons. Portal del coneixement obert de la UPC