Speedup of Elastic–Plastic Analysis of Large-scale Model with Crack Using Partitioned Coup ling Method with Subcycling Technique

Yasunori Yusa and Shinobu Yoshimura

1　Introduction

Fracture mechanics plays a significant role in evaluating structural integrity with maintenance codes such as ASME Boiler and Pressure Vessel Code Section XI(Rules for Inservice Inspection of Nuclear Power Plant Components)and the JSME Code for Nuclear Power Generation Facilities(Rules on Fitness-for-service for Nuclear Power Plants).Although linear elastic fracture mechanics(LEFM)is utilized for fatigue and stress corrosion cracking(SCC)in such codes,elastic–plastic cracks are not generally considered at present.Recently,mainly due to disasters such as earthquakes and tsunamis,it has become important to study the behavior of a cracked structure in a plastic state.There is an issue in analyzing a cracked structure with plasticity.Elastic–plastic analysis with a nonlinear finite element method(FEM)requires much more computation time than linear elastic analysis.Plasticity is,however,observed near the crack in the case of small strain,whereas a large portion far from the crack can be regarded as an elastic body.A structure of complicated shape in the real world tends to be modeled as a largescale finite element model.For this issue,we have studied a partitioned coupling method[Yusa and Yoshimura(2013)],which was imported from the field of fluid–structure interaction coupling[M inami and Yoshimura(2010)].In this method,the entire analysis model is decomposed into two non-overlapped domains,i.e.,global and local domains.The two domains are analyzed separately and iteratively with an iterative method to satisfy both geometrical compatibility and force equilibrium on the global–local interface.

