Acessibilidade / Reportar erro

Global/Local Non-Intrusive Crack Analysis Using Element Free Galerkin and Linear Galerkin Finite Volume Methods

Abstract

Although various numerical methods are proposed to solve a wide range of problems so far, each has its own advantages and disadvantages. Therefore, it would be highly desirable to develop a computationally efficient combined approach, which makes use of their advantages and reduces their drawbacks. In this regard, a non-intrusive iterative global/local (IGL) algorithm is proposed to evaluate the structures with local discontinuity. Firstly, in each iteration, the entire continuum domain is solved using the Galerkin Finite Volume (GFV) method, which consumes light computational workload for predicting accurate results of linear structural response. Then, the Element Free Galerkin (EFG) meshless method is locally employed as an accurate but computationally expensive solver to cover the shortcomings of the GFV solver in the crack analysis. The specific form of Quasi-Newtonian and dynamic accelerators are adopted to increase the convergence rate in each computational increment. Finally, the proposed method's accuracy and speed are examined by computing stress intensity factors (SIFs) for several distinct 2D cracked models.

Keywords:
Iterative Global/Local Method; Element Free Galerkin Method; Galerkin Finite Volume Method; Non-intrusive Coupling; Crack Analysis; Localized Multi-Grid Patch

Graphical Abstract

1. Introduction

The crack analysis and its numerical simulation is among the most complex problems in computational mechanics, which needs accurate models and substantial amount of computational resources. Despite the rapidly increasing abilities and power of the supercomputers, the calculation costs could be still unaffordable in some complex cases. In this regard, it is essential to focus on analysis efficiency and accuracy, instead of computational power, to develop a flexible, simple to implement, and reliable solution. In this point of view, numerical techniques can be sort into two main categories: enhanced single-models and coupled multi-models.

Single-models such as generalized finite element (GFEM) (Kim et al. 2012Kim, D. J., Duarte, C. A. and Proenca, S. P. (2012) ‘A generalized finite element method with global-local enrichment functions for confined plasticity problems’, Computational Mechanics, 50(5), pp. 563-578. doi: 10.1007/s00466-012-0689-7.
https://doi.org/10.1007/s00466-012-0689-...
), extended finite element (XFEM (Fries and Belytschko, 2010Fries, T. P. and Belytschko, T. (2010) ‘The extended / generalized finite element method : An overview of the method and its applications’, Review Literature And Arts Of The Americas, (April), pp. 253-304. doi: 10.1002/nme.
https://doi.org/10.1002/nme...
)), and hybrid analytical and extended finite element method (HAX-FEM (Réthoré et al. 2012Réthoré, J., Roux, S. and Hild, F. (2012) ‘Hybrid analytical and extended finite element method (HAX-FEM): A new enrichment procedure for cracked solids’, International Journal for Numerical Methods in Engineering, 81(3), pp. 269-285. doi: 10.1002/nme.
https://doi.org/10.1002/nme...
)) were fundamentally introduced to reduce the high intrusiveness of adaptive methods. These methods solve the displacement discontinuity and the crack tip singularity, by implementing analytical enriched functions. Since the industrial sources of discontinuities are usually localized in extremely small regions, even if the discontinuous enrichments are well-defined, the original mesh size may be too coarse for precise crack simulation (Rangarajan and Lew, 2012Rangarajan, R. and Lew, A. J. (2012) ‘A local multigrid X-FEM strategy for 3-D crack propagation’, (August 2008), pp. 581-600. doi: 10.1002/nme.
https://doi.org/10.1002/nme...
). On the other hand, due to the single configuration of these models, the tiny local complexity may have global consequences and badly affects the overall efficiency. Therefore, although these enhanced techniques have conceptually proved their promising performances, they had practical challenges in industrial problems. In this sense, the coupling techniques can be classified as applied multi-model answers. These techniques are able to locally employ specific models, which could be fundamentally different from the basic model in terms of formulation, governing equations, or discretization techniques. Recently, some researchers began to combine methods in order to benefit from their advantages by trying to eliminate their disadvantages. Many details are provided by Zienkiewicz et al. (Zienkiewicz et al. 1977Zienkiewicz, O. C., Kelly, D. W. and Bettess, P. (1977) ‘The coupling of the finite element method and boundary solution procedures’, International Journal for Numerical Methods in Engineering, 11(2), pp. 355-375. doi: 10.1002/nme.1620110210.
https://doi.org/10.1002/nme.1620110210...
), von Estorff (Estorff and Antes, 1991Estorff, O. Von and Antes, H. (1991) ‘On FEM-BEM coupling for fluid-structure interaction analyses in the time domain’, International Journal for Numerical Methods in Engineering, 31(6), pp. 1151-1168. doi: 10.1002/nme.1620310609.
https://doi.org/10.1002/nme.1620310609...
), and Belytschko and Lu (Lu et al. 1991Lu, Y. Y., Belytschko, T. and Liu, W. K. (1991) ‘a variational coupled FE-BE method for elasticity and fracture mechanics’, 85, pp. 21-37.). The coupling algorithms are generally based on domain decomposition (DD) (El-Gebeily et al. 2002El-Gebeily, M., Elleithy, W. M. and Al-Gahtani, H. J. (2002) ‘Convergence of the domain decomposition finite element-boundary element coupling methods’, Computer Methods in Applied Mechanics and Engineering, 191(43), pp. 4851-4867. doi: 10.1016/S0045-7825(02)00405-X.
https://doi.org/10.1016/S0045-7825(02)00...
), on multi-grid (MG) algorithms (Henshaw, 2005Henshaw, W. D. (2005) ‘On Multigrid for Overlapping Grids’, SIAM Journal on Scientific Computing, 26(5), pp. 1547-1572. doi: 10.1137/040603735.
https://doi.org/10.1137/040603735...
), on extended/generalized finite elements (Gupta et al., 2015Gupta, V. et al. (2015) ‘Stable GFEM (SGFEM): Improved conditioning and accuracy of GFEM/XFEM for three-dimensional fracture mechanics’, Computer Methods in Applied Mechanics and Engineering, 289, pp. 355-386. doi: 10.1016/j.cma.2015.01.014.
https://doi.org/10.1016/j.cma.2015.01.01...
), on multi-scale projection (Holl et al. 2012Holl, M., Loehnert, S. and P. Wriggers (2012) ‘An adaptive multiscale method for crack propagation and crack coalescence’. doi: 10.1002/nme.
https://doi.org/10.1002/nme...
), on structural zooming (Hirai et al. 1984Hirai, I., Wang, B. P. and Pilkey, W. D. (1984) ‘An efficient zooming method for finite element analysis’, International Journal for Numerical Methods in Engineering, 20(9), pp. 1671-1683. doi: 10.1002/nme.1620200910.
https://doi.org/10.1002/nme.1620200910...
), on patch methods through Arlequin coupling (Ben Dhia, 1998Ben Dhia, H. (1998) ‘Multiscale mechanical problems: the Arlequin method’, Comptes Rendus De L Academie Des Sciences Serie Ii Fascicule B-Mecanique Physique Astronomie, 326(12), pp. 899-904. doi: doi:10.1016/S1251-8069(99)80046-5.
https://doi.org/doi:10.1016/S1251-8069(9...
) or on integral matching (Chahine et al. 2011Chahine, E., Laborde, P. and Renard, Y. (2011) ‘A non-conformal eXtended Finite Element approach: Integral matching Xfem’, Applied Numerical Mathematics, 61(3), pp. 322-343. doi: 10.1016/j.apnum.2010.10.009.
https://doi.org/10.1016/j.apnum.2010.10....
). Some coupling approaches mentioned so far were quite intrusive, which represents a restriction to use in large scale applications. Some others were expressed in an implicit direct coupling approach, which requires setting a monolithic system of equation. In this regard, the following drawbacks should be concerned: 1) Due to the direct coupling of the system matrices, the resultant equations are not sparsely populated anymore, which means the local optimization technique, cannot be generally used anymore. (2) In conventional transient direct coupling, the time step duration must be the same at each sub-domain, which causes numerical challenges, especially in the coupling of the different materials with completely different wave speeds. (3) In cases with discontinuous local sub-region, the rather large system of equations needs to be coupled and solved at each sub-iteration, which is extremely time-consuming.

