A novel X-FEM based fast computational method for crack propagation
11 A novel X-FEM based fast computational method for crack propagation
Zhenxing Cheng a , Hu Wang a* a State Key Laboratory of Advanced Design and Manufacturing for Vehicle Body, Hunan University, Changsha, 410082, P.R. China
Abstract
This study suggests a fast computational method for crack propagation, which is based on the extended finite element method (X-FEM). It is well known that the X-FEM might be the most popular numerical method for crack propagation. However, with the increase of complexity of the given problem, the size of FE model and the number of iterative steps are increased correspondingly. To improve the efficiency of X-FEM, an efficient computational method termed decomposed updating reanalysis (DUR) method is suggested. For most of X-FEM simulation procedures, the change of each iterative step is small and it will only lead a local change of stiffness matrix. Therefore, the DUR method is proposed to predict the modified response by only calculating the changed part of equilibrium equations.
Compared with other fast computational methods, the distinctive characteristic of the proposed method is to update the modified stiffness matrix with a local updating strategy, which only the changed part of stiffness matrix needs to be updated. To verify the performance of the DUR method, several typical numerical examples have been analyzed and the results demonstrate that this method is a highly efficient method with high accuracy.
Keywords
Crack propagation, Fast computational method, Decomposed updating reanalysis method, Extended finite element method * Corresponding author Tel: +86 0731 88655012; fax: +86 0731 88822051. E-mail address: [email protected] (Hu Wang)
A great amount of engineering practice indicates that the quality and stability of engineering structures are closely related to the internal crack propagation[1, 2]. Therefore, prediction of the path of crack propagation and analysis of the stability of crack are significant for estimating the safety and the reliability of engineering structures. There are many numerical methods have been developed to simulate crack propagation process , such as finite element method (FEM) [3], boundary element method (BEM) [4], meshless method [5, 6], edge-based finite element method (ES-FEM) [7, 8], numerical manifold method (NMM) [9, 10], extended finite element method (X-FEM) [11-13] and so on[14]. Compared with above methods, the X-FEM might be the most popular numerical method for the crack propagation simulation due to its superiority of modeling both strong and weak discontinuities within a standard FE framework. The X-FEM was first developed by Belytschko and Black [15]. They analyzed the crack propagation problem with minimal re-meshing. Then, Dolbow et al [16]. and Moës et al. [17] improved this method by using a Heaviside function to enrichment function, and it also has been extended to 3D static crack modeling by Sukumar et al.[18]. Sequentially, the X-FEM was significantly improved by coupling with the level set method (LSM) which is used to track both the crack position and tips [19]. Moreover, the X-FEM has been applied to multiple engineering fields, such as dynamic crack propagation or branching [20, 21], crack propagation in composites [22] or shells [23-25], multi-field problems [26], multi-material problems [27, 28], solidification [29], shear bands [30], dislocations [31] and so on. More details of the development of X-FEM can be found in Refs [12, 32-34]. Generally, in order to improve the accuracy of simulation, a very refined mesh with a very small increment of crack propagation or fatigue cycles should be engaged. Correspondingly, the computational cost is expensive. Therefore, reanalysis algorithm is used to improve the efficiency. Reanalysis, as a fast computational method, is used to predict the response of modified structures efficiently without full analysis, and reanalysis method can be divided into two categories: direct methods (DMs) and approximate methods. DMs can update the inverse of modified stiffness matrix quickly by Sherman-Morrison-Woodbury lemma [35, 36] and obtain the exact response of the modified structure, but usually it can only solve the problems of local or low-rank modifications. In recent decades, many DMs have been suggested. Such as, Song et al. suggested a novel direct reanalysis algorithm based on the binary tree characteristic to update the triangular factorization in sparse matrix solution [37]. Liu et al. applied Cholesky factorization to structural reanalysis [38]. Huang and Wang suggested an independent coefficient (IC) method for large-scale problems with local modification [39]. Compared with DMs, approximate methods can solve the high-rank modifications, but the exact response usually cannot be obtained. The approximate methods mainly include local approximations (LA), global approximations (GA) [40], iterative approximations (IA) [41] and combined approximations (CA) [42-46]. Moreover, many other reanalysis methods have been proposed in recent years. Zuo et al. combined reanalysis method with genetic algorithm (GA) [47]. Sun et al . extended the reanalysis method into a structural optimization process [48]. To improve the efficiency of reanalysis method, He et al. developed a multiple-GPU based parallel IC reanalysis method [49]. Materna et al. applied the reanalysis method to nonlinear problems [50]. In this study, a fast computational method is proposed to model crack propagation under the framework of X-FEM. It is observed that modeling crack propagation by the X-FEM will bring the additional DOFs in each iterative step and it will lead to a local change of stiffness matrix. Considering this characteristic, the decomposed updating reanalysis (DUR) method theoretically can improve the efficiency of X-FEM significantly. Moreover, a local updating stiffness matrix strategy is suggested to improve the efficiency of stiffness matrix assembling. Furthermore, in order to guarantee the accuracy and efficiency of the DUR method, a local strategy has been introduced to update the Cholesky factorization of stiffness matrix efficiently.
The rest of this paper is organized as follows. The basic theories of X-FEM are briefly introduced in Section 2. The details of DUR method are described in Section 3. Then, several numerical examples are tested in Section 4 to investigate the performance of the DUR method. Finally, conclusions are summarized in Section 5.
In the X-FEM, the standard FEM shape function should be enriched by the enrichment function. Assume that the enriched displacement approximation of X-FEM can be defined as: tan ( ) ( ) ( ) ( ) E h I I J JI J s dard enrich N N u u u x x u x x q , (1) where I N and I u denote the standard FEM shape function and nodal degrees of freedom (DOF), respectively. The enriched displacement approximation should be divided into two parts: the standard finite element approximation and partition of unity enrichment. The x ( ) means enrichment function while the J q is the additional nodal degrees of freedom. Fig. 1 An arbitrary crack line in a structured mesh Consider an arbitrary crack in a structured mesh as shown in Fig. 1, then Eq.(1) can be rewritten as ( ) ( ) ( ) ( ) ( ) ( ) S T h I I I I I I I II I I
N H N N u x x u x x a x x b , (2) where is the solution domain, S is the domain cut by crack, T is the domain which crack tip located, ( ) H x is the shifted Heaviside enrichment and ( ) x is the shifted crack tip enrichment. The details of ( ) H x and ( ) x are given as following:
1( ) 1
Above crackH Below crack x , (3) { ( )} sin , cos , sin sin , sin cos2 2 2 2 x r . (4) The discrete X-FEM equations are obtained by substituting Eq.(2) into the principle of virtual work. Assume that the discrete equations can be defined as uu ua ub uTua aa ab aT Tub ab bb b K K K u FK K K FK K K b F a , (5) where uu K is the traditional finite element stiffness matrix, , , ua aa ab K K K are components with Heaviside enrichment and , , ub ab bb
K K K are components with crack tip enrichment.
Generally, the direction and magnitude of crack propagation at each iterative step are used to determine how the crack will propagate. The direction of crack propagation is found from the maximum circumferential stress criterion and the crack will propagate in the direction where is maximum [51]. The angle of crack propagation is defined as
12 arctan 84
I IIIII II
K Ksign KK K , (6) where is defined in the crack tip coordinate system, I K and II K are the mixed-mode stress intensity factors. The details are given in the reference [51]. There are two main quasi-static manners when modeling crack growth. The first one assumes a constant increment of crack growth at each cycle [16]while the other option is to assume a constant number of cycles and apply a fatigue crack growth law to predict the crack growth increment for the fixed number of cycles [52]. In this study, a fixed increment of crack growth a is considered. The DUR method is proposed to model quasi-static crack propagation under the framework of X-FEM and the framework of the DUR method is presented in Fig. 2.
Start
EndCrack stop extending?Separate equilibrium equationsInitial equilibrium equation solution
Calculate the coefficients vector y Calculate modified displacements:
Local stiffness matrix updating strategyLocal Cholesky factorization updating strategy
Extract unbalanced equations
Reanalysis algorithm B a l a n ce d b l o c k Unbalanced block
NoYes
Yes No Yes
NoExtract balanced equationsCalculate fundamental solution system B Fig. 2 The framework of the DUR method It can be found that the DUR method includes three parts: local updating stiffness matrix strategy, local Cholesky factorization updating strategy and reanalysis algorithm. The local stiffness matrix updating strategy is suggested to improve the efficiency of stiffness matrix assembling according to the characteristic of X-FEM. Moreover, considering the local change of stiffness matrix, the specific reanalysis algorithm is used to predict the modified response and improve the efficiency of the X-FEM. Furthermore, in order to guarantee the accuracy and the efficiency of the suggested reanalysis method, a local strategy has been introduced to update the Cholesky factorization of stiffness matrix efficiently. More details can be found in Section 3.2, 3.3 and 3.4.
As mentioned above, the stiffness matrix can be given as the following form, and each component has been marked by different colors associated with Fig. 1. Fig. 3 Color-marked stiffness matrix Consider a small increment of crack based on Fig. 1, and the comparison is shown in Fig. 4. Fig. 4 A small increment of crack propagation It can be found that the stiffness component uu K associated with the standard finite element approximation will be constant at each iterative step of crack propagation. It means that the change part of stiffness matrix is the enriched part, which should be much smaller than the un-enriched part. Furthermore, it can be found that once an element had become a split element, it will be a split element during entire computational procedure. It implies that once an element has been enriched by Heaviside function, the value of stiffness matrix will be a constant in the subsequent iterative steps. Therefore, the changed part of the stiffness matrix in each iterative step will be a small part as shown in Fig. 5. Fig. 5 Changed and unchanged part of stiffness matrix Assume that the unchanged part is a
N A N A matrix and the changed part contains a
B B matrix with N A B , B N A coupling matrices. If
N A B , only a small part of stiffness matrix needs to be updated in each iterative steps. Obviously, it should be high efficient and the matrix assembling cost should be significantly reduced.
The reanalysis algorithm is used to predict the response of the current iterative step by using the information of the first iterative step. In this study, a reanalysis algorithm is proposed for X-FEM according to the characteristic of X-FEM. The reanalysis algorithm avoids the full analysis after the first iterative step, and the response of the subsequent iterative steps can be efficiently obtained. The details of the reanalysis algorithm are described as following. Assume that the equilibrium equation of the i-th iterative step is i i i K U F , (7) where i U is the displacement in the i-th iterative step, and the equilibrium equation of the first iterative step can be given as K U F . (8) Assume that the solution of the equation in the i-th iterative step can be defined as i i U U U , (9) then substitute Eq.(9) into Eq.(7),0
The reanalysis algorithm is used to predict the response of the current iterative step by using the information of the first iterative step. In this study, a reanalysis algorithm is proposed for X-FEM according to the characteristic of X-FEM. The reanalysis algorithm avoids the full analysis after the first iterative step, and the response of the subsequent iterative steps can be efficiently obtained. The details of the reanalysis algorithm are described as following. Assume that the equilibrium equation of the i-th iterative step is i i i K U F , (7) where i U is the displacement in the i-th iterative step, and the equilibrium equation of the first iterative step can be given as K U F . (8) Assume that the solution of the equation in the i-th iterative step can be defined as i i U U U , (9) then substitute Eq.(9) into Eq.(7),0 i i i K U U F , (10) or i i i i K U F K U . (11) Define the residual value of displacement δ as i i i δ F K U , (12) then Eq.(11) can be written as i K U δ . (13) Consider that only a small part of stiffness matrix will change in every iterative step, the most part of δ should be zero. Based on this property, the i U can be divided into two blocks: unbalanced and balanced blocks, according to Eq.(14): i sum Δ K K δ . (14) If j Δ , the j-th DOF is unbalanced, otherwise the j-th
DOF is balanced. Accord to this, the Eq.(13) can be rewritten as i imm mn mi i n nnm nn
K K U 0U δK K , (15) where m is the number of balanced DOFs, and n is the number of unbalanced DOFs. Equation (15) can be rewritten as i imm m mn n K U K U 0 (16) and i inm m nn n n
K U K U δ . (17) Obviously, Eq.(16) is a homogeneous equation set which has infinite solutions, and the solution U is one of them. Let Tn U (where 1 is the k-th member of n U , and
1, 2, , k n ), then the fundamental solution system of Eq. (16) can be obtained as i imm mnnn K KB E , (18) where nn E is a rank- n unit matrix. Define the general solution of Eq.(16) is U By , (19) where y is a dimension- n vector. Then, substitute Eq.(19) into Eq.(17) to find a unique solution of Eq.(13), obtain i i i inn nm mm mn n K K K K y δ . (20) Solve Eq.(20), y can be obtained, and then U can be obtained by Eq.(19). Sequentially, the i U can be obtained by Eq.(9). Because only a small part of stiffness matrix will change in each iterative step, the number n should not be too large. Therefore, Eq.(20) is a small scale problem. The key point is how to obtain imm K from the K . Assume that mm mnnm nn K KK K K , (21) and then the Cholesky factorization of K can be defined as T TmmT mm nmTnm nn nnT Tmm mm mm nmT T Tnm mm nm nm nn nn
L 0 L LK LL L L 0 LL L L LL L L L L L , (22) where mm L and nn L are lower triangular matrices. Compared with Eq.(21), it can be found that i Tmm mm mm mm K K L L . (23) Therefore, the fundamental solution system B can be calculated by Eq.(18), and the Cholesky factorization of stiffness matrix in the first iterative step can directly re-used. In describe this method clearly, the imm K can be associated with Fig. 4. Assume that the left figure of Fig. 4 shows the crack position in the first iterative step, and the right figure of Fig. 4 shows the crack position in the second iterative step. Then the stiffness matrices of X-FEM in the first and the second iterative steps can be written as Fig. 6 and Fig. 7 respectively. Fig. 6 The stiffness matrix of the first iterative step Fig. 7 The stiffness matrix of the second iterative step It can be found that the red-marked part is constant as shown in Fig. 6 and Fig. 7. Compared with Eq.(15), the mm K can be give as uu uamm Tua aa K KK K K , (24) and the mn K , nm K and nn K can be give as ua ubmn K KK 0 0 , (25)
22 2
Tuanm Tub
K 0K K 0 , (26) aa abnn Tab bb
K KK K K . (27) Generally, the accuracy and the efficiency of the reanalysis methods are depended on the percentage of changed part, and if the percentage of changed part is too large, the accuracy and the efficiency should be unavailable for engineering problems. Therefore, a suitable critical value of changed percentage should be stipulated, and usually 5% is chosen for the critical value [53]. In the DUR method, the changed part is composed by unbalanced DOFs. Then the changed percentage can be defined as the following form: nN , (28) where n is the number of unbalanced DOFs and N is the number of total DOFs. In order to guarantee the accuracy and efficiency of the reanalysis method, a suitable critical value of should be set. In this study, 5% is also chose as the critical value of , so the initial information needs to be updating when the and the key issue is how to obtain the Cholesky factorization of the initial stiffness matrix efficiently. Therefore, a local updating Cholesky factorization strategy has been suggested. This strategy is based on the property that only a small part of stiffness matrix will change in each iterative step, so that the Cholesky factorization of modified stiffness matrix can be updating based on the Cholesky factorization of initial stiffness matrix. An example is given to explain this strategy, which is based on Fig. 6 and Fig. 7. Assume that Fig. 6 means the initial stiffness matrix and Fig. 7 means the modified stiffness matrix. Recall the Cholesky factorization of initial stiffness matrix takes the form T T TT T
K K L 0 L LK L LK K L L 0 L , (29) where / T cholchol L KL L KL K L L . (30) Compared with Fig. 6, K can be considered equivalent to uu uaTua aa K KK K K , (31) and K , K can be written as
11 1 112 221 , ub bbab KK K KK . (32) Then the Cholesky factorization of modified stiffness matrix L can be obtained by
22 112 221 22
L 0L L L , (33) where / T chol cholchol L K K LL L KL K L L , (34) and compared with Fig. 7, , aa abua ub Tab bb K KK KK K0 0 K K . (35) Then L can be reused at subsequent iterative steps of crack propagation, and usually the size of L is much large than other parts. Therefore, the strategy should save much computational cost than calculate the Cholesky factorization of the entire stiffness matrix directly. In order to test the accuracy and efficiency of the DUR method, three examples are tested by the proposed methods. These three cases involve edge and center crack propagation, concentrated and uniformed load problems, thus the performance of the DUR method could be verified thoroughly. In this study, the comparison has been made between the DUR and full analysis, and the errors of displacement, Von Mises stress and Von Mises strain are defined by the following formulas:
PFR FAu FA E U UU (36)
PFR FAFA E σ σσ (37) where PFR U , PFR σ mean the results of DUR method, and FA U , FA σ mean the results of full analysis. Moreover, in order to investigate the performance of the DUR method, the CPU running time has been recorded and all the simulations were performed on an Intel(R) Core(TM) i7-5820K 3.30GHz CPU with 32GB of memory within MATLAB R2016b in x64 Windows 7. The problem shown in Fig. 8 is an adaptation of an example presented in reference [54]. The initial crack length is a mm , the force F N , and linear elastic material behavior is assumed. The material is aluminum 7075-T6 with E MPa , and a plane strain state is considered. The increment of propagation a mm . For this problem, the result of experimental test with specimens 16 mm thick was given by Giner et al. [54] as shown in Fig. 9.
13 Ø13 . FF Ø2028.5
Fig. 8 The geometry of edge crack in a plate with a hole Fig. 9 The result of edge crack in a plate with a hole in the reference [54] Then solve this case by the DUR method and only the first iterative step need to be calculated by full analysis method while other iterative steps should be predicted efficiently by the DUR method. In order to investigate the accuracy of the DUR method, The comparisons of displacement and stress between the DUR and full analysis are illustrated in Fig. 10 and Fig. 11, respectively. It is obviously that the results of DUR method and full analysis are almost the same. Compared with Fig. 9, it is proved that the DUR method is accurate. Moreover, the errors of each iterative step which is defined by Eq.(36) and Eq.(37) are shown in Fig. 12. It can be found that the DUR method is accurate. The computational costs of the DUR and full analysis are also listed in Tab.
1. It shows that the computational cost of the DUR method is much cheaper than the full analysis method. Moreover, the computational results of crack tip coordinates are shown in Tab. 2 and the results of the DUR and full analysis are also the same. Fig. 10 Displacement contour of edge crack in a plate with a hole Fig. 11 Von Mises stress contour of edge crack in a plate with a hole Fig. 12 The errors of each iterative step Tab. 1 Performance comparison of edge crack in a plate with a hole
Computational Time/ s Average Errors DUR Full analysis Displacement Von Mises Stress 0.968 18.953 4.8916e-13 1.7356e-12
Tab. 2 Crack tip coordinates of edge crack in a plate with a hole
Iterative steps Crack tip coordinates of DUR Method Crack tip coordinates of full analysis X-coordinate Y-coordinate X-coordinate Y-coordinate 0 0.000 67.500 0.000 67.500 10.000 67.500 10.000 67.500 1 10.998 67.559 10.998 67.559 5 14.994 67.727 14.994 67.727 10 19.993 67.805 19.993 67.805 15 24.986 67.560 24.986 67.560 20 29.918 66.766 29.918 66.766 25 34.605 65.059 34.605 65.059 30 37.695 61.523 37.695 61.523
As shown in Fig. 13, a plate with a circular inclusion is considered. Assume that there is an edge crack in the left edge of the plate. The aluminum 7075-T6 is used as the material of plate in the previous case, the Young’s modulus E MPa , the Poisson’s ratio , and a plane strain state is assuming. The material of inclusion is considered as carbon fiber reinforced composite and assume that the Young’s modulus E MPa , the Poisson’s is the same as the plate. The initial crack length a mm , and the uniformed load
50 / q N mm . The increment of propagation a mm . Then the DUR method is used to solve this problem and only the first iterative step needs to be calculated by full analysis method while others iterative steps should be predicted efficiently by the DUR method. qq R15
Inclusion
Fig. 13 The geometry of edge crack in a plate with a circular inclusion The comparisons of displacement and stress are shown as Fig. 14 and Fig. 15, Fig. 15 presents the Von Mises stress results of iterative step 1, 20 and 40. It is obviously that the result of DUR method is very close to the full analysis method. Furthermore, an error list of some selected iterative steps are shown as Fig. 16. It can be found that the DUR method is highly accurate because the maximum of error is . Moreover, the computational costs of the DUR and full analysis are also listed in Tab. 3. It shows that the efficiency of the DUR method is much higher than the full analysis method. Moreover, the crack tip coordinates are also listed in Tab. 4 and the results of the DUR and full analysis are matched.0
Fig. 13 The geometry of edge crack in a plate with a circular inclusion The comparisons of displacement and stress are shown as Fig. 14 and Fig. 15, Fig. 15 presents the Von Mises stress results of iterative step 1, 20 and 40. It is obviously that the result of DUR method is very close to the full analysis method. Furthermore, an error list of some selected iterative steps are shown as Fig. 16. It can be found that the DUR method is highly accurate because the maximum of error is . Moreover, the computational costs of the DUR and full analysis are also listed in Tab. 3. It shows that the efficiency of the DUR method is much higher than the full analysis method. Moreover, the crack tip coordinates are also listed in Tab. 4 and the results of the DUR and full analysis are matched.0 Fig. 14 Displacement contour of edge crack in a plate with a circular inclusion Fig. 15 Von Mises stress contour of edge crack in a plate with a circular inclusion Fig. 16 The errors of each iterative step Tab. 3 Performance comparison of edge crack in a plate with a circular inclusion
Computational Time/ s Average Errors DUR Full analysis Displacement Von Mises Stress 4.078 17.203 1.0534e-11 1.0951e-11
Tab. 4 Crack tip coordinates of edge crack in a plate with a circular inclusion
Iterative steps Crack tip coordinates of DUR Method Crack tip coordinates of full analysis X-coordinate Y-coordinate X-coordinate Y-coordinate 0 0.000 60.000 0.000 60.000 10.000 60.000 10.000 60.000 1 10.987 59.842 10.987 59.842 5 14.941 59.240 14.941 59.240 10 19.828 58.188 19.828 58.188 15 24.647 56.857 24.647 56.857 20 29.401 55.308 29.401 55.308 25 34.135 53.698 34.135 53.698 30 38.956 52.383 38.956 52.383 35 43.911 51.754 43.911 51.754 40 48.908 51.594 48.908 51.594 45 53.908 51.552 53.908 51.552 49 55.908 51.545 55.908 51.545 A center crack in a plate with a circular inclusion and a hole is considered as shown in Fig. 17. The material of plate is aluminum 7075-T6, the Young’s modulus E MPa , the Poisson’s ratio , and a plane strain state is assuming. The material of inclusion is considered as carbon fiber reinforced composite and assume that the Young’s modulus E MPa , the Poisson’s is the same as the plate. The initial crack length a mm , and the uniformed load
50 / q N mm . The increment of propagation a mm . Then, the DUR method is used to solve this problem and only the first iterative step needs to be calculated by the full analysis method while other iterative steps should be predicted efficiently by the DUR method. qq R10 R10
Inclusion
Fig. 17 The geometry of center crack in a plate with a circular inclusion and a hole The displacement and Von Mises stress contours of iterative step 1, 10 and 24 are shown as Fig. 18 and Fig. 19. It is obviously that the result of the DUR method is very close to the full analysis method. Furthermore, an error list of some selected iterative steps are shown as Fig. 20. It can be found that the DUR method is accurate because the maximum of error is . Moreover, the computational costs of the DUR and full analysis are also listed in Tab. 5. It shows that the efficiency of the DUR method is much higher than full analysis method. Besides this, the crack tip coordinates are also listed in Tab. 6 and the result of DUR and full analysis are matched. Fig. 18 Displacement contour of center crack in a plate with a circular inclusion and a hole Fig. 19 Von Mises stress contour of center crack in a plate with a circular inclusion and a hole Fig. 20 The errors of each iterative step Tab. 5 Performance comparison of center crack in a plate with a circular inclusion and a hole
Computational Time/ s Average Errors DUR Full analysis Displacement Von Mises Stress 2.3438 8.4375 4.9236e-13 5.6032e-13
Tab. 6 Crack tip coordinates of center crack in a plate with a circular inclusion and a hole
Iterative steps Crack tip coordinates of DUR Method Crack tip coordinates of full analysis Crack tip 1 Crack tip 2 Crack tip 1 Crack tip 2
0 (25.000, 60.000) (35.000, 60.000) (25.000, 60.000) (35.000, 60.000) 1 (24.012,
Three numerical examples have been tested in this section and it can be found that the DUR is an accurate and efficient and method. Moreover, a representative case has been calculated by the DUR method under different computational scales from 1000 to 100,000 to fully investigate the efficiency of the DUR method. The log-log plots of comparison results are shown in Fig. 21, Fig. 22 and the error analysis is also shown. Fig. 21 The comparison of computational cost used in solving equilibrium equations Fig. 22 The comparison of computational cost used in stiffness matrix updating It can be found that the efficiency of the DUR method is much higher than the full analysis method, and the advantage is more obvious for large scale problems. It is obvious that the accuracy should be improved significantly with the increase of DOFs. Accord to Fig. 21, it can be observed that the accuracy of the DUR method is very high, and the DUR can be regarded as an exact method. Moreover, the computational cost of stiffness matrix updating was also plotted, and Fig. 22 shows that the local updating strategy largely reduces the computational cost of stiffness matrix updating than traditional global updating strategy.
In this study, the DUR method is proposed for crack propagation. The DUR method consists of three strategies: local stiffness matrix updating, decomposed updating reanalysis, local Cholesky factorization updating strategies. Considering the characteristic of local change of stiffness matrix during X-FEM iterative procedure, the local stiffness matrix updating method can achieve the modified stiffness matrix quickly, and the local Cholesky factorization updating strategy is used to guarantee the accuracy and efficiency of the DUR method. More importantly, the decomposed updating reanalysis strategy is suggested to improve the efficiency of solving the equilibrium equations significantly. Therefore, the DUR method not only reduces the computational cost of solving equilibrium equations but also saves computational cost of stiffness matrix assembling. Numerical examples show that the DUR method is accurate for the crack propagation. For both the edge and center crack propagation, the accuracy of DUR method is very high. The log-log plots show that the efficiency of DUR method is much more higher than full analysis, and the advantage is more obvious for large scale problems. Moreover, compared with other reanalysis method, the comparisons of stress between DUR and full analysis are made, and the stress can be obtained accurately and efficiently. Acknowledgements
This work has been supported by Project of the Key Program of National Natural Science Foundation of China under the Grant Numbers 11572120, National Key Research and Development Program of China 2017YFB0203701.
References [1] Kanaun S. Discrete model of hydraulic fracture crack propagation in homogeneous isotropic elastic media. International Journal of Engineering Science. 2017;110:1-14. [2] Areias P, Msekh MA, Rabczuk T. Damage and fracture algorithm using the screened Poisson equation and local remeshing. Engineering Fracture Mechanics. 2016;158:116-43. [3] Branco R, Antunes FV, Costa JD. A review on 3D-FE adaptive remeshing techniques for crack growth modelling. Engineering Fracture Mechanics. 2015;141:170-95. [4] Santana E, Portela A. Dual boundary element analysis of fatigue crack growth, interaction and linkup. Engineering Analysis with Boundary Elements. 2016;64:176-95. [5] Tanaka S, Suzuki H, Sadamoto S, Imachi M, Bui TQ. Analysis of cracked shear deformable plates by an effective meshfree plate formulation. Engineering Fracture Mechanics. 2015;144:142-57. [6] Khosravifard A, Hematiyan MR, Bui TQ, Do TV. Accurate and efficient analysis of stationary and propagating crack problems by meshless methods. Theoretical and Applied Fracture Mechanics. 2017;87:21-34. [7] Nguyen-Xuan H, Liu GR, Bordas S, Natarajan S, Rabczuk T. An adaptive singular ES-FEM for mechanics problems with singular field of arbitrary order. Computer Methods in Applied Mechanics and Engineering. 2013;253:252-73. [8] Liu GR, Nourbakhshnia N, Zhang YW. A novel singular ES-FEM method for simulating singular stress fields near the crack tips for linear fracture problems. Engineering Fracture Mechanics. 2011;78(6):863-76. [9] Liu Z, Zheng H, Sun C. A domain decomposition based method for two-dimensional linear elastic fractures. Engineering Analysis with Boundary Elements. 2016;66:34-48. [10] Zhang HH, Li LX, An XM, Ma GW. Numerical analysis of 2-D crack propagation problems using the numerical manifold method. Engineering Analysis with Boundary Elements. 2010;34(1):41-50. [11] Zeng Q, Liu Z, Xu D, Wang H, Zhuang Z. Modeling arbitrary crack propagation in coupled shell/solid structures with X-FEM. International Journal for Numerical Methods in Engineering. 2016;106(12):1018-40. [12] Belytschko T, Gracie R, Ventura G. A review of extended/generalized finite element methods for material modeling. Modelling and Simulation in Materials Science and Engineering. 2009;17(4):043001. [13] Chopp DL, Sukumar N. Fatigue crack propagation of multiple coplanar cracks with the coupled extended finite element/fast marching method. International Journal of Engineering Science. 2003;41(8):845-69. [14] Cuomo M, Contrafatto L, Greco L. A variational model based on isogeometric interpolation for the analysis of cracked bodies. International Journal of Engineering Science. 2014;80:173-88. [15] Belytschko T, Black T. Elastic crack growth in finite elements with minimal remeshing. International journal for numerical methods in engineering. 1999;45(5):601-20. [16] Dolbow J, Belytschko T. A finite element method for crack growth without remeshing. International journal for numerical methods in engineering. 1999;46(1):131-50. [17] Belytschko T, Moës N, Usui S, Parimi C. Arbitrary discontinuities in finite elements. International Journal for Numerical Methods in Engineering. 2001;50(4):993-1013. [18] Sukumar N, Moës N, Moran B, Belytschko T. Extended finite element method for three-dimensional crack modelling. International Journal for Numerical Methods in Engineering. 2000;48(11):1549-70. [19] Stolarska M, Chopp D, Moës N, Belytschko T. Modelling crack growth by level sets in the extended finite element method. International journal for numerical methods in Engineering. 2001;51(8):943-60. [20] Xu D, Liu Z, Liu X, Zeng Q, Zhuang Z. Modeling of dynamic crack branching by enhanced extended finite element method. Computational Mechanics. 2014;54(2):489-502. [21] Song J-H, Areias PMA, Belytschko T. A method for dynamic crack and shear band propagation with phantom nodes. International Journal for Numerical Methods in Engineering. 2006;67(6):868-93. [22] Huynh D, Belytschko T. The extended finite element method for fracture in composite materials. International Journal for Numerical Methods in Engineering. 2009;77(2):214-39. [23] Ahmed A, van der Meer FP, Sluys LJ. A geometrically nonlinear discontinuous solid-like shell element (DSLS) for thin shell structures. Computer Methods in Applied Mechanics and Engineering. 2012;201-204:191-207. [24] Zhuang Z, Cheng B. A novel enriched CB shell element method for simulating arbitrary crack growth in pipes. Science China Physics, Mechanics and Astronomy. 2011;54(8):1520-31. [25] Areias P, Belytschko T. Non‐linear analysis of shells with arbitrary evolving cracks using XFEM. International Journal for Numerical Methods in Engineering. 2005;62(3):384-415. [26] Zilian A, Legay A. The enriched space–time finite element method (EST) for simultaneous solution of fluid–structure interaction. International Journal for Numerical Methods in Biomedical Engineering. 2008;75(3). [27] Sukumar N, Chopp DL, Moës N, Belytschko T. Modeling holes and inclusions by level sets in the extended finite-element method. Computer methods in applied mechanics and engineering. 2001;190(46):6183-200. [28] Zhuang Z, Cheng BB. Equilibrium state of mode-I sub-interfacial crack growth in bi-materials. International Journal of Fracture. 2011;170(1):27-36. [29] Chessa J, Smolinski P, Belytschko T. The extended finite element method (XFEM) for solidification problems. International Journal for Numerical Methods in Engineering. 2002;53(8):1959-77. [30] Areias P, Belytschko T. Two‐scale shear band evolution by local partition of unity. International Journal for Numerical Methods in Engineering. 2006;66(5):878-910. [31] Belytschko T, Gracie R. On XFEM applications to dislocations and interfaces. International Journal of Plasticity. 2007;23(10):1721-38. [32] Yazid A, Abdelkader N, Abdelmadjid H. A state-of-the-art review of the X-FEM for computational fracture mechanics. Applied Mathematical Modelling. 2009;33(12):4269-82. [33] Abdelaziz Y, Hamouine A. A survey of the extended finite element. Computers & Structures. 2008;86(11-12):1141-51. [34] Fries TP, Belytschko T. The extended/generalized finite element method: an overview of the method and its applications. International Journal for Numerical Methods in Engineering. 2010;84(3):253-304. [35] Woodbury MA. Inverting modified matrices. Memorandum report. 1950;42(106):336. [36] Sherman J, Morrison WJ. Adjustment of an Inverse Matrix Corresponding to a Change in One Element of a Given Matrix. The Annals of Mathematical Statistics. 1950;21(1):124-7. [37] Song Q, Chen P, Sun S. An exact reanalysis algorithm for local non-topological high-rank structural modifications in finite element analysis. Computers & Structures. 2014;143:60-72. [38] Liu HF, Wu BS, Li ZG. Method of Updating the Cholesky Factorization for Structural Reanalysis with Added Degrees of Freedom. Journal of Engineering Mechanics. 2014;140(2):384-92. [39] Huang G, Wang H, Li G. A reanalysis method for local modification and the application in large-scale problems. Structural and Multidisciplinary Optimization. 2013;49(6):915-30. [40] Haftka RT, Nachlas JA, Watson LT, Rizzo T, Desai R. Two-point constraint approximation in structural optimization. Computer Methods in Applied Mechanics and Engineering. 1987;60(3):289-301. [41] Li Z, Wu B. A preconditioned conjugate gradient approach to structural reanalysis for general layout modifications. International journal for numerical methods in engineering. 2007;70(5):505-22. [42] Kirsch U. A unified reanalysis approach for structural analysis, design, and optimization. Structural and Multidisciplinary Optimization. 2003;25(2):67-85. [43] Chen SH, Yang ZJ. A universal method for structural static reanalysis of topological modifications. International Journal for Numerical Methods in Engineering. 2004;61(5):673-86. [44] Wu BS, Lim CW, Li ZG. A finite element algorithm for reanalysis of structures with added degrees of freedom. Finite Elements in Analysis and Design. 2004;40(13-14):1791-801.0