Various methods that utilize two meshes or a mesh with a non-numerical solution have been proposed.S-version FEM(SFEM)[Fish(1992)]was proposed to reduce the number of finite elements required for a sufficiently accurate solution and to reduce the human effort in meshing.In SFEM,a fine local mesh that may contain a crack is superimposed on a coarse global mesh,and the two meshes are connected with Lagrange multipliers.Elastic–plastic problems were analyzed explicitly using SFEM[Nakasumi,Suzuki,Fujii,and Ohtsubo(2002);Okada,Endoh,and Kikuchi(2007)].Nakasumi,Suzuki,Fujii,and Ohtsubo(2002)analytically gave that the elastic constants on the local mesh can be different from those on the global mesh.Okada,Endoh,and Kikuchi(2007)applied the SFEM for an elastic–plastic fracture problem and evaluated theJ-integral.In their studies,incremental analysis was employed without an implicit iterative method such as the New ton–Raphson method.Nodal unbalance evaluation,which is necessary for implicit iterative methods,is troublesome in SFEM since the effect of the local reaction force is complicatedly added to that of the global reaction force.Suzuki,Ohtsubo,Nakasumi,and Shinmura(2002)utilized iterative methods with SFEM and analyzed benchmark problems with a general-purpose commercial code.The iterative substructure method[Nishikawa,Serizawa,and Murakawa(2007)]is a similar method to SFEM with regard to meshing.In this method,an iterative method is used to obtain a converged solution.Thermal elastic–plastic problems in welding were analyzed by Nishikawa,Serizawa,and Murakawa(2007).The global–local method[Whitcomb(1991)],which is also known as the zooming method,is a popular approach in structural analysis.In this method,the numerical result of a global analysis is given to a local model as the boundary conditions,and then,the local model is analyzed in detail.Since the result of the local analysis cannot affect the global analysis,in general,either geometrical compatibility or force equilibriumis not satisfied.Iterative global–local analysis was studied by Whitcomb(1991)to solve this issue.These three iterative methods[Suzuki,Ohtsubo,Nakasumi,and Shinmura(2002);Nishikawa,Serizawa,and Murakawa(2007);Whitcomb(1991)]described above are not sophisticated with regard to it-erative algorithms.In the present study,a quasi-New ton method,which showed the best performance of several iterative algorithms with several line search algorithms in fluid–structure interaction benchmarks[M inami and Yoshimura(2010)],is used.In addition,in the case of elastic–plastic problems or generally nonlinear problems,methods based on the principle of superposition would not be appropriate since the principle does not hold.Extended FEM(XFEM)[Mo?s,Dolbow,and Belytschko(1999)]was proposed to avoid meshing efforts concerning a crack.In the XFEM,a Heaviside function representing the crack discontinuity is superimposed on the finite element interpolation functions.The XFEM has been applied for various problems[Abdelaziz and Hamouine(2008)],which include elastic–plastic fracture problems.Elguedj,Gravouil,and Combescure(2006)employed the Hutchinson–Rice–Rosengren(HRR) field[Hutchinson(1968);Rice and Rosengren(1968)]to represent the singularity of an elastic–plastic crack.Prabel,Combescure,Gravouil,and Marie(2007)solved a dynamic crack propagation problemin elastic–plastic media using XFEM linear function approximation.Samaniego and Belytschko(2005)modeled shear bands in plasticity using XFEM.The elastic finite element alternating method(FEAM)[Nishioka and Atluri(1983)]was proposed to connect the analytical solutions of the stress intensity factors(SIFs)with the uncracked structure mesh.In FEAM,the analytical solution and the mesh are solved iteratively,and a converged solution is obtained.The elastic FEAM was then extended to elastic–plastic FEAM[Pyo,Okada,and Atluri(1995)]by coupling the FEAM with the initial stress algorithm[Nikishkov and Atluri(1994)].Compared to XFEM and FEAM,our method can be applied for any plasticity models as far as the model can be applied for conventional FEM.FEAM was extended to SGBEM–FEM alternating method that utilizes symmetric Galerkin boundary element method(SGBEM)instead of a theoretical solution[Nikishkov,Park,and Atluri(2001)].In SGBEM–FEM,a global uncracked FEM mesh and a local cracked SGBEM mesh are solved alternately to satisfy force equilibrium,using forces on a mesh interface and for cesonacrack surface.Various two-and three-dimensional problems of fracture and fatigue crack grow th were analyzed by SGBEM–FEM[Dong and Atluri(2013a,b)].The domain decomposition method(DDM)based on the preconditioned conjugate gradient(PCG)methods[Yoshimura,Shioya,Noguchi,and Miyamura(2002);Ogino,Shioya,Kawai,and Yoshimura(2005);Miyamura,Noguchi, Shioya,Yoshimura,and Yagawa(2002)]has been used for parallel finite element analysis(FEA).In DDM,the entire analysis model is decomposed into multiple subdomains,and the PCG solver is employed for the subdomain interface.Since multiple FEAs are performed for every subdomain at every PCG iteration,this approach gives high parallelism.The ADVENTURE System by Yoshimura,Shioya,Noguchi,and Miyamura(2002)includes a parallel solid mechanics solver based on DDM.Ogino,Shioya,Kawai,and Yoshimura(2005)solved an elastodynamic problem of a large-scale structure using the solver on a supercomputer.Miyamura,Noguchi,Shioya,Yoshimura,and Yagawa(2002)applied DDM for a large-scale elastic–plastic problem.Plasticity is observed everywhere in their study,whereas it is modeled to be observed locally near the crack in the present study.The latter approach can lead to an efficient simulation since a large portion far from the crack is modeled as an elastic body.

In this paper,the partitioned coupling method,which has already been applied for a linear elastic fracture problem[Yusa and Yoshimura(2013)],is applied for an elastic–plastic problem with a crack.The cracked local domain is modeled as an elastic–plastic body,whereas the large-scale global domain is modeled as an elastic body.The subcycling technique[Farhat,Lesoinne,and Maman(1995)]is utilized for incremental analysis to reduce the number of global elastic analyses.The methodology is presented in detail in Section 2.A benchmark problemis then analyzed using the present method and conventional FEM in Section 3,and the computational performance of these methods is carefully examined.A conclusion is finally given in Section 4.

2 Method