In view of the first aspect, Whitcomb (Whitcomb, 1991Whitcomb, J. D. (1991) ‘Iterative global/local finite element analysis’, Computers and Structures, 40(4), pp. 1027-1031. doi: 10.1016/0045-7949(91)90334-I.
https://doi.org/10.1016/0045-7949(91)903...
) experienced an iterative coupling through Global/Local (IGL) configuration to take into account local effects over global background mesh. Whitcomb, for instance, modeled a global rectangular plate expected to remain elastic with a local circular patch assumed to behave nonlinearly by means of the FEM (Whitcomb and Woo, 1993a; 1993b). Elleithy et al. (Elleithy and Al-Gahtani, 2000Elleithy, W. M. and Al-Gahtani, H. J. (2000) ‘An overlapping domain decomposition approach for coupling the finite and boundary element methods’, Engineering Analysis with Boundary Elements, 24(5), pp. 391-398. doi: 10.1016/S0955-7997(00)00014-X.
https://doi.org/10.1016/S0955-7997(00)00...
) used the relaxation parameters on transfer region and represents various accelerators to converge the coupling iterations (Elleithy and Tanaka, 2002, 2003). Since they have struggled in elastostatic problems, it was very efficient concerning just the first and second viewpoints. Soares et al. (Soares, 2011Soares, D. (2011) ‘Coupled numerical methods to analyze interacting acoustic-dynamic models by multidomain decomposition techniques’, Mathematical Problems in Engineering, 2011. doi: 10.1155/2011/245170.
https://doi.org/10.1155/2011/245170...
; Godinho et al. 2013Godinho, L. et al. (2013) ‘Iterative coupling between the MFS and Kansa’s method for acoustic problems’, WIT Transactions on Modelling and Simulation, 54(April 2016), pp. 123-132. doi: 10.2495/BEM130121.
https://doi.org/10.2495/BEM130121...
) have successfully performed the iterative solver in the analysis of acoustic-dynamic interaction through the relaxation coefficients by multi-domain decomposition techniques (Soares et al. 2004). Their work is known as Neumann-Dirichlet algorithm in engineering literature.

Other techniques such as S-version (Fish, 1993Fish, J. (1993) ‘Adaptive S-method for Linear Elastostatics’, 104, pp. 363-396.), Spectral Overlay (Belytschko et al. 1990Belytschko, T., Fish, J. and Bayliss, A. (1990) ‘The Spectral Overlay on finite elements for problems with high gradients’, 81, pp. 71-89.), Mortar method (Bouclier et al. 2017Bouclier, R., Passieux, J. C. and Salaün, M. (2017) ‘Development of a new, more regular, mortar method for the coupling of NURBS subdomains within a NURBS patch: Application to a non-intrusive local enrichment of NURBS patches’, Computer Methods in Applied Mechanics and Engineering, 316, pp. 123-150. doi: 10.1016/j.cma.2016.05.037.
https://doi.org/10.1016/j.cma.2016.05.03...
), and Arlequin method (energy averaging scheme) (Ben Dhia, 1998Ben Dhia, H. (1998) ‘Multiscale mechanical problems: the Arlequin method’, Comptes Rendus De L Academie Des Sciences Serie Ii Fascicule B-Mecanique Physique Astronomie, 326(12), pp. 899-904. doi: doi:10.1016/S1251-8069(99)80046-5.
https://doi.org/doi:10.1016/S1251-8069(9...
; Bauman et al., 2008Bauman, P. T. et al. (2008) ‘On the application of the Arlequin method to the coupling of particle and continuum models’, Computational Mechanics, 42(4), pp. 511-530. doi: 10.1007/s00466-008-0291-1.
https://doi.org/10.1007/s00466-008-0291-...
) are also multi-model approaches. These techniques are concerning the interfacial boundary values in the global/local framework. Winh phu Nguyen (Nguyen-Xuan et al., 2012Nguyen-Xuan, H. et al. (2012) ‘A cell-based smoothed finite element method for three dimensional solid structures’, KSCE Journal of Civil Engineering.) has recently carried out some researches on the molecular dynamics using Arlequin framework. Park and Felippa (2000Park, K. C. and Felippa, C. A. (2000) ‘A variational principle for the formulation of partitioned structural systems’, International Journal for Numerical Methods in Engineering, 47(1-3), pp. 395-418. doi: 10.1002/(SICI)1097-0207(20000110/30)47:1/3<395::AID-NME777>3.0.CO;2-9.
https://doi.org/10.1002/(SICI)1097-0207(...
) used Partitioned Analysis which was a kind of multi-level and hierarchical DD method, using Lagrange multipliers to satisfy the compatibility condition for interface boundary deformations. Hagen and Von Estorff (Von Estorff and Hagen, 2006) were conducted extensive parametric studies to achieve an ideal relaxation parameter based on the interface displacements and unbalanced forces.

Despite all these efforts and progress, the computational efficiency and accuracy of coupling techniques are still an open issue that needs more investigations to reach an efficient model with lower intrusiveness (this is the point we will focus on in this paper). In this regard, a new class of non-intrusive coupling techniques has recently been established by Gendre et al. (Gendre et al., 2009Gendre, L. et al. (2009) ‘Non-intrusive and exact global / local techniques for structural problems with local plasticity’, 44, pp. 233-245. doi: 10.1007/s00466-009-0372-9.
https://doi.org/10.1007/s00466-009-0372-...
). Later, it has been developed by Duval and Passieux (Duval et al., 2014Duval, M. et al. (2014) ‘Local/Global Non-Intrusive Parallel Coupling for Large Scale Mechanical Analysis’, 343(July), pp. 1210-1211. doi: 10.1126/science.1249625.
https://doi.org/10.1126/science.1249625...
; 2016b) to solve crack problems. This method is generally dedicated to large-scale finite element models, encountered the local phenomenon (such as local mesh refinement, local nonlinearity, or local discontinuity), and supported by the local module due to the lack of competence in the global model. The purpose is to generate a non-intrusive global/local coupling, which allows us to implement the necessary modifications without any changes carrying out to the corresponding mesh or solver. This technique seems restrictive, nevertheless enables one to couple numerous software such as Abaqus and Ls-Dyna or research black box codes together without any fundamental modifications.

In this paper, the non-intrusive localized IGL coupling scheme is presented to improve the overall efficiency in the mixed-mode crack analysis. The Galerkin finite volume (GFV) (Alkhamis et al., 2008Alkhamis, M. T. et al. (2008) ‘Utilizing NASIR galerkin finite volume analyzer for 2D plane strain problems under static and vibrating concentrated loads’, Jordan Journal of Civil Engineering, 2(4), pp. 335-343.; Sabbagh-Yazdi et al. 2012Sabbagh-Yazdi, S. R., Ali-Mohammadi, S. and Pipelzadeh, M. K. (2012) ‘Unstructured finite volume method for matrix free explicit solution of stress-strain fields in two dimensional problems with curved boundaries in equilibrium condition’, Applied Mathematical Modelling, 36(5), pp. 2224-2236. doi: 10.1016/j.apm.2011.08.001.
https://doi.org/10.1016/j.apm.2011.08.00...
; 2019; Sabbagh-Yazdi and Najar-Nobari, 2020) is globally chosen as the fast and robust technique to linearly discretize the entire domain, while the small critical zone is served with the EFG (Belytschko et al. 1994Belytschko, T., Lu, Y. Y. and Gu, L. (1994) ‘Element-Free Galerkin methods’, International Journal, 37(April 1993), pp. 229-256. doi: 10.1002/nme.1620370205.
https://doi.org/10.1002/nme.1620370205...
; 1995). Both methodologies have relative benefits and drawbacks. The GFV, for instance, is appropriate for linear analysis of the homogeneous and isotropic materials. Nevertheless, the EFG is noticeably more advantageous regarding the nonlinearity, stress concentration, or discontinuity. Here, we have explained how to embed the EFG module within a global linear GFV solver, both considered as black-box models. Afterward, an enhanced version of the non-intrusive approach is explain that the relaxation parameters would match with the material constitutive matrices of each global element. We will continue to examine the efficiency and performance of the proposed coupling scheme through two numerical examples, namely an infinite tensile plate with one straight center crack, and a finite plate with one straight edge crack under uniaxial load. Considering linear fracture mechanics, it is shown that the proposed methodology provides acceptable outcomes and is both fast and robust.

2. The General Principle of Explicit Overlapping Coupling

The proposed explicit overlapping scheme begins with the GFV model, which is called the global model and is supposed to completely behave linearly. In order to consider the discontinuity, the EFG has been locally considered as a more reliable and specific method. Since the GFV is unable to explain the discontinuity, this method is just valid in the complementary area. The local model can be a part of the global model (Fig. 1a) or be a separate set of nodes (Fig. 1b), which may result in more accurate solutions. Since two separate models may fundamentally non-matching, robust mapping operators must extensively be performed in the interface called ΩG|L . These mapping operators, called prolongation and restriction operators, must be well chosen interpolating functions in order to prevent the validation loss.

Fig. 1
Schematic overlapping global/local configuration

The solution algorithm contains the global and local discrete formulations, the non-intrusive coupling techniques, the mapping operators, and the enhanced accelerators regarding the global/local data sets. In the following sections, we tried to explain each part briefly.

2.1. Galerkin finite volume (GFV) solver

