Parallel 3D limit analysis via the alternating direction method of multipliers

Parallel Computing TechniquesParallel 3D limit analysis via the alternating direction method of multipliers XIII International Conference on Computational Plasticity. Fundamentals and Applications COMPLAS XIII E. O˜nate, D.R.J. Owen, D. Peric & M. Chiumenti (Eds) PARALLEL 3D LIMIT ANALYSIS VIA THE ALTERNATING DIRECTION METHOD OF MULTIPLIERS NUNO M. DEUSDADO∗, M ´ ARIO V. DA SILVA AND ARMANDO N. ANT ˜ AO UNIC, Department of Civil Engineering Faculty of Science and Technology (FCT) Nova University of Lisbon (UNL) Quinta da Torre, 2829-516 Monte da Caparica, Portugal ∗e-mail: [email protected] Key words: Slope, Limit Analysis, Upper Bound, Lower Bound, Finite Element Method, Parallel Processing Abstract. This document presents a parallelized implementation of a FEM formulation concerning both the upper and lower bound limit analysis theorems. By applying these methods to evaluate the stability of a vertical cut (i.e. vertical slope) in drained and undrained conditions, it is possible to discuss the efficiency of this strategy. 1 INTRODUCTION The prediction of the maximum load that can be supported by a structure, known as collapse load, is a fundamental problem in civil engineering. For the analysis of mechanical structures subject to noncyclic loading and presenting a perfect plastic behaviour, the limit analysis theorems have proven to give good results and to be an efficient and competitive strategy comparing with other methods. The static and kinematical limit analysis theorems allow the computation of both the lower and upper bounds for the collapse loads of these mechanical structures. Nowadays, and in spite of the remarkable evolution of computers performance, the determination of accurate collapse load estimates can still represent a significant computational effort. In fact, when strict high quality bounds solutions are searched, it is required a high degree refinement of the mesh, and thus substantial CPU and RAM resources, particularly for 3D problems. The alternating direction method of multipliers (ADMM) technique has been used by the authors to solve these problems, due to its iterative solution scheme based on an operator splitting algorithm, which is not only easy to implement but also suitable to efficiently solve large-scale variational problems with parallel processing. In this work, after a brief description of these numerical formulations, the results of its application to evaluating the stability of vertical cut are presented and some conclusions are drawn. 1077 Nuno M. Deusdado, M´ario V. da Silva and Armando N. Ant˜ao 2 NUMERICAL FORMULATION As it is usual in limit analysis, an optimization problem is formulated. In this paper, two complementary finite element analysis problems based on the kinematical and the static theorems are presented in order to determine an upper and a lower bound of the true collapse load multiplier of a mechanical structure, respectively. These problems have a significant number of decision variables and constraints. To obtain its solutions the ADMM algorithm has been used by the authors, because it is a very versatile and robust (practically never fails to converge) technique and is also inherently parallelizable [1]. This algorithm has an iterative solution scheme based on an operator splitting algorithm leading to two subproblems. The first results in a quadratic unconstrained problem, easily solvable through a linear system of equations. Conveniently, the governing system matrix remains unaltered during the entire iterative process. In fact, all the modifications appearing during the iterative process affect only the right hand side term. Therefore, a single matrix factorization procedure is needed throughout the whole iterative process. The second subproblem collects all the nonlinear parts of the original problem and consists in a sum of small independent optimization problems. In most cases of practical interest, using nonsmooth techniques to deal with not everywhere differentiable function, it is possible to derive closed-form solutions. Otherwise, special-purposed algorithms can be applied efficiently. In summary, the first subproblem (global minimization) takes into account the contribution of all finite elements simultaneously, while the second subproblem (local minimization) is carried out at element level, meaning that it is solved independently for each element [2]. 2.1 Kinematical limit analysis theorem Based on the kinematical theorem, the search for an upper bound (UB) of the collapse load multiplier, λ(λ∈R+), of a mechanical structure can be formulated as the following mathematical minimization problem: minimize λ(˙u, ˙ε)=Ω D(˙ε)dΩ− Π(˙u) subject to u=0 in Γu ˙ε=B˙uin Ω ˙ε∈Cc Πλ(˙u)=Ω bT λ˙udΩ+Γσ tT λ˙udΓσ=1 (1) where udenotes the displacement field, ˙uthe velocity field, ˙εthe plastic strain rate, Bthe differential compatibility operator and Dthe plastic energy dissipation rate per unit volume (expressed in terms of kinematic fields only [3]). The structure is submitted to a given distribution of a set of constant loads: body loads, bon the domain Ωand external surface forces, ton the static boundary Γσ. These loads are divided in constant loads (defined with the overline “∼” symbol) and the loads affected by the load multiplier 1078 Nuno M. Deusdado, M´ario V. da Silva and Armando N. Ant˜ao (defined with the subscript “λ” symbol). The first two equations from the problem (Eq. 1) establish, respectively, the compatibility conditions in the kinematic boundary, Γu, and on the domain. The next equation, enforces the normality plastic flow rule by defining Cc as the space formed by all plastically admissible strain states (orthogonal to the yielding surface in at least one point [3]). Finally, the last constraint scales the external work rate for the variable loads, Πλ(˙u). For each finite element i(of the total number of finite elements, nE), two independent and simultaneous approximations are adopted for the velocity field and for the plastic strain rate (Eq .2). For the velocity field approximation (d(i)), quadratic nodal shape functions (Q(i)) are adopted, using 10-node tetrahedral elements. For the plastic strain field approximation (e(i)), linear nodal shape functions (L(i)) using 4-node tetrahedral elements are adopted [4]. ˙u(i)=Q(i)d(i);˙ε(i)=L(i)e(i)(2) The ADMM algorithm [1] is used to find the optimal solution of the problem expressed by equation 1. This scheme, as shown in Table 1, is a cyclic kiteration with a two steps minimization procedure followed by a third stage where the Lagrange multiplier, µ, is updated. In Table 1, µλrepresent the Lagrange multiplier (of the equation that scales the external work rate), rdenotes the penalty parameter and A0is an diagonal matrix where the entries are equal to 1 except those affecting the shear strain components with a value equal to 1 2. The assemblage procedure of A,f(ek,µ k) and Fλ(nodal force vector) follows the procedure normally applied to the finite element method (FEM), f(ek,µ k)=  F+ nE  i=1 Ω(i) (BQ(i))TA0(re(i) k−µ(i) k)dΩ A= nE  i=1 Ω(i) (BQ(i))TA0(BQ(i))dΩ Fλ= nE  i=1 Ω(i) (Q(i))Tb(i) λdΩ+Γ(i) σ (Q(i))Tt(i) λdΓσ (3) The local minimization process for each point jof the total number of points of the finite element mesh, np, starts by computing the principal axes of s(j)(that are coincident with the principal axes of the tensor e(j)[5]) using an eigenvalue decomposition. In the principal space the stationary point of the objective function without constraints is computed. If the solution belongs to the feasible domain (if it complies with the nonlinear constraint), it matches the optimal solution, otherwise it is necessary to project the stationary solution onto the feasible domain. In [6, 7] it can be found the description in detail of this projection and the assemblage process of all the operators presented in table 1. 1079 Nuno M. Deusdado, M´ario V. da Silva and Armando N. Ant˜ao Table 1: ADMM solution scheme for kinematical theorem Global minimization minimize LG(dk,µ λk)=µλk(1 −FT λdk)−f(ek,µ k)Tdk+r 2dT kAdk solution: rA Fλ FT λ0dk µλk=f(ek,µ k) −1 Local minimization minimize LL(e(j) k)=D(e(j) k)−s(j) k TA0e(j) k+r 2e(j) k TA0e(j) k subject to s(j) k=µ(j) k+rBQ(j)d(j) k e(j) k∈Cc Update Lagrange multiplier µ(j) k=µ(j) k−1+r(BQ(j)d(j) k−e(j) k) 2.2 Static limit analysis theorem Based on the static theorem, the search for a lower bound (LB) of the collapse load multiplier, can also be formulated as a mathematical maximization problem with three constraints. Two equations establish the equilibrium conditions on the domain and in the static boundary and one inequality enforces the yield condition. Like in the kinematic theorem, for each finite element itwo independent and simultaneous approximations are adopted. The ADMM algorithm [1] finds the optimal solution of this problem through a cyclic kiteration scheme defined with a two step minimization procedure followed by a third stage where the Lagrange multiplier is updated. The local minimization process follows the same scheme presented earlier, but in this case the minimization is performed in order to the stresses and the feasible domain is defined by f(σ)≤0. 2.3 Parallelization technique The parallelization of the solution schemes (table 1) have been implemented as follows: - Local minimization and update lagrange multiplier: both these stages of the algorithm are entirely suitable for parallelization, as mentioned previously. In short, these processes are performed independently for each element, meaning that each 1080 1081 Nuno M. Deusdado, M´ario V. da Silva and Armando N. Ant˜ao 3.5 4 4.5 5 5.5 6 6.5 7 1 1.5 2 2.5 3 3.5 4 4.5 5 γH cu L H UB [8] UB [9] LB [9] UB [10] UB present work LB present work Figure 2: Undrained vertical cut problem: Comparison with previous results visible) and symmetric faces. For the lower bound analsysis, the applied stress (normal σn and shear τ) boundary condition are given. The figure 1 presents a coarse finite element partition. The volume is divided in an orthogonal grid where each division is discretized by 24 tetrahedrons. In the present work, tetrahedra with ten nodes define the mesh for the kinematic theorem and for the static theorem the mesh is composed by tetrahedra with four nodes. 4 RESULTS The calculations presented in this work were performed in a cluster of 48 nodes, with a maximum 4Gb of RAM per node. The available resources were used at their maximum capacity in all the calculations. This implies that the size of the mesh elements depends on the dimensions of the model treated, meaning that, as the value of Lgrows, the meshes get coarser. The discretization is similar for the UB and LB cases, the difference residing in the fact that the UB calculation use ten nodes tetrahedra and the LB use four node tetrahedra. The results obtained are presented in figures 2 and 3. In these figures are also presented the results from UB and LB approaches of other authors ([9] and [10]) and the results issuing from a UB approach obtained previously by some of the authors of the present work [8]. As referred previously, two types of calculations, as the behaviour of the soil is concerned, were performed: undrained conditions and drained conditions. The stability of the undrained vertical cut problem is represented in figure 2. As expected the UB obtained are very close to the previously obtained by [8]. They represent a clear improvement regarding the UB results of [9] and [10]. Concerning the results of the LB approach there is also an improvement of the best result known in the literature. However this improve1082 Nuno M. Deusdado, M´ario V. da Silva and Armando N. Ant˜ao 6.5 7 7.5 8 8.5 9 9.5 10 10.5 11 11.5 1 1.5 2 2.5 3 3.5 4 4.5 5 γH c′ L H UB [10] UB present work LB present work Figure 3: Drained vertical cut problem with φ′= 30◦: Comparison with previous results ment is, in this case, a small one. The results obtained allow to say that, for practical purpose from an engineering point of view, the stability number is known. In fact, the error of the calculations is at most equal to 2% which encloses the collapse solution with great accuracy. The calculations for the drained conditions consider a friction angle of 30◦and a c′ equal to 1 KPa. The results obtained are presented in figure 3. It can be observed that the UB result of the present work improves considerably the best result known in the literature, except for L H= 5. As in the undrained case, the UB and LB results are almost coincident, but in this case the variation is less than 4% which also correspond to a good accuracy. In both cases (undrained and drained conditions) the stability number γH cwill reduce with increasing L/H ratio. For the undrained problem, the meshes for L/H=1 and L/H=5 cases are shown in figures 4 and 5, respectively. For the drained case analyzed, the meshes for L/H=1 and L/H=5 cases are shown in figures 6 and 7, respectively. The meshes resulting from the UB computations are represented in its deformed shape and include the plastic dissipation pattern. For the static theorem, the figures present the yield function pattern. From these figures it can be seen that: - The deformed mesh shows that the mechanism is formed by three zones: The first zone the soil rest in place with no dissipation, the second zone is definedbya high plastic dissipation (shear surface), and finally the third zone the soil has no dissipation but move out of the soil mass like a rigid body; - The shear surface (failure surface) has the shape of a curvilinear cone [11]; 1083 Nuno M. Deusdado, M´ario V. da Silva and Armando N. Ant˜ao (a) Deformed mesh and plastic dissipation (UB) (b) Yield function (LB) Figure 4: Undrained vertical cut problem with L/H=1 (a) Deformed mesh and plastic dissipation (UB) (b) Yield function (LB) Figure 5: Undrained vertical cut problem with L/H=5 - For the same L/H, the plastic dissipation (UB) and the yield function (LB) are very similar; - As expected, the UB meshes results of the purely cohesive soil (figures 4a and 5a) coincide with the graphical results reported previously by our team [8]; - A slope with frictional soil have a more vertical sliding surface than a slope with cohesive soil, so the mobilized mechanism is smaller. Collecting both information, it can be concluded that the appearance of the slipe surface (second zone of the mechanism) influence significantly the stability number. With the increasing of the L/H ratio the influence of the proximity of the end surface decreases, so the contribution of the two curved ends relative to the straight slip surface will decrease, leading to more unstable slopes. 1084 Nuno M. Deusdado, M´ario V. da Silva and Armando N. Ant˜ao (a) Deformed mesh and plastic dissipation (UB) (b) Yield function (LB) Figure 6: Drained vertical cut problem with L/H=1 (a) Deformed mesh and plastic dissipation (UB) (b) Yield function (LB) Figure 7: Drained vertical cut problem with L/H=5 5 CONCLUSIONS In this document, the fundamentals of the numerical parallel implementation of a FEM formulation concerning both the upper and lower bound limit analysis theorems were briefly presented. The vertical slope in drained (with φ′= 30◦) and undrained condition examples were analyzed. Due to the small gap between the upper and lower bound limits a good estimate of the exact solution can be inferred. The graphical representation for both theorems also shows a great agreement between the plasticity zone obtained by the two approaches. These results confirm the effectiveness of the parallel limit analysis strategy advocated in this work and encourage its use in the study of further geomechanical problems. 1085