2.1　Partitioned Coupling Method

First of all,the entire analysis model is decomposed into two non-overlapped domains,as shown in Fig.1.Since the two domains are analyzed separately in a partitioned coupling method,in general,either geometrical compatibility or force equilibrium on the global–local interface is not satisfied.However,they become satisfied by an iterative solution technique.Here,the analysis on the local domain,?L,is represented as

and the analysis on the global domain,?G,is represented as

Here,kis the iteration step,uΓis the displacement vector on the global–local interface,Γ,?uΓis the predicted displacement vector onΓ,andfΓis the force vector onΓ.Lis a function in which local analysis is performed under an enforced displacement boundary condition,uΓ,onΓ,and the negative of the reaction force vector,fΓ,is returned.Gis also a function in which global analysis is performed under an external force boundary condition,fΓ,onΓ,and the displacement vector,uΓ,is returned.These functions themselves would be nonlinear ow ing to plasticity.From Eq.1 and Eq.2,

is obtained.A residual vector,R,is now de fined as a function of

When geometrical compatibility on the global–local interface is satis fied,the residual should vanish as

This nonlinear equation is to be solved by an iterative solution algorithm described in the next paragraph.Equation 5 is checked numerically with a tolerance,τ,as

in the present study.In addition,force equilibrium also becomes satisfied when geometrical compatibility is satisfied at thek-th iteration step as

From this equation,

is obtained.The function ofLis applied to both sides,and then Eq.1 is used.Finally,

is derived.The left side of this equation is the force-based residual vector,whose formis similar to Eq.4.

In the present study,the limited-memory Broyden method[Minami and Yoshimura(2010);Kelley(2003)]is adopted to solve Eq.5.The limited-memory Broyden method is a quasi-New ton method and satisfies the secant condition

Figure 1:Decomposed analysis model and assumed boundary conditions on global–local interface.

Here,Bis an approximate Jacobian matrix,

is the search direction vector,and

is the residual vector.The Broyden updating formula is here de fined as

Using the Sherman–Morrison formula,

whereMis an invertible matrix,anduandvare vectors whose size is the same as the matrix size,one can derive

and

A dense matrix,B,is eliminated in the above equations so that this method does not require a large amount of memory.The algorithm of the limited-memory Broyden method is summarized as follows.

end while

An initial inverse approximate Jacobian matrix,B(0)-1,is a user-defined parameter,which is determined to be a diagonal matrix whose diagonal entries are 0.1 in the present study.

2.2　Subcycling Technique

Owing to strain path dependence,elastic–plastic analysis with nonlinear FEM generally requires incremental analysis as follows:

forloop of incremental stepsdo

Analysis with New ton–Raphson method

end for

One can select two approaches for incremental analysis with the partitioned coupling method.

The first approach is the incremental approach,in which partitioned coupling iterations are conducted at every incremental step.The algorithm of the incremental approach is described as follows.

forloop of incremental stepsdo

whileloop of partitioned coupling iterationsdo

Local analysis with New ton–Raphson method

Global analysis

end while

end for

An initial guess,,of the Broyden method described in the previous subsection is zero filled at the first incremental step or filled with a previous converged solution,at other incremental steps.However,in this approach,the global domain is analyzed with incremental steps,even though it is an elastic body.

The second approach is the subcycling approach,in which incremental steps are conducted at each partitioned coupling iteration.This approach enables the reduction of the number of global analyses.The subcycling technique is occasionally used in the field of fluid–structure interaction coupling[Farhat,Lesoinne,and Maman(1995)],especially with staggered algorithms.This technique is preferred because the time scale required for fluid analysis is smaller than that for structural analysis.This situation is very similar to the present study.The algorithm of the subcycling approach is described in the following:

forloop of partitioned coupling iterationsdo

Decision of number of incremental steps

whileloop of incremental stepsdo

Local analysis with New ton–Raphson method

end while

Global analysis

Turning back incremental time(restoring internal variables such as equivalent plastic strain and yield stress)