As mentioned before, the GFV is considered as the global model which was first introduced by Sabbagh-Yazdi et al. (Alkhamis et al., 2008Alkhamis, M. T. et al. (2008) ‘Utilizing NASIR galerkin finite volume analyzer for 2D plane strain problems under static and vibrating concentrated loads’, Jordan Journal of Civil Engineering, 2(4), pp. 335-343.; Sabbagh-Yazdi et al. 2012; 2019; Sabbagh-Yazdi and Najar-Nobari, 2020) as an explicit matrix free module, which is an iterative cell-vertex FV solver. GFV is well-examined, and its promising results are proven in various mechanical problems such as dynamic analysis (Sabbagh-Yazdi et al. 2019), crack growth (Sabbagh-Yazdi and Bayatlou, 2012), dynamical failure (Amraei and Fallah 2016Amraei, A. and Fallah, N. (2016) ‘A development in the finite volume method for the crack growth analysis without global remeshing’, International Journal of Engineering, Transactions A: Basics, 29(7), pp. 890-900. doi: 10.5829/idosi.ije.2016.29.07a.03.
https://doi.org/10.5829/idosi.ije.2016.2...
), and thermal elasticity (Sabbagh-Yazdi and Amiri-SaadatAbadi, 2013). The governing equation of the GFV is returned to the Cauchy’s equilibrium equations as follows:

ρ 2 u i t 2 = σ x i + b i (1)

Through the Variational method, by multiplying the residual of the above equation by the weight function (щ), and performing the partial integration, the Galerkin weak form obtains as (Sabbagh-Yazdi and Najar-Nobari, 2020Sabbagh-Yazdi, S. R. and Najar-Nobari, H. (2020) ‘Coupling Nonlinear Element Free Galerkin and Linear Galerkin Finite Volume Solver for 2D Modeling of Local Plasticity in Structural Material’, International Journal of Engineering, 33(3), pp. 387-400. doi: 10.5829/ije.2020.33.03c.03.
https://doi.org/10.5829/ije.2020.33.03c....
):

Ω ω . ρ 2 u i t 2 d Ω = [ ω . F i ] Γ Ω ( F i . ω ) d Ω + Ω ω . b i d Ω (2)

Since the weight function (ω) and the linear interpolation function (ϕ) are the same, the [ω.Fi]Γ in equation (2), would be zero-valued on the boundary edges of each control volume (Sabbagh-Yazdi et al. 2018Sabbagh-Yazdi, S. R., Farhoud, A. and Asil Gharebaghi, S. (2018) ‘Simulation of 2D linear crack growth under constant load using GFVM and two-point displacement extrapolation method’, Applied Mathematical Modelling. Elsevier Inc., 61, pp. 650-667. doi: 10.1016/j.apm.2018.05.022.
https://doi.org/10.1016/j.apm.2018.05.02...
). In addition, the second term in this equation is discretized as follows:

Ω ( F i . ω ) d Ω 1 2 k = 1 N ( F i . Δ l ) k = k = 1 N ( σ ˜ i 1 Δ x 2 σ ˜ i 2 Δ x 1 ) k (3)

Finally, applying the finite difference method according to the GFV algorithm (Sabbagh-Yazdi and Amiri-SaadatAbadi, 2013Sabbagh-Yazdi, S. R. and Amiri-SaadatAbadi, T. (2013) ‘GFV solution on UTE mesh for transient modeling of concrete aging effects on thermal plane strains during construction of gravity dam’, Applied Mathematical Modelling, 37(1-2), pp. 82-101. doi: 10.1016/j.apm.2012.01.045.
https://doi.org/10.1016/j.apm.2012.01.04...
), and by putting equation (3) into (2), we have the following explicit form:

( u i ) n k + 1 = 2 ( u i ) n k ( u i ) n k 1 + 3 ( Δ t n ) 2 2 ρ Щ n [ k = 1 N ( σ ˜ i 1 Δ x 2 σ ˜ i 2 Δ x 1 + F i ) k + 2 b i Щ n 3 ] t (4)

σ i j = C u i , j 1 A n m = 1 3 C { u x Δ y u y Δ x ( u x Δ y u y Δ x ) } m (5)

In which, (ui)nk+1 is the displacement of the n th node in the i th direction at (k+1)th iteration, and An (as shown in Fig. 2) is the cumulative area of the surrounding cells, associated to node n.

Fig. 2
The area of the linear triangular elements at each control volume

It is necessary to notice that the Δtn is a virtual parameter without any physical interpretation, which is essential for equilibrium in iterative scheme. To enhance the convergence speed and stabilize the GFV explicit solution, one should evaluate this virtual time increment by the following value (Sabbagh-Yazdi et al. 2012Sabbagh-Yazdi, S. R., Ali-Mohammadi, S. and Pipelzadeh, M. K. (2012) ‘Unstructured finite volume method for matrix free explicit solution of stress-strain fields in two dimensional problems with curved boundaries in equilibrium condition’, Applied Mathematical Modelling, 36(5), pp. 2224-2236. doi: 10.1016/j.apm.2011.08.001.
https://doi.org/10.1016/j.apm.2011.08.00...
).

Δ t n r n C w a v e (6)

where rn is the proportion of area to perimeter in each control volume (Sabbagh-Yazdi et al. 2018Sabbagh-Yazdi, S. R., Farhoud, A. and Asil Gharebaghi, S. (2018) ‘Simulation of 2D linear crack growth under constant load using GFVM and two-point displacement extrapolation method’, Applied Mathematical Modelling. Elsevier Inc., 61, pp. 650-667. doi: 10.1016/j.apm.2018.05.022.
https://doi.org/10.1016/j.apm.2018.05.02...
).

2.2. Element free Galerkin (EFG) as a local patch

The equation of equilibrium for two-dimensional, homogeneous, isotropic, and linear-elastic solids is written as follows:

L T σ + b = 0 i n Щ σ n = t ¯ o n Г t u = u ¯ o n Г u (7)

where L is the differential operator, t¯ and u¯ are the prescribed traction and displacement vectors on the Neumann and Dirichlet boundary edges, respectively. The weak form of the equation (7) can be assembled as:

Ω ( L δ u ) T ( D L u ) d Ω Ω δ u T b d Ω Γ t δ u T t ¯ d Ω = 0 (8)

where D is the matrix of elastic constants. In EFG method, the MLS (Lu et al. 1994Lu, Y. Y., Belytschko, T. and Gu, L. (1994) ‘A new implementation of the element free Galerkin method’, Computer Methods in Applied Mechanics and Engineering, 113(3-4), pp. 397-414. doi: 10.1016/0045-7825(94)90056-6.
https://doi.org/10.1016/0045-7825(94)900...
) shape functions are used to approximate the field variables at any point of interest. To simulate the crack problems in EFG, the XFEM enrichment functions are implemented to efficiently analyze discontinuities and singularities. In this regard, a few higher order Gaussian points are added to the near crack points which contribute to the overall approximation through the use of enrichment functions. Therefore, the approximated displacement (uh(x)) of a point of interest x could be written as (Ghorashi et al. 2010Ghorashi, S. S., Sabbagh-Yazdi, S. R. and Mohammadi, S. (2010) ‘Element free Galerkin method for crack analysis of orthotropic plates’, 1(1), pp. 1-13.):

u h ( x ) = i = 1 n φ i ( x ) u i + j = 1 n c φ j ( x ) [ H ( x ) a i + k = 1 4 Q k ( x ) b i k ] (9)

where ai and bik are the additional nodal dof near the crack edge and tip, associated to the enrichment functions; and ui is the vector of regular nodal dof. After manipulating equation (8) by equation (9), the final discretized system equations for the developed discontinuous enriched EFG has been derived as:

K U = F (10)

where K and F are assembled as follows:

K i j n = [ K i j u u K i j u b K i j b u K i j b b ] (11)

F i n = [ F i u F i b α F i b α F i b α F i b 4 ] (12)

which

K i j u b = Ω ( B i u ) T D B j b d Ω (13)

F i u = Ω φ i T b d Ω + Ω φ i T t ¯ d Г (14)

F i b α = Ω φ i T Q a b d Ω + Ω φ i T Q a t ¯ d Г (15)

Since the Kronecker delta properties are not satisfied by the MLS shape functions, in this study, the Lagrange multiplier method (Belytschko et al. 1994Belytschko, T., Lu, Y. Y. and Gu, L. (1994) ‘Element-Free Galerkin methods’, International Journal, 37(April 1993), pp. 229-256. doi: 10.1002/nme.1620370205.
https://doi.org/10.1002/nme.1620370205...
) is applied to effectively enforce the essential boundary conditions. Furthermore, the diffraction method (Organ et al., 1996Organ, D. et al. (1996) ‘Continuous meshless approximations for nonconvex bodies by diffraction and transparency’, Computational Mechanics, 18(3), pp. 225-235. doi: 10.1007/BF00369940.
https://doi.org/10.1007/BF00369940...
) has been used to detect the support domain for points on opposite sides of the crack as shown in Fig. 3.

Fig. 3
Support domain detection

2.2.1. Stress Intensity Factors (SIFs)

The stress intensity factor (SIF) is one of the most significant parameters for introducing the fracture properties of a crack tip field. The standard path independent J -integral for a crack is defined as:

J = Г [ W s δ i j σ i j u i x i ] n j d Г (16)

where Г is an arbitrary domain surrounding the crack-tip, and n j is the j th component of the outward unit normal to Г. In this paper, the technique established by Kim and Paulino (Kim and Paulino, 2002Kim, J.-H. and Paulino, G. H. (2002) ‘Finite element evaluation of mixed mode stress intensity factors in functionally graded materials’, International Journal for Numerical Methods in Engineering, 53(8), pp. 1903-1935. doi: 10.1002/nme.364.
https://doi.org/10.1002/nme.364...
) to evaluate the mixed-mode SIF is applied and briefly reviewed. Many details are given, e.g., by Ghorashi and Sabbagh-Yazdi (Ghorashi et al. 2010Ghorashi, S. S., Sabbagh-Yazdi, S. R. and Mohammadi, S. (2010) ‘Element free Galerkin method for crack analysis of orthotropic plates’, 1(1), pp. 1-13.; 2011). In this technique, two equilibrium states (actual state and auxiliary state) have been employed to obtain more accurate results. The actual and auxiliary states are obtained by near tip asymptotic fields as given by numerical and analytical solutions, respectively. By combining an actual and auxiliary solutions, the equation (17) would be obtained as follows (Ghorashi et al. 2011):

J ( A c t + A u x ) = J ( A c t ) + J ( A u x ) + M (17)

where J(Act+Aux) is the mixed-mode J-integral value, J(act) and J(aux) are the J-integral values for actual and auxiliary states, as follow:

J ( a c t ) = K I a c t 2 / E + K I I a c t 2 / E J ( a u x ) = K I a u x 2 / E + K I I a u x 2 / E (18)

In equation (17), M is equal to:

M = Γ [ W s ( A c t + A u x ) δ i j σ i j u i a u x x 1 σ i j a u x u i a c t x 1 ] n j d s (19)

Using the auxiliary fields along the analytical field, the J-integral could be written as follows:

J ( A c t + A u x ) = Γ [ W s ( A c t + A u x ) δ i j ( σ i j A c t + σ i j A u x ) ( u i A c t + u i A u x ) ] n j d s (20)

In which Ws(Act+Aux) is the strain energy density for linear elastic conditions as follows:

W s ( A c t + A u x ) = 1 2 ( σ i j A c t ε i j A u x + σ i j A u x ε i j A c t ) (21)

For mixed-mode problems in linear elastic fracture mechanics (LEFM), the J-integral is related to SIFs derived from the Irwin relation as follows:

J = 1 E * [ ( K I a c t + K I a u x ) 2 + ( K I I a c t + K I I a u x ) 2 ] = 1 E * ( K I a c t 2 + K I I a c t 2 ) + 1 E * ( K I a u x 2 + K I I a u x 2 ) + 2 E * ( K I a c t K I a u x + K I I a c t K I I a u x ) = J ( a c t ) + J ( a u x ) + M (22)

2.3. Iterative Global/Local Coupling

Following the non-intrusive framework, the global/local algorithm just depends on the global and local responses without any direct matrix manipulations. The basic force equilibrium and the displacement compatibility to pair ΩL and ΩG, are given as follows:

u = { u L o c a l i n c r a c k e d a r e a u G l o b a l e v e r y w h e r e e l s e u L = u G o n Г u σ L n l + σ G n g = 0 o n Г t (23)

where nl and ng are normal boundary vectors on local and global boundary nodes. To solve the above equations both accurately and non-intrusively, implementing an iterative scheme could be the appropriate choice, as suggested by Whitcomb and Passieux (Whitcomb and Woo, 1993; Passieux et al. 2013). In this regard, the following four iterative steps must be fulfilled on the interfaces (Passieux et al., 2013):

a) Global initial solution: Solve the problem assuming the entire global domain, except the near crack area, is in the elastic condition. The modulus of elasticity is considered to be zero in cracked global elements for two first steps. Intuitively, it represents an appropriate initial estimate for the following steps. According to this step, the global solution is appeared in the following form:

K G U G k = F G k 1 + Re s k 1 Re s k 1 = [ R K L U L k 1 2 ] [ K G | L U G k 1 ] (24)

where FGk1 and UGk1 are the global force and displacement vectors in (k1)th step. The Res contains the F G and two additional terms to inactivate the KG|LUGk1 by the restriction of the KLULk1/2. Multiplying KG1 by both sides of Equation (24), we get:

U G k = U G k 1 + K G 1 Re s k 1 (25)

In equation (24), K L and K G are the local and global stiffness matrices. In commercial software and explicit solvers such as Matrix free GFV, practically the stiffness matrix is unachievable. Here, based on the GFV explicit solution in equation (4), the displacements are the function of the last two sub-iterations as below:

u t + Δ t = m = 1 N U ˜ ( u t , u t Δ t , σ , f , b ) (26)

Then, according to the Hooke’s law, equation (25) can be revised as follows:

U G k = U G k 1 + m = 1 N U ˜ ( { [ P T U L k 1 2 ] U G k 1 } , σ G k 1 , F i , b i ) (27)

If the local stresses are mapped onto the global domain as boundary tractions:

F G k = F G k 1 + m = 1 N F ˜ ( U G k 1 , { [ P T K L u L k 1 2 ] F G k 1 } , b i ) (28)

b) Solution on local domain: The local cracked problem is solved with prescribed values, which are extracted from the current global solution. The displacement compatibility can better be guaranteed by the essential boundary values (prescribed displacements); however, it is not necessarily the most reliable one. Assuming the displacement prolongation, the local solution ULk+1/2is improved as follows (Sabbagh-Yazdi and Najar-Nobari, 2020Sabbagh-Yazdi, S. R. and Najar-Nobari, H. (2020) ‘Coupling Nonlinear Element Free Galerkin and Linear Galerkin Finite Volume Solver for 2D Modeling of Local Plasticity in Structural Material’, International Journal of Engineering, 33(3), pp. 387-400. doi: 10.5829/ije.2020.33.03c.03.
https://doi.org/10.5829/ije.2020.33.03c....
):

U L k + 1 2 = U L k 1 2 + P Δ U G k Δ U G k = U G k R U L k 1 2 (29)

In this article, (k12) always refers to local step between (k1)th and (k)th global steps, and (k+12) refers to local step between (k)th and (k+1)th global steps. Multiplying KL by both sides of Equation (29), we get:

F L k + 1 2 = F L k 1 2 + P Δ T G k Δ T G k = [ P σ G k σ L k 1 2 ] | Г t . n l i (30)

c) Residual computation: Convergence fundamentally depends on the interface equilibrium residual of the overlapping domains. In order to reach convergence, the interfacial residual’s norm (η) should be minimized as follows:

η = U G k + 1 U G k U G k + 1 i f η < ε : { s t o p } ; e l s e : { k k + 1 } (31)

d) Global modification: If the residual norm were not small enough, it would be supplemented with the global values as a corrective interfacial term. Since this algorithm is similar to the modified Newton iteration, the efficiency of this step can be easily improved by using some accelerators such as SR1 (Allix et al., 2011Allix, O. et al. (2011) ‘Non-intrusive coupling: An attempt to merge industrial and research software capabilities’, Recent Developments and Innovative Applications in Computational Mechanics, pp. 125-133. doi: 10.1007/978-3-642-17484-1_15.
https://doi.org/10.1007/978-3-642-17484-...
) or Aitken (Küttler and Wall, 2008Küttler, U. and Wall, W. A. (2008) ‘Fixed-point fluid-structure interaction solvers with dynamic relaxation’, Computational Mechanics, 43(1), pp. 61-72. doi: 10.1007/s00466-008-0255-5.
https://doi.org/10.1007/s00466-008-0255-...
) as below:

{ U G | L k + 1 } = U G k + 1 | Г u = ( 1 α ) P T U L k + 1 2 | Г u + α U G k | Г u (32)

{ F G | L k + 1 } = F G k + 1 | Г t = ( 1 β ) P T σ L k + 1 2 | Г t . n j i G + β F G k | Г t (33)

where α and β are relaxation parameters to accelerate the convergence. {UG|Lk+1} and {FG|Lk+1} are the modified local interfacial values on global domain.

The above iterative scheme has some noticeable preferences. Firstly, it is reliable because the accuracy can directly be enhanced by introducing the desired extent for the residual norms. Secondly, the algorithm is considerably non-intrusive due to the overlapping sub-domains act like the black-boxes. Thirdly, it is light and easy to implement, because it only requires to send and receive the supplementary interfacial nodal displacement or force residuals, which are conventional numerical outputs, between global and local solvers.