end for

In this algorithm,the number of incremental steps is determined at each partitioned coupling iteration by a simple formula:

Figure 2:Dimension parameters and boundary conditions of pressure vessel model.

is calculated as and ?εcharis a user-defined parameter.uΓidenotes an assumed enforced displacement boundary condition on the global–local interface,andxirepresents the nodal coordinates.?εcharis determined to be 0.01%in the present study.

3 Benchmark

A cracked pressure vessel model with 6 million degrees offreedom(DOFs)was analyzed by the developed partitioned coupling solver as well as the developed conventional FEM solver.The dimension parameters and boundary conditions of the model are shown in Fig.2.An 18-deg through-wall crack is introduced at the nozzlepart.The decomposed mesh is visualized in Fig.3 by using ahandy graphics and GUI library,AutoGL[Kawai(2006)].The numbers of elements,nodes,and DOFs in the global mesh are 1,308,720,2,117,000,and 6,351,000,respectively.

Those in the local mesh are 54,710,80,066,and 240,198,respectively.The ratio of the number of global DOFs to the number of local DOFs is 26:1.The numbers of nodes and DOFs in the global–local interface are 960 and 2,880,respectively.The employed material parameters are a Young’s modulus of 210 GPa,a Poisson’s ratio of 0.3,and an initial yield stress,σy0,of 250 MPa.The employed stress–strain curve is

The von Mises’equivalent stress distribution computed by the subcycling partitioned coupling solver is visualized in Fig.4.Its deformation is magnified by 100.The computed yielding zone is visualized in Fig.5.The result of the partitioned coupling analysis seems almost the same as that of the conventional FEA.

Figure 3:Mesh of pressure vessel model.Entire mesh(top left),mesh near skirt(bottom left),sectioned mesh near skirt(top middle),sectioned mesh near nozzle(bottom middle),decomposed mesh near nozzle(top right),and local mesh(bottom right).

Figure 4:Computed equivalent stress distribution of pressure vessel model.Entire model(top left),model near skirt(bottom left),sectioned model near skirt(top middle),sectioned model near nozzle(bottom middle),decomposed model near nozzle(top right),and local model(bottom right).

The computational performance was measured on a personal computer with an Intel Core i7-3930K Sandy Bridge CPU,64 GB of DDR3 SDRAM PC3-12800,and a Debian GNU/Linux 6.0 squeeze operating system.Intel C Compiler 13.0 was used with the-fast flag,and the PARDISO linear system solver in Intel Math K-ernel Library(MKL)10.2 was used for matrix LDL factorization and the triangular solution,which is also called forward and backward substitutions.The incremental partitioned coupling analysis was 3.40 times faster than conventional FEA,and the subcycling analysis was 3.34 times faster.Detailed computation time and memory usage are described in Tab.4.Speedup of the incremental partitioned coupling method was caused by a constant stiffness matrix in the global domain.The large stiffness matrix in the global elastic domain remained constant throughout the analysis.The matrix was factorized only once,and then,multiple triangular solutions were conducted.On the other hand,speedup of the subcycling method was caused by a reduction in the number of global analyses. Eighteen global analyses were conducted in the present method,whereas conventional FEA required 47.If a linear system solver based on the preconditioned conjugate gradient method,which has been used for very-large-scale problems[Yoshimura,Shioya,Noguchi,and Miya-mura(2002);Ogino,Shioya,Kawai,and Yoshimura(2005);M iyamura,Noguchi,Shioya,Yoshimura,and Yagawa(2002)],is selected for the global analysis,the subcycling partitioned coupling solver would still remain fast,but the incremental solver would become slow.This is because the constant matrix does not contribute to speedup in the PCG methods.

Table 1:Given and measured number of partitioned coupling iterations and incremental steps in pressure vessel model analysis

Table 2:Measured number of New ton–Raphson iterations in pressure vessel model analysis

Table 3:Measured number of linear system solutions in pressure vessel model analysis

4　Conclusion