Regarding the mentioned iterative scheme, the convergence algorithm based on the Neumann-Dirichlet procedure (El-Gebeily et al. 2002El-Gebeily, M., Elleithy, W. M. and Al-Gahtani, H. J. (2002) ‘Convergence of the domain decomposition finite element-boundary element coupling methods’, Computer Methods in Applied Mechanics and Engineering, 191(43), pp. 4851-4867. doi: 10.1016/S0045-7825(02)00405-X.
https://doi.org/10.1016/S0045-7825(02)00...
) can be summarized as Table 1. This algorithm solves the GFV and the EFG sub-system separately, which means that separate solving procedures can be employed to solve each system of equations. Thus, the symmetry and sparsity of the EFG matrices can easily be taken into account, which results in a more efficient methodology.

By solving the GFV and the EFG apart, one also has an optimized system of equations, which is significant regarding the analytical accuracy and performance. It should be also pointed out that when a proposed algorithm is considered, which includes a renewal of the EFG matrices at different iterative steps, the GFV system of equations is not affected by this renewal. In this way, a substantial amount of calculations is avoided.

In this procedure, it is essential to initially consider the zero value for modulus of elasticity of cracked global elements in continuous global models, to enhance the convergence and avoid interfacial distortion of the first two steps.

Table 1
The Iterative global/local algorithm based on the Neumann-Dirichlet procedure

2.3.1. Convergence Properties - Relaxation

The convergence of the iterative processes is the most significant part, in terms of speed and accuracy. In Equations (32) and (33), applying the well-chosen relaxation parameters (α or β), remarkably enhances the stability of the iterative scheme and ensures the convergence of it, even if ΩG is stiffer than ΩL. The enhanced global value consists in an appropriate portion of the predicted value (PTuLk+1/2) and the previous computed global value (uGk), by implementing the μ(k). One of the most efficient techniques is the quasi-Newtonian accelerator formulated based on the Symmetric Rank One (SR1) formula, extensively employed in iterative schemes (Allix et al., 2011Allix, O. et al. (2011) ‘Non-intrusive coupling: An attempt to merge industrial and research software capabilities’, Recent Developments and Innovative Applications in Computational Mechanics, pp. 125-133. doi: 10.1007/978-3-642-17484-1_15.
https://doi.org/10.1007/978-3-642-17484-...
). Although it enhances the convergence speed, the updating process of the global stiffness is complicated. Another category of acceleration techniques that modifies interfacial terms without updating the global stiffness matrix is Aitken's Δ2 accelerator (Küttler and Wall 2008Küttler, U. and Wall, W. A. (2008) ‘Fixed-point fluid-structure interaction solvers with dynamic relaxation’, Computational Mechanics, 43(1), pp. 61-72. doi: 10.1007/s00466-008-0255-5.
https://doi.org/10.1007/s00466-008-0255-...
), which is one of the most effective ones.

2.3.2. Dynamic relaxation: Aitken’s Δ2 acceleration

This method is most beneficial for speeding up the iterations that are converging linearly according to the results of the two previous steps. This parameter is utilized, as follows:

u G k + 1 = μ ( k ) P T u L k + 1 2 + ( 1 μ ( k ) ) u G k (34)

where μ(k) is the relaxation coefficient in the kth step, obtained from the results of the two previous steps, uGk+1 is the displacement component used in the next global step and PTuLk+1/2 is the processed local displacement at the end of k th step. In order to estimate the optimum relaxation factor, the following functional (Godinho et al., 2013Godinho, L. et al. (2013) ‘Iterative coupling between the MFS and Kansa’s method for acoustic problems’, WIT Transactions on Modelling and Simulation, 54(April 2016), pp. 123-132. doi: 10.2495/BEM130121.
https://doi.org/10.2495/BEM130121...
) is considered:

f ( μ ( k ) ) = u G k + 1 u G k 2 (35)

Assuming PTuLk+1/2PTuLk1/2=Δ(PTuLk+1/2) and uGkuGk1=ΔuGk then by manipulating equation (34) by equation (35), the functional is extracted as follows:

f ( μ ( k ) ) = μ ( k ) Δ ( P T u L k + 1 / 2 ) + ( 1 μ ( k ) ) Δ u G k 2 (36)

The relaxation parameter (μ(k)) is optimized by deriving from μ(k) as follow:

f ( μ ( k ) ) μ ( k ) = 2 μ ( k ) Δ ( P T u L k + 1 2 ) + ( 1 μ ( k ) ) Δ u G k Δ ( P T u L k + 1 2 ) Δ u G k = 0 (37)

Then μ(k) is updated through the recurrent formula as follows:

μ ( k + 1 ) = μ ( k ) ( Δ u G k . ( Δ u G k Δ ( P T u L k + 1 2 ) ) ) Δ u G k Δ ( P T u L k + 1 2 ) 2 (38)

Since the Aitken’s Δ2 method employs just interfacial displacement or stress components from two preceding iterations, the resultant parameter is computationally efficient and fast.

2.3.3. Quasi-Newtonian Acceleration

Another feasible approach to accelerate the convergence is the Quasi-Newtonian method to renew the stiffness matrix. The SR1 is an adaptable and practical updating scheme which generates a series of updated stiffness matrices as follows (Duval et al., 2016aDuval, M. et al. (2016a) ‘Non-intrusive Coupling: Recent Advances and Scalable Nonlinear Domain Decomposition’, Archives of Computational Methods in Engineering, 23(1), pp. 17-38. doi: 10.1007/s11831-014-9132-x.
https://doi.org/10.1007/s11831-014-9132-...
):

K G k + 1 = K G k + ( Δ f G k K G k Δ u G k ) ( Δ f G k K G k Δ u G k ) T ( Δ u G k ) T ( Δ f G k K G k Δ u G k ) K G k + ( f G k + 1 ) ( f G k + 1 ) T ( Δ u G k ) T ( f G k + 1 ) (39)

where ΔuGk=uGk+1uGk and ΔfGk=fGk+1fGk. Since the stiffness matrix is ​​fundamentally unavailable in the GFV, due to the matrix-free structure, one can use the above equation to update the material constitutive matrix instead of the stiffness matrix as follow:

D G k + 1 = D G k + ( Δ σ G k D G k Δ ε G k ) ( Δ σ G k D G k Δ ε G k ) T ( Δ ε G k ) T ( Δ σ G k D G k Δ ε G k ) D G k + ( σ G k + 1 ) ( σ G k + 1 ) T ( Δ ε G k ) T ( σ G k + 1 ) (40)

The DGk+1 obtained from Equation (40) could be implemented to compute global stress and displacement based on Hooke’s law on (k+1)th step.

2.3.4. Nodal Conformity and Data Transition Between Non-Matching Grids

In cases, which local and global nodes are not matched on the interface, the data-mapping operators between these two sets of non-matching nodes should be implemented. These mapping tools are called prolongation and restriction operators, which are applied to many multi-model and multi-domain methods. In this research, the MLS approximation is extensively used, both to construct EFG shape functions and data transmission between local and global interfaces as prolongation and restriction operators. Fig. 4 shows how MLS can effectively map displacement data between two sets of nodes.

Fig. 4
Nodal conformity by MLS approximation

In this regard, the transition region is not necessary to be a line and could be a wide strip of boundary nodes.

3. Numerical Examples

In this section, two crack problems are examined to express the precision and performance of the proposed method. Each example is investigated based on the interfacial parameters through six different models, as shown in Table 2. These models are classified according to their global/local configurations and the accelerator types. In order to control the accuracy of the proposed methods, various criteria such as SIFs and the L2 error of nodal displacements and stresses, are used and compared with exact solutions. In this paper, the unstructured triangle meshes are generated with the program “Gmsh”, whereas the local triangle meshes are created by an open-source Matlab code “DistMesh” from the University of California, Berkeley. All post-processing and numerical operations, e.g. displacement and stress contours are performed in Matlab.

Table 2
Non-intrusive classification based on the interface parameters

3.1. Center crack under far field uniform tensile loading

In this example, the infinite plate with a central crack is considered. The crack length is 2a, and the uniform far-field stress of у0 is applied vertically as illustrated in Fig. 5. In this case, just the finite closed square of 10×10 mm (ABCD) is simulated. The center of this area is located at the crack tip. The material and geometric specifications of this example are as follows: υ = 0.3 E = 200 GPa a = 100 mm e = 5mm

Fig. 5
Geometric specifications of infinite plate with central crack

In this case, the analytical solution is available. Considering the polar coordinates at the crack tip, the analytical displacement fields (Ghorashi et al. 2011Ghorashi, S. S., Mohammadi, S. and Sabbagh-Yazdi, S. R. (2011) ‘Orthotropic enriched element free Galerkin method for fracture analysis of composites’, Engineering Fracture Mechanics. Elsevier Ltd, 78(9), pp. 1906-1927. doi: 10.1016/j.engfracmech.2011.03.011.
https://doi.org/10.1016/j.engfracmech.20...
) are computed as follows:

u x ( r , θ ) = 2 ( 1 + ν ) 2 π K I E r cos θ 2 ( 2 2 ν cos 2 θ 2 ) u y ( r , θ ) = 2 ( 1 + ν ) 2 π K I E r sin θ 2 ( 2 2 ν cos 2 θ 2 ) K I = у 0 π a (41)

And the analytical stress fields are:

у x x ( r , θ ) = K I 2 π r cos θ 2 ( 1 sin θ 2 sin 3 θ 2 ) у y y ( r , θ ) = K I 2 π r cos θ 2 ( 1 + sin θ 2 sin 3 θ 2 ) у x y ( r , θ ) = K I 2 π r sin θ 2 cos θ 2 cos 3 θ 2 (42)

In this case, due to the availability of the analytical solution, the finite square (ABCD) is modeled with prescribed boundary conditions according to the analytical values. Therefore, the displacement fields are enforced at the right, top, and bottom edges, and the prescribed force is applied at the left edge as Dirichlet and Neumann boundary conditions, respectively.

The mode I SIF errors for various J-integral radiuses are computed and given in Fig. 6. It is obtained by a local circular shape with a radius of 2 mm, support domain with radius of rs, and various J-integral radiuses of rj. A sensitivity analysis proved that circular support domains with radius of 1.7Д lead to the best results in this problem, where Д is the local average nodal spacing.

Fig. 6
Mode I SIFs errors for various rj/rL values after 3 iterations (Local domain is a circle with a radius of 2mm)

From this figure, it is quite clear that when rj/rL0.8, the error is dramatically increased, due to the interference of the local and J-integral domains. Also at lower range of rj/rL, the large amount of error is observed. It is due to the vicinity of the J-integral radius (rj) and the enrichment radius (renr). In fact, good results are obtained with renr<rj<0.8rL.

Generally, it is found that the proposed method has agreeable results in comparison with the exact solution.

3.2. Finite Plate with an Edge Crack Subjected to Uniform Tension

In the following example, a square plate with a single edge crack and under a uniform tensile load is considered. The geometry of the plate and the crack are shown in in Fig. 7. The axial tensile load is uniformly applied on top and bottom edges. The bottom right corner of the plate is fixed along the x-axis, and the top right corner is hinged.

Fig. 7
Geometry and loading of the plate with an edge crack

The force and material properties are as follows: υ = 0.3 E = 107 MPa у0=103KPa

The analytical stress field for this example is same as the previous example, but the reference SIFs should be corrected by the following factor depends on a/b given by Tada, Paris, and Irwin (Tada et al. 1973Tada, H., Paris, P. C. and Irwin, G. (1973) ‘The Stress Analysis of Cracks Handbook’, Del Research Corporation, Hellertown, PA.) as:

K I r e f = у 0 π a [ 1.12 0.231 ( a b ) + 10.55 ( a b ) 2 21.72 ( a b ) 3 + 30.39 ( a b ) 4 ] K I I r e f = 0 M P a m (43)

where a and b are crack length and plate width, respectively. The normalized SIF for mode one for the present problem is equal to K¯I=KIref/у0рa=2.8264.

In this example, the initial global model has an unstructured mesh made of 1058 linear triangular elements and 580 nodes, and four different local models are considered for M1 to M6 as depicted in Fig. 8.

Fig. 8
regular and irregular local geometries with a width of 30 mm: a,b) circular and polygonal shape at crack tip, c,d) circular and polygonal shape along the crack path

The von Mises stress contours after convergence are plotted on deformed shapes in Fig. 9 and Fig. 10. In these figures, the “a” and “b” are related to M1 and M2, “c” and “d” are related to M3 and M4, “e” is related to M5(a) and M6(a), and “f” is related to M5(b) and M6(b). As it is shown, the global and local displacements and stresses are quite consistent together. Concerning regular and irregular local boundaries (see Fig. 8b and 8d), the irregular one is less than 1% faster, due to its fully matched interfacial nodes and absence of MLS approximation. So from now on, we regard the regular one due to its simplicity.

The performance evaluation of the proposed coupling scheme is examined by the relative SIF error estimation. These factors are sensitive to global/local interfacial boundary conditions, and one can accurately monitor the performance by computing them. Therefore, the relative errors with respect to the reference analytical solution are considered. In this regard, these relative errors are depicted in Fig. 11, for different models at each iteration.

Fig. 9
von Mises stresses on global domain

Fig. 10
von Mises stresses on local domain

Fig. 11
SIF mode I convergence for various models

It seems that the accelerated models (with subscripts) are converged very quickly after three or four iterations with an error of 10−3. The procrustean opening of the global domain associated with the crack opening is the origin of several difficulties in M5 and M6 models, which is the greatest of the challenges in these two models in terms of localization, instability, and convergence.

Fig. 12 shows the SIFs (mode I and II) which is calculated with different rj/a (a = crack length) in M5 model. From this figure, same as the previous example, in general aspect of view, it can be seen that, the results have an acceptable accuracy in range of renrrj0.8rL, compared with EFG and XFEM.

Fig. 12
SIFs (mode I and II) for various rJ/a values in M5

To compare the coupling speed, we considered the local domain size as the most time-consuming part. In this regard, the M1 results with circular local domain are illustrated in Fig. 13.

Fig. 13
Iteration time and numbers to achieve convergence in M1 model, for each local area.

This figure shows the Comparison of the cumulative solution time for various local domain sizes. From this figure, one can observe that, the excessively small local domains could be a source of inaccuracy. In contrast, the bigger one would cause more iteration numbers, while it does not necessarily improve the accuracy. Therefore, the larger local domain causes an increase not only in computational time for each iteration but also in iteration numbers.

Finally, the computational speed and accuracy for all models, after first convergence (3 or 4 iterations), are shown in Fig. 14. It seems that, among the proposed models, M1 to M6, the most effective one is the M1 model, which is computationally fast and accurate, compared with EFG and XFEM.

From Fig. 14, one can observe that both accelerators had excellent performance in all situations, nevertheless, the SR1 was slightly more efficient but a little more time-consuming in comparison with Aitken Δ2.

Fig. 14
Comparison of the precision and speed (local area width is 3 cm)

4. Conclusion

In this paper, we proposed a global/local non-intrusive coupling approach to simulate a local discontinuity in a global linear structure. Since the global model is considered a plain model, it could not directly simulate the crack properties and its analysis. Therefore, this part of the analysis is performed by the more accurate and more specialized local model, overlapped over that global model, and ran through the IGL non-intrusive framework. In this technique, thanks to the non-intrusive approach, neither the connectivity nor the global solver are changed during the analysis. Therefore, various industrial FEA software and individual research codes can be linked smoothly.

In this research, the GFV and EFG methods are chosen respectively as global and local models due to their remarkable advantages. Six separate models (M1 to M6) are considered based on global/local configuration. Each of these models is analyzed with regular and irregular local boundary edges. Performance is improved by applying the MLS approximation on non-conforming (nested) interface between global coarse mesh and local fine nodes. Since the convergence rate is relatively slow, especially in fully non-intrusive models (M5 and M6), the efficient acceleration techniques based on Aitken’s ∆2 and Quasi-Newtonian SR1 methods are employed to improve the convergence speed. We also proposed a novel accelerator through the material constitutive matrix based upon the current SR1 method, intended to explicit global solvers.

To summarize the achievements of proposed GFV/EFG method for 2D modeling of local discontinuity, the following main conclusions can be stated:

  • ✓ The algorithm is robust and reliable because the accuracy can directly be increased as much as required by introducing the desired extent for the residual norms.

  • ✓ Considering the efficiency of the accelerators, both accelerators had excellent performance in all situations, nevertheless, the Quasi-Newtonian SR1 appears to be slightly more efficient but a little more time-consuming in comparison with Aitken Δ2. It seems that one may acquire a more reliable performance by applying the mixed acceleration.

  • ✓ The global data sets are remained intact during the analysis, which cause the method to be considerably non-intrusive, regarding the overlapping sub-domains act like the black-boxes. Therefore, the commercial software could be easily implemented during this process.

  • ✓ The algorithm is light and easy to implement, because it only requires to send and receive the supplementary interfacial nodal displacement or force residuals, which are conventional numerical outputs, between global and local solvers.

  • ✓ The SIFs are investigated for various J-integral radiuses, crack length and local area diameters. It was concluded that the results were considerably more accurate in the range of renrrj0.8rL.

  • ✓ In M5 and M6, due to the weak global/local geometric consistency, the crack path is successfully simulated in global domains by initially implementing the idea of zero value for modulus of elasticity on elements intersected with crack.

  • ✓ Due to the compatibility of the MLS operator, the data transition between two non-matching sub-domains is effectively handled, which supports the flexible modeling for designers.

    Nomenclature
  • A k   Area of triangular element
  • Δ l k   Normal vector of the k th surface in GFV
  • b i   Body force vector
  • ρ   Material mass density
  • E   Young’s modulus
  • υ   Poisson’s ratio
  • F i   Force vector
  • ΩG , Ω L   Global and local domain
  • Ω G|L   Intersection of global and local domain
  • Γ t , Γ u   Local natural and essential boundary edges
  • C w a v e   Wave velocity
  • P, R   Prolongation and restriction operators
  • C   Material constitutive matrix (D -1)
  • φ i   MLS shape function associated to node i
  • H ( x )   Heaviside step function
  • Q k ( x )   Enrichment function
  • KI , K II   Stress intensity factors (mode I and II)
  • uL , u G   Local and global displacement values
  • σL , σ G   Local and global stress values
  • W s   Strain energy density
  • δ i j   Kronecker delta function
  • F G | L   Interfacial nodal forces
  • u G | L   Interfacial nodal displacements
  • μ ( k )   Relaxation coefficient in the kth step
  • K L   Local stiffness matrix
  • B   Shape function derivatives matrix.
  • Η   L2-norm error
  • d max   Dimensionless size of the support domain
  • Re s k   The interfacial residual norm
  • r j   J-integral radius
  • r s   Support domain radius