In this paper,the partitioned coupling method was applied for elastic–plastic analysis of a large-scale model with a crack.In the method,the entire analysis model is decomposed into two non-overlapped domains(i.e.,global and local domains),and the two domains are analyzed with an iterative method.The cracked local domain was modeled as an elastic–plastic body,whereas the large-scale global domain wasmodeled as an elastic body.The subcycling technique was utilized for incremental analysis to reduce the number of global elastic analyses.In Section 2,the modeling of the partitioned coupling method with the limited-memory Broyden method and the subcycling technique with a simple formula to determine the number of incremental steps were presented in detail.A cracked pressure vessel model with 6 million degrees offreedom was analyzed in Section 3.The number of partitioned coupling iterations,incremental steps,New ton–Raphson iterations,and linear system solutions were carefully examined,and the computational performance was also measured.The incremental partitioned coupling solver was 3.40 times faster than the conventional FEM solver,and the subcycling partitioned coupling solver was 3.34 times faster.This is because the number of linear system solutions was reduced by the subcycling technique,and the large stiffness matrix in the global domain remained constant throughout the analysis.In the partitioned coupling method with the subcycling technique,this feature would remain if a linear system solver based on a preconditioned conjugate gradient method were selected instead of the direct LDL solver used in the present study.

Table 4:Measured computation time and memory usage in pressure vessel model analysis

In the future,modeling of more nonlinear and complicated phenomena such as large strain,crack surface contact,and crack propagation will become possible in the local domain,as will effective analysis with the partitioned coupling method,in which two different solvers work together.

Acknowledgement:This work was supported by a MEXT grant forHPCI S-trategic Program Field 4:Next-generation Industrial Innovations.The authors wish to thank Dr.M inamiat Chuo University for advising on constructing the Broyden algorithm and the members ofADVENTURE Projectfor engaging in helpful discussions.

Abdelaziz,Y.;Hamouine,A.(2008):A survey of the extended finite element.Computers and Structures,vol.86,no.11–12,pp.1141–1151.

Dong,L.;Atluri,S.N.(2013):Fracture&fatigue analyses:SGBEM–FEM or XFEM?part 1:2D structures.Computer Modeling in Engineering and Sciences,vol.90,no.2,pp.91–146.

Dong,L.;Atluri,S.N.(2013):Fracture&fatigue analyses:SGBEM–FEM or XFEM?part 2:3D solids.Computer Modeling in Engineering and Sciences,vol.90,no.5,pp.379–413.

Elgued j,T.;Gravouil,A.;Combescure,A.(2006):Appropriate extended functions for X-FEM simulation of plastic fracture mechanics.Computer Methods in Applied Mechanics and Engineering,vol.195,no.7–8,pp.501–515.

Farhat,C.;Lesoinne,M.;Maman,N.(1995):Mixed explicit/implicit time integration of coupled aeroelastic problems:Three- field formulation,geometric conservation and distributed solution.International Journal for Numerical Methods in Fluids,vol.21,no.10,pp.807–835.

Fish,J.(1992):The s-version of the finite element method.Computers and Structures,vol.43,no.3,pp.539–547.

Hutchinson,J.W.(1968):Singular behaviour at the end of a tensile crack in a hardening material.Journal of the Mechanics and Physics of Solids,vol.16,no.1,pp.13–31.

Kawai,H.(2006):ADVENTURE AutoGL:A handy graphics and GUI library for researchers and developers of numerical simulations.Computer Modeling in Engineering and Sciences,vol.11,no.3,pp.111–120.

Kelley,C.T.(2003):Solving Nonlinear Equations with Newton’s Method.Society for Industrial and Applied Mathematics,Philadelphia.

Minami,S.;Yoshimura,S.(2010):Performance evaluation of nonlinear algorithms with line-search for partitioned coupling techniques for fluid–structure interactions.International Journal for Numerical Methods in Fluids,vol.64,no.10–12,pp.1129–1147.