References

  • Alkhamis, M. T. et al. (2008) ‘Utilizing NASIR galerkin finite volume analyzer for 2D plane strain problems under static and vibrating concentrated loads’, Jordan Journal of Civil Engineering, 2(4), pp. 335-343.
  • Allix, O. et al. (2011) ‘Non-intrusive coupling: An attempt to merge industrial and research software capabilities’, Recent Developments and Innovative Applications in Computational Mechanics, pp. 125-133. doi: 10.1007/978-3-642-17484-1_15.
    » https://doi.org/10.1007/978-3-642-17484-1_15
  • Amraei, A. and Fallah, N. (2016) ‘A development in the finite volume method for the crack growth analysis without global remeshing’, International Journal of Engineering, Transactions A: Basics, 29(7), pp. 890-900. doi: 10.5829/idosi.ije.2016.29.07a.03.
    » https://doi.org/10.5829/idosi.ije.2016.29.07a.03
  • Bauman, P. T. et al. (2008) ‘On the application of the Arlequin method to the coupling of particle and continuum models’, Computational Mechanics, 42(4), pp. 511-530. doi: 10.1007/s00466-008-0291-1.
    » https://doi.org/10.1007/s00466-008-0291-1
  • Belytschko, T., Fish, J. and Bayliss, A. (1990) ‘The Spectral Overlay on finite elements for problems with high gradients’, 81, pp. 71-89.
  • Belytschko, T., Lu, Y. Y. and Gu, L. (1994) ‘Element-Free Galerkin methods’, International Journal, 37(April 1993), pp. 229-256. doi: 10.1002/nme.1620370205.
    » https://doi.org/10.1002/nme.1620370205
  • Belytschko, T., Lu, Y. Y. and Gu, L. (1995) ‘Crack propagation by element-free Galerkin methods’, Engineering Fracture Mechanics, 51(2), pp. 295-315. doi: 10.1016/0013-7944(94)00153-9.
    » https://doi.org/10.1016/0013-7944(94)00153-9
  • Ben Dhia, H. (1998) ‘Multiscale mechanical problems: the Arlequin method’, Comptes Rendus De L Academie Des Sciences Serie Ii Fascicule B-Mecanique Physique Astronomie, 326(12), pp. 899-904. doi: doi:10.1016/S1251-8069(99)80046-5.
    » https://doi.org/doi:10.1016/S1251-8069(99)80046-5
  • Bouclier, R., Passieux, J. C. and Salaün, M. (2017) ‘Development of a new, more regular, mortar method for the coupling of NURBS subdomains within a NURBS patch: Application to a non-intrusive local enrichment of NURBS patches’, Computer Methods in Applied Mechanics and Engineering, 316, pp. 123-150. doi: 10.1016/j.cma.2016.05.037.
    » https://doi.org/10.1016/j.cma.2016.05.037
  • Chahine, E., Laborde, P. and Renard, Y. (2011) ‘A non-conformal eXtended Finite Element approach: Integral matching Xfem’, Applied Numerical Mathematics, 61(3), pp. 322-343. doi: 10.1016/j.apnum.2010.10.009.
    » https://doi.org/10.1016/j.apnum.2010.10.009
  • Duval, M. et al. (2014) ‘Local/Global Non-Intrusive Parallel Coupling for Large Scale Mechanical Analysis’, 343(July), pp. 1210-1211. doi: 10.1126/science.1249625.
    » https://doi.org/10.1126/science.1249625
  • Duval, M. et al. (2016a) ‘Non-intrusive Coupling: Recent Advances and Scalable Nonlinear Domain Decomposition’, Archives of Computational Methods in Engineering, 23(1), pp. 17-38. doi: 10.1007/s11831-014-9132-x.
    » https://doi.org/10.1007/s11831-014-9132-x
  • Duval, M. et al. (2016b) ‘Non-intrusive domain decomposition algorithm for solving nonlinear problems’, p. 3687.
  • El-Gebeily, M., Elleithy, W. M. and Al-Gahtani, H. J. (2002) ‘Convergence of the domain decomposition finite element-boundary element coupling methods’, Computer Methods in Applied Mechanics and Engineering, 191(43), pp. 4851-4867. doi: 10.1016/S0045-7825(02)00405-X.
    » https://doi.org/10.1016/S0045-7825(02)00405-X
  • Elleithy, W. M. and Al-Gahtani, H. J. (2000) ‘An overlapping domain decomposition approach for coupling the finite and boundary element methods’, Engineering Analysis with Boundary Elements, 24(5), pp. 391-398. doi: 10.1016/S0955-7997(00)00014-X.
    » https://doi.org/10.1016/S0955-7997(00)00014-X
  • Elleithy, W. M. and Tanaka, M. (2002) ‘Interface relaxation algorithms for coupling the FEM and BEM’, in International Series on Advances in Boundary Elements, pp. 721-730. doi: 10.2495/BE020671.
    » https://doi.org/10.2495/BE020671
  • Elleithy, W. M. and Tanaka, M. (2003) ‘Interface relaxation algorithms for BEM-BEM coupling and FEM-BEM coupling’, Computer Methods in Applied Mechanics and Engineering, 192(26-27), pp. 2977-2992. doi: 10.1016/S0045-7825(03)00312-8.
    » https://doi.org/10.1016/S0045-7825(03)00312-8
  • Estorff, O. Von and Antes, H. (1991) ‘On FEM-BEM coupling for fluid-structure interaction analyses in the time domain’, International Journal for Numerical Methods in Engineering, 31(6), pp. 1151-1168. doi: 10.1002/nme.1620310609.
    » https://doi.org/10.1002/nme.1620310609
  • Fish, J. (1993) ‘Adaptive S-method for Linear Elastostatics’, 104, pp. 363-396.
  • Fries, T. P. and Belytschko, T. (2010) ‘The extended / generalized finite element method : An overview of the method and its applications’, Review Literature And Arts Of The Americas, (April), pp. 253-304. doi: 10.1002/nme.
    » https://doi.org/10.1002/nme
  • Gendre, L. et al. (2009) ‘Non-intrusive and exact global / local techniques for structural problems with local plasticity’, 44, pp. 233-245. doi: 10.1007/s00466-009-0372-9.
    » https://doi.org/10.1007/s00466-009-0372-9
  • Ghorashi, S. S., Mohammadi, S. and Sabbagh-Yazdi, S. R. (2011) ‘Orthotropic enriched element free Galerkin method for fracture analysis of composites’, Engineering Fracture Mechanics. Elsevier Ltd, 78(9), pp. 1906-1927. doi: 10.1016/j.engfracmech.2011.03.011.
    » https://doi.org/10.1016/j.engfracmech.2011.03.011
  • Ghorashi, S. S., Sabbagh-Yazdi, S. R. and Mohammadi, S. (2010) ‘Element free Galerkin method for crack analysis of orthotropic plates’, 1(1), pp. 1-13.
  • Godinho, L. et al. (2013) ‘Iterative coupling between the MFS and Kansa’s method for acoustic problems’, WIT Transactions on Modelling and Simulation, 54(April 2016), pp. 123-132. doi: 10.2495/BEM130121.
    » https://doi.org/10.2495/BEM130121
  • Gupta, V. et al. (2015) ‘Stable GFEM (SGFEM): Improved conditioning and accuracy of GFEM/XFEM for three-dimensional fracture mechanics’, Computer Methods in Applied Mechanics and Engineering, 289, pp. 355-386. doi: 10.1016/j.cma.2015.01.014.
    » https://doi.org/10.1016/j.cma.2015.01.014
  • Henshaw, W. D. (2005) ‘On Multigrid for Overlapping Grids’, SIAM Journal on Scientific Computing, 26(5), pp. 1547-1572. doi: 10.1137/040603735.
    » https://doi.org/10.1137/040603735
  • Hirai, I., Wang, B. P. and Pilkey, W. D. (1984) ‘An efficient zooming method for finite element analysis’, International Journal for Numerical Methods in Engineering, 20(9), pp. 1671-1683. doi: 10.1002/nme.1620200910.
    » https://doi.org/10.1002/nme.1620200910
  • Holl, M., Loehnert, S. and P. Wriggers (2012) ‘An adaptive multiscale method for crack propagation and crack coalescence’. doi: 10.1002/nme.
    » https://doi.org/10.1002/nme
  • Kim, D. J., Duarte, C. A. and Proenca, S. P. (2012) ‘A generalized finite element method with global-local enrichment functions for confined plasticity problems’, Computational Mechanics, 50(5), pp. 563-578. doi: 10.1007/s00466-012-0689-7.
    » https://doi.org/10.1007/s00466-012-0689-7
  • Kim, J.-H. and Paulino, G. H. (2002) ‘Finite element evaluation of mixed mode stress intensity factors in functionally graded materials’, International Journal for Numerical Methods in Engineering, 53(8), pp. 1903-1935. doi: 10.1002/nme.364.
    » https://doi.org/10.1002/nme.364
  • Küttler, U. and Wall, W. A. (2008) ‘Fixed-point fluid-structure interaction solvers with dynamic relaxation’, Computational Mechanics, 43(1), pp. 61-72. doi: 10.1007/s00466-008-0255-5.
    » https://doi.org/10.1007/s00466-008-0255-5
  • Lu, Y. Y., Belytschko, T. and Gu, L. (1994) ‘A new implementation of the element free Galerkin method’, Computer Methods in Applied Mechanics and Engineering, 113(3-4), pp. 397-414. doi: 10.1016/0045-7825(94)90056-6.
    » https://doi.org/10.1016/0045-7825(94)90056-6
  • Lu, Y. Y., Belytschko, T. and Liu, W. K. (1991) ‘a variational coupled FE-BE method for elasticity and fracture mechanics’, 85, pp. 21-37.
  • Nguyen-Xuan, H. et al. (2012) ‘A cell-based smoothed finite element method for three dimensional solid structures’, KSCE Journal of Civil Engineering.
  • Organ, D. et al. (1996) ‘Continuous meshless approximations for nonconvex bodies by diffraction and transparency’, Computational Mechanics, 18(3), pp. 225-235. doi: 10.1007/BF00369940.
    » https://doi.org/10.1007/BF00369940
  • Park, K. C. and Felippa, C. A. (2000) ‘A variational principle for the formulation of partitioned structural systems’, International Journal for Numerical Methods in Engineering, 47(1-3), pp. 395-418. doi: 10.1002/(SICI)1097-0207(20000110/30)47:1/3<395::AID-NME777>3.0.CO;2-9.
    » https://doi.org/10.1002/(SICI)1097-0207(20000110/30)47:1/3<395::AID-NME777>3.0.CO;2-9
  • Passieux, J. et al. (2013) ‘Local / global non-intrusive crack propagation simulation using a multigrid X-FEM solver’, Computational Mechanics, Springer Verlag, 52(6), pp. 1381-1393.
  • Rangarajan, R. and Lew, A. J. (2012) ‘A local multigrid X-FEM strategy for 3-D crack propagation’, (August 2008), pp. 581-600. doi: 10.1002/nme.
    » https://doi.org/10.1002/nme
  • Réthoré, J., Roux, S. and Hild, F. (2012) ‘Hybrid analytical and extended finite element method (HAX-FEM): A new enrichment procedure for cracked solids’, International Journal for Numerical Methods in Engineering, 81(3), pp. 269-285. doi: 10.1002/nme.
    » https://doi.org/10.1002/nme
  • Sabbagh-Yazdi, S. R. and Amiri-SaadatAbadi, T. (2013) ‘GFV solution on UTE mesh for transient modeling of concrete aging effects on thermal plane strains during construction of gravity dam’, Applied Mathematical Modelling, 37(1-2), pp. 82-101. doi: 10.1016/j.apm.2012.01.045.
    » https://doi.org/10.1016/j.apm.2012.01.045
  • Sabbagh-Yazdi, S. R. and Bayatlou, M. (2012) ‘Equilibrium condition nonlinear modeling of a cracked concrete beam using a 2D Galerkin finite volume solver’, Meth. Civil Eng, 3, pp. 63-76.
  • Sabbagh-Yazdi, S. R. and Najar-Nobari, H. (2020) ‘Coupling Nonlinear Element Free Galerkin and Linear Galerkin Finite Volume Solver for 2D Modeling of Local Plasticity in Structural Material’, International Journal of Engineering, 33(3), pp. 387-400. doi: 10.5829/ije.2020.33.03c.03.
    » https://doi.org/10.5829/ije.2020.33.03c.03
  • Sabbagh-Yazdi, S. R., Ali-Mohammadi, S. and Pipelzadeh, M. K. (2012) ‘Unstructured finite volume method for matrix free explicit solution of stress-strain fields in two dimensional problems with curved boundaries in equilibrium condition’, Applied Mathematical Modelling, 36(5), pp. 2224-2236. doi: 10.1016/j.apm.2011.08.001.
    » https://doi.org/10.1016/j.apm.2011.08.001
  • Sabbagh-Yazdi, S. R., Amiri, T. and Gharebaghi, S. A. (2019) ‘A Proposed Damping Coefficient of Quick Adaptive Galerkin Finite Volume Solver for Elasticity Problems’, 13(1), pp. 56-79. doi: 10.24874/jsscm.2019.13.01.04.
    » https://doi.org/10.24874/jsscm.2019.13.01.04
  • Sabbagh-Yazdi, S. R., Farhoud, A. and Asil Gharebaghi, S. (2018) ‘Simulation of 2D linear crack growth under constant load using GFVM and two-point displacement extrapolation method’, Applied Mathematical Modelling. Elsevier Inc., 61, pp. 650-667. doi: 10.1016/j.apm.2018.05.022.
    » https://doi.org/10.1016/j.apm.2018.05.022
  • Soares Jr, D., Estorff, O. von and Mansur, W. J. (2004) ‘Iterative coupling of BEM and FEM for nonlinear dynamic analyses’, Computational Mechanics, 34(1), pp. 67-73. doi: 10.1007/s00466-004-0554-4.
    » https://doi.org/10.1007/s00466-004-0554-4
  • Soares, D. (2011) ‘Coupled numerical methods to analyze interacting acoustic-dynamic models by multidomain decomposition techniques’, Mathematical Problems in Engineering, 2011. doi: 10.1155/2011/245170.
    » https://doi.org/10.1155/2011/245170
  • Tada, H., Paris, P. C. and Irwin, G. (1973) ‘The Stress Analysis of Cracks Handbook’, Del Research Corporation, Hellertown, PA.
  • von Estorff, O. and Hagen, C. (2006) ‘Iterative coupling of FEM and BEM in 3D transient elastodynamics’, Engineering Analysis with Boundary Elements, 30(7), pp. 611-622. doi: 10.1016/j.enganabound.2006.01.007.
    » https://doi.org/10.1016/j.enganabound.2006.01.007
  • Whitcomb, J. D. (1991) ‘Iterative global/local finite element analysis’, Computers and Structures, 40(4), pp. 1027-1031. doi: 10.1016/0045-7949(91)90334-I.
    » https://doi.org/10.1016/0045-7949(91)90334-I
  • Whitcomb, J. D. and Woo, K. (1993a) ‘Application of Iterative Global/Local Finite Element Analysis. Part 1: Linear Analysis’, Communications in Numerical Methods in Engineering, 9(9), pp. 745-756. doi: 10.1002/cnm.1640090905.
    » https://doi.org/10.1002/cnm.1640090905
  • Whitcomb, J. D. and Woo, K. (1993b) ‘Application of Iterative Global/Local Finite Element Analysis. Part 2: Geometrically Non-Linear Analysis’, 9(December 1991), pp. 757-766.
  • Zienkiewicz, O. C., Kelly, D. W. and Bettess, P. (1977) ‘The coupling of the finite element method and boundary solution procedures’, International Journal for Numerical Methods in Engineering, 11(2), pp. 355-375. doi: 10.1002/nme.1620110210.
    » https://doi.org/10.1002/nme.1620110210

Edited by

Editor:

Marco L. Bittencourt.

Publication Dates

  • Publication in this collection
    16 Oct 2020
  • Date of issue
    2020

History

  • Received
    01 Aug 2020
  • Reviewed
    11 Aug 2020
  • Accepted
    31 Aug 2020
  • Published
    08 Sept 2020
Individual owner www.lajss.org - São Paulo - SP - Brazil
E-mail: lajsssecretary@gmsie.usp.br