Miyamura,T.;Noguchi,H.;Shioya,R.;Yoshimura,S.;Yagawa,G.(2002):Elastic–plastic analysis of nuclear structures with millions of DOFs using the hierarchical domain decomposition method.Nuclear Engineering and Design,vol.212,no.1–3,pp.335–355.

Mo?s,N.;Dolbow,J.;Belytschko,T.(1999):A finite element method for crack grow th without remeshing.International Journal for Numerical Methods in Engineering,vol.46,pp.131–150.

Nakasumi,S.;Suzuki,K.;Fujii,D.;Ohtsubo,H.(2002):An elastic and elasto–plastic mixed analysis using overlaying mesh method[in Japanese].Transactions of the Japan Society of Mechanical Engineers,A,vol.68,no.668,pp.603–610.

Nikishkov,G.P.;Atluri,S.N.(1994):An analytical–numerical alternating method for elastic–plastic analysis of cracks.Computational Mechanics,vol.13,no.6,pp.427–442.

Nikishkov,G.P.;Park,J.H.;Atluri,S.N.(2001):SGBEM–FEM alternating method for analyzing 3D non-planar cracks and their grow th in structural components.Computer Modeling in Engineering and Sciences,vol.2,no.3,pp.401–422.

Nishikawa,H.;Serizawa,H.;Murakawa,H.(2007):Actualapplication of FEM to analysis of large scale mechanical problems in welding.Science and Technology of Welding and Joining,vol.12,no.2,pp.147–152.

Nishioka,T.;Atluri,S.N.(1983):Analytical solution for embedded elliptical cracks,and finite element alternating method for elliptical surface cracks,subjected toarbitrary loadings.Engineering Fracture Mechanics,vol.17,no.3,pp.247–268.

Ogino,M.;Shioya,R.;Kawai,H.;Yoshimura,S.(2005):Seismic response analysis of nuclear pressure vessel model with ADVENTURE system on the Earth Simulator.Journal of the Earth Simulator,vol.2,pp.41–54.

Okada,H.;Endoh,S.;K ikuchi,M.(2007):Application of s-version finite element method to two-dimensional fracture mechanics problems.Journal of Solid Mechanics of Material in Engineering,vol.1,no.5,pp.699–710.

Prabel,B.;Combescure,A.;Gravouil,A.;Marie,S.(2007):Level set X-FEM non-matching meshes:application to dynamic crack propagation in elastic–plastic media.International Journal for Numerical Methods in Engineering,vol.69,no.8,pp.1553–1569.

Pyo,C.R.;Okada,H.;Atluri,S.N.(1995):An elastic–plastic finite element alternating method for analyzing wide-spread fatigue damage in aircraft structures.Computational Mechanics,vol.16,no.1,pp.62–68.

Rice,J.R.;Rosengren,G.F.(1968):Plane strain deformation near a crack tip in a power-law hardening material.Journal of the Mechanics and Physics of Solids,vol.16,no.1,pp.1–12.

Samaniego,E.;Belytschko,T.(2005):Continuum–discontinuum modelling of shear bands.International Journal for Numerical Methods in Engineering,vol.62,no.13,pp.1857–1872.

Suzuki,K.;Ohtsubo,H.;Nakasumi,S.;Shinmura,D.(2002):Global local iterative analysis using overlaying mesh method[in Japanese].Journal of the Society of Naval Architects of Japan,vol.192,pp.691–696.

W hitcomb,J.D.(1991):Iterative global/local finite element analysis.Computers and Structures,vol.40,no.4,pp.1027–1031.

Yoshimura,S.;Shioya,R.;Noguchi,H.;Miyamura,T.(2002):Advanced general-purpose computational mechanics system for large-scale analysis and design.Journal of Computational and Applied Mathematics,vol.149,no.1,pp.279–296.

Yusa,Y.;Yoshimura,S.(2013):Mixed-mode fracture mechanics analysis of large-scale cracked structures using partitioned iterative coupling method.Computer Modeling in Engineering and Sciences,vol.91,pp.445–461.