Abstract
The eXtended Finite Element Method (XFEM) is a versatile method for solving crack propagation problems. Meanwhile, XFEM predictions for crack onset and propagation rely on the stress field which tends to converge at a slower rate than that of displacements, making it challenging to capture critical load at crack onset accurately. Furthermore, identifying the critical region(s) for XFEM nodal enrichments is user-dependent. The identification process can be straightforward for small-scale test specimen while in other cases such as complex structures it can be unmanageable. In this work a novel approach is proposed with three major objectives; alleviate user-dependency; enhance predictions accuracy; minimize computational effort. An automatic critical region(s) identification based on material selected failure criterion is developed. Moreover, the approach enables the selection of optimized mesh necessary for accurate prediction of failure loads at crack initiation. Also, optimal enrichment zone size determination is automated. The proposed approach was developed as an iterative algorithm and implemented in ABAQUS using Python scripting. The proposed algorithm was validated against our test data of unnotched specimens and relevant test data from the literature. The results of the predicted loads/displacements at failure are in excellent agreement with measurements. Crack onset locations were in very good agreement with observations from testing. Finally, the proposed algorithm has shown a significant enhancement in the overall computational efficiency compared to the conventional XFEM. The proposed algorithm can be easily implemented into user-built or commercial finite element codes.
1. Introduction
Stress concentrations may cause the initiation of a fatigue crack in a structure [1]. It arises mainly by either a concentrated force acting on a body or a geometrical discontinuity such as holes, a sharp geometrical change, or a cracked surface. Finite Element Method (FEM) has proven to be a versatile analysis technique in solving structural engineering problems. Meanwhile, problems involving discontinuities or singularities, for instance, crack onset and propagation problems, are quite problematic [2]. One major challenge in the conventional FEM is the need for mesh regeneration [1] to align the mesh with crack boundaries. Moreover, stress concentration at the crack-tip requires mesh refinement for accurate representation [3]. Mesh alignment and refinement operations negatively affect computational efficiency and accuracy of predictions [4]. The eXtended Finite Element Method (XFEM) initially was proposed by Belytschko and his collaborators [3, 5]. Belytschko and Black [5] introduced a new technique for solving crack growth problems with minimal remeshing. Their methodology was developed based on the work by Melenk and Babuška [6, 7]. The method provided an application of the partition of unity theorem to conventional FEM. XFEM works by enriching the nodes of the traditional finite element mesh by special shape functions to account for the displacement field discontinuities which are present in the case of a macrocrack. Hence, in XFEM domain, remeshing is no longer needed to account for the crack presence and its propagation. The method is well suited for crack propagation [8], and it was implemented in commercial Finite Element Analysis (FEA) codes such as ABAQUS and ANSYS [9, 10]. Meanwhile, the current implementation of XFEM in commercial finite element code has its challenges [11]. The primary problem is related to the need of embedding an initial crack a priori into the finite element mesh to trigger crack propagation. Inserting a crack raises the need for an expert user to identify the crack location. Alternatively, critical region(s) which are more likely to fail have to be determined for XFEM nodal enrichments. Hence without an expert user, enriching nodes of the entire domain of a finite element mesh becomes a general practice, which in turn results in more runtime and reduced computational efficiency. Heavily tied to the primary challenge comes a second one resulting from the dependency of prediction accuracy on mesh quality; more specifically the prediction of crack onset being dependent on mesh quality and density. XFEM predictions applied to composites has been on rapid development by researchers. For example, Grogan et al. [12] proposed a methodology for simulating thermal fatigue delamination in FRP composites using the framework of XFEM. Their methodology was used for the safe design of composite cryogenic fuel tanks. The proposed model predictions were validated with measurements from static and fatigue test methods. In their models, initial cracks were embedded to trigger crack propagation till failure. Duarte et al. [13] provided a comparative study between XFEM and Hashin’s damage criterion applied to the failure of composites. They concluded that the XFEM solution of the problem dealing with laminated composites overestimated failure loads at smaller deformations, suggesting that further investigation is required since XFEM predicted stiffer behavior. Petrov et al. [14] presented a parametric study for assessing the performance of XFEM applied to cross-ply composite laminates cracking. In their conclusions, XFEM predictions overestimated the number of cracks as well as corresponding strains and stresses. Bobiński and Tejchman [15, 16] used XFEM to study cracked concrete elements. In their work, they compared XFEM predictions to testing measurements of double-notched specimens reported by Nooru-Mohamed [14]. Load-deflection curves in [15, 16] reflected overpredictions of approximately 20% when compared to measurements. In review article on recent developments in damage modeling methodologies for composites [17], the multiscale modeling technique was highlighted as promising technique and promoted further investigation. Meanwhile Liu and Zheng [17] identified multiscale modeling computational challenges in incorporating microscopic modeling and macroscale failure mechanisms. As an example on multiscale modeling, Unger and Eckardt [18] presented multiscale modeling of concrete by coupling a homogeneous macroscale model with a heterogeneous mesoscale one. Their results were in a good agreement compared to testing results, but the multiscale modeling approach required long computational time (45000~64000 sec) for solving 2D problems. In the same review article Liu and Zheng [17] pointed out the advantages as well as the promising outcomes of XFEM which is the focus of the current research.
Therefore, in the current work, a novel approach is developed with the following objectives: (1) automatic identification of critical region(s) to alleviate user-dependency, (2) rigorous automatic mesh refinement based on stress convergence within these region(s) to ensure accurate predictions; and (3) automatically enriching critical region(s) for optimal XFEM execution. The main aim of current work is to eliminate the reliance on an expert user in identifying the critical region(s) and mesh refinement to enhance predictions accuracy at failure onset (damage initiation). For this purpose, an automation algorithm is developed to enable automatic identification of critical region(s) location/size and performing optimal mesh refinement procedure. The algorithm enriches only necessary nodes corresponding to the critical region(s) for XFEM modeling to predict the crack onset failure load and location. For the purpose of validating the current methodology, notched specimens were excluded since they are mainly used to study crack propagation rather than its onset. It is noteworthy to mention that the scope of the current work is predicting failure onset location together with failure loads/displacements with minimal user intervention to allow analyzing a real-life structure. To this end, a set of six unnotched concrete prismatic specimens were prepared and tested for validation. In addition further comparisons were established with relevant test data of unnotched specimens from the literature.
2. Research Significance
In the current work, some of XFEM implementation challenges are considered. The first challenge is the method dependency on user skills for critical region(s) identification for nodal enrichments. The second challenge stems from the high dependency of predictions accuracy on mesh quality and density. The proposed approach overcomes both challenges by automating the XFEM modeling process to arrive at a convergent mesh as well as potential (crack onset) zone without user intervention, hence allowing regular users to predict failure onset accurately. Another chief advantage of the proposed methodology is eliminating the need to embed initial cracks when analyzing crack propagation problems. This provides further advantages when the analysis is not limited to propagation and damage onset prediction is of primary importance. Critical load predictions at crack initiation facilitated by the current approach proved to be in excellent agreement with measurements obtained when testing unnotched specimens as well as relevant test data from the literature. Therefore, the proposed method allows accurate prediction of failure onset and eliminates the need for biasing specimens by introducing notches or first cracks. Furthermore, the proposed approach enables efficient mesh optimization and optimal enrichments which in turn enhances the overall computational efficiency. In conclusion, applying the proposed approach has significant effects on providing accurate and computationally efficient analysis of complex structures where critical region(s) identification can be challenging even for an expert user.
3. XFEM Fundamentals and ABAQUS Implementation
3.1. Mathematical Formulation
In the current section, the mathematical aspects of XFEM are briefly presented as in [3]. Consider a finite element mesh of a two-dimensional cracked body as shown in Figure 1(a), in which a crack is shown using a dashed line, finite element mesh is represented by solid lines, and nodes are represented by black filled circles (set I). As can be seen, in XFEM, the crack is surrounded by two types of nodes. The first type is the “Heaviside” enrichment nodes which are illustrated using square shapes (set J). The second type is the “Crack-tip” nodes which are presented using red circles (set K). The formulation of the XFEM problem is similar to that of conventional FEM formulation where both are based on standard Galerkin’s formulation with a slight difference in the former [9]. XFEM utilizes the global shape functions of conventional FEM throughout all nodes in the mesh (set I). A subset of the elemental nodes is referred to as “enriched nodes” (sets J and K) where extra terms are added to the global shape function to account for discontinuities in the displacement field (e.g., crack jump) in addition to the crack-tip singularities.
(a)
(b)
Therefore general shape function for all nodes in the domain takes the form ofwhere presents the global coordinates, are shape functions of , and are the degrees of freedom of node . is the Heaviside function or the jump function and are the shape functions related to the discontinuity at node , while are the additional degrees of freedom associated with the jump function. are the crack-tip enrichment functions, are the shape functions related to the crack-tip functions at node , and are the additional degrees of freedom related to the elastic asymptotic crack-tip enrichment functions. For 2D elasticity problem, the crack-tip enrichment functions are given bywhere represent the polar coordinate system with the origin at the crack-tip as shown in Figure 2. The tangent is at , and the outward normal at is denoted by .
Throughout the domain of the problem, nodes which are not enriched by a Heaviside function nor a crack-tip asymptotic function are associated with the conventional shape functions of FEM. Hence, (1) can be simplified to include only the first summation term on RHS leading to the traditional formulation of FEM reading as For the region which is censored by the crack (crack domain), the displacement approximation function of XFEM can be reduced to only include both first and second summation terms of (1) and can be written asFinally, to account for the crack-tip singularities as well as its propagation, (1) can be reduced to only include first and third summation terms on the RHS which takes the form of (5) as follows:As can be observed from the previous equations, the computational effort required for the solution of XFEM is higher than that needed for conventional FEM, because XFEM accounts for more degrees of freedom to capture crack behavior using the same number of nodes leading to an increased problem size, demanding more computational effort.
3.2. Enrichment Zone Sizing
In a 2D linear elastic boundary value problem shown in Figure 1(b), the crack domain is denoted by while the boundary is presented by . The boundary conditions are composed of three sets, namely, , , and such that . The displacements are imposed on , while tractions are imposed on . The crack surface is represented by which is assumed to be traction-free. The equilibrium equations and the constitutive relationships are given bywhere is the gradient operator, is Cauchy’s stress tensor, represents Hooke’s tensor, and is the strain tensor. The prescribed tractions arewhere is the outward unit normal vector to . Consequently for traction-free crack surface , . Under the assumptions of small strains and displacements, the kinematic equations readwhere is the linearized strains and is the displacement field. Equations (6) and (7) represent the strong form of the governing equations. In order to transform strong formulation of the problem into the weak form which is better suited for finite element computations [19], the displacement must belong to a set of kinematically admissible displacement fields [2, 5–7]. The weak form of the equilibrium equations is given by (9), which is solved using Galerkin’s method.Based on the weak form in (9), the error in the energy norm can be calculated [19, 20]. Gupta and Duarte [21] utilized the expression of the error in the energy norm to develop a priori estimate for the enrichment zone size. In their work, they explained that the error in the element located immediately outside a distance from the crack-tip is to remain bounded for optimal convergence. Under this condition and accounting for displacement field near the crack-tip, they developed an expression for optimal zone enrichment reading as follows:where is the minimal enrichment zone size for optimal convergence, is the polynomial degree of interpolation functions, and is the characteristic length of the element. To attain the optimal rate of convergence, the left-hand side of (10) is to remain less than or equal to constant . This equation was extensively studied in [21, 22] proving its validity.
3.3. XFEM in ABAQUS
Modeling process using XFEM for crack initiation in ABAQUS involves several steps: first, selection of regions that are more likely to fail or initiate a crack [9]. As an illustration, let us consider the Koyna dam problem shown in Figure 3: the Koyna dam survived an earthquake of magnitude 6.5 on the Richter scale on December 11, 1967 [23]. The quake did not cause significant damage, but it triggered the initiation of some cracks in the dam.
There exist two approaches to analyze this problem using XFEM. The first approach is to embed an initial crack in the critical region to trigger crack propagation analysis as shown in Figure 4(a). Merely this is done to study crack propagation rather than its onset. Embedding initial cracks reduces the accuracy of predicted failure limits, which is detrimental in the case of brittle materials as their failure is said to be catastrophic; once the damage is initiated, it will rapidly propagate under various types of loading till fracture. In the second approach, the user is required to identify potential failure region(s) which might be a straightforward task for a small structure (e.g., specimen). On the other hand, it might become a very challenging task dealing with large structures. Upon identifying the critical region(s), the user would select this region(s) as shown in Figure 4(b). Corresponding nodes for this region(s) are enriched with XFEM unique shape functions to account for crack onset. Failing to identify the critical region(s) will potentially lead to the enrichment of the entire domain of the problem. Full domain enrichment results in a drastic increase of required computational requirements and there is a possibility of an ill-conditioned system of equations that may cause convergence problems.
Initiation of a crack depends on a selected failure criterion for damage, typically stress or strain-based failure criterion. Accurate and efficient evaluation of the stress/strain fields for critical regions is dominant to precisely encounter crack initiation criterion. It is known that stresses in finite elements tend to converge at a slower rate than that of displacements which prompts the need of using refined and optimal mesh throughout the critical region(s) to precisely capture stress/strain fields [9]. Mesh refinement and optimization process is a user-dependent process relying on the user’s experience. Therefore, the lack of expert user insight entails enriching the entire domain of the problem and raises the potential of using less than optimal mesh, leading to a drastic increase in computational requirements along with lower predictions accuracy. Hence, overcoming the previously mentioned challenges is the primary objective of the current work. The following section devoted to present and discuss the proposed approach by emphasizing the role of each module.
4. The Proposed Approach
In the current work, a novel method is proposed to overcome XFEM modeling challenges in ABAQUS. The main aim is to automate the modeling process to predict failure onset (damage initiation) while maintaining optimal computational efficiency. It worth noting that the primary interest of current work is predicting failure onset without biasing or introducing any initial cracks a priori. The proposed algorithm was implemented using Python scripting in ABAQUS with three main tasks to be performed: first, automatic identification of critical region(s); second, mesh optimization for precise predictions; third, optimal XFEM execution for predicting crack onset location together with the corresponding load. A typical four-point bending problem was used as an example to demonstrate the algorithm procedures. A two-dimensional model was used throughout this section. The flowchart of the proposed algorithm is shown in Figure 5.
The main algorithm starts with reading geometric parameters, material properties, and loading conditions as user input parameters. The first primary step is an automatic determination of mesh size. For this purpose, a subroutine was implemented into the main algorithm to correlate the model geometry to the initial mesh size. A stress convergence condition was utilized for mesh refinement and optimization to ensure accurate predictions. Figure 6 illustrates the stress convergence performed by the optimal mesh subroutine for the problem in hand.
The second step is the essential one, where the algorithm is to check the model to identify the potential region(s) for crack onset. Predictions accuracy relies mainly on the proper selection of failure criterion which is strongly related to the material behavior. In the current work, brittle and quasi-brittle materials are the ones of interest. Hence, the failure criterion was chosen based on the third stress invariant . The work of Papadopoulos [24] concluded that fracture would occur when the stress tensor determinant reaches a critical value. The evaluation of a stress tensor will result in the same value independent of the material orientation, which can be considered as a material property [24]. The third stress invariant equation assuming zero out of plane stresses can be expressed aswhere and are the normal stress components, while is the shear stress component. The subroutine initiates iterations by incrementing loading magnitude while checking if failure criterion was encountered. During each iteration, the algorithm extracts data related to the most possible regions to fail. Figure 7 presents the third invariant stress field from the finite element simulation showing the critical region experiencing the maximum values.
Once the failure criterion is encountered, the algorithm automatically highlights critical zone(s) for crack initiation as depicted in Figure 8.
As can be seen from Figure 8, the algorithm highlighted a small zone of the beam’s lower mid-segment. By comparison with Figure 7, this zone is the one suffering maximum third invariant stress. The location is identified assuming perfect loading and boundary conditions. The critical region identification subroutine can also account for imperfect scenarios in load application/boundary conditions resulting from misalignments. In this case, the subroutine expands the critical zone to include its neighborhood. Critical zone(s) for crack onset are expanded to include regions of potential crack propagation direction. It worth noting that the crack propagation process is proven to be accurate in the standard XFEM method and consequently current implementation in ABAQUS. Meanwhile, without identifying a crack onset zone, the entire domain of the problem becomes a potential region for enrichment. It can be concluded that the region identification subroutine predicts the critical region(s) based on loading levels and imperfections. Figure 9 provides the final critical zone(s) as identified by the subroutine.
In the final step of the algorithm for optimal XFEM execution, a dedicated subroutine performs a mesh refinement to the automatically identified critical region(s). Figure 10 shows the refined mesh for the problem in hand. Subsequently, critical region’s nodes are enriched, and the algorithm solves for crack onset.
Figure 11 shows the enriched nodes of the model highlighted in red. The algorithm submits a new job to the ABAQUS solver, in which the critical region(s) have already been identified for nodal enrichments and the critical region mesh is optimized without user intervention.
5. Numerical Modeling
As discussed in Section 3.2, embedding a notch (crack) in the finite element model or the tested specimens may bias crack onset predictions. Therefore, to test and validate the proposed algorithm, a problem of an unnotched prism under four-point bending was selected. The problem geometry along with the loading conditions as per the American Standard for Testing and Materials (ASTM) designation for the selected problem is shown in Figure 12. For the detailed specifications of the problem, the reader is referred to ASTM C78 [25].
The total length of the beam is 400 mm centrally positioned on two supports spanned 300 mm apart. The beam is of square cross section with a side length of 100 mm. The load consists of two concentrated forces applied at the edges of the mid-segment of supported span allowing the beam to experience a uniform bending moment in this segment. It is noteworthy to mention that there is no need for a prior definition of the crack location as in the case of a notched beam. The crack location is to be identified automatically using the proposed algorithm.
The current test is conducted on a beam of brittle material, i.e., concrete. Brittle materials are said to undergo catastrophic failure due to their low strain to failure capacity. In other words, once a crack is initiated, it propagates under various levels of loading until fracture. The material model for concrete is chosen to be linearly elastic in compression and tension until failure. For failure initiation, a traction separation law based on maximum principal stress is adopted taking into consideration the ease of determining the maximum tensile strength for concrete by testing. The damage evolution was selected based on the fracture energy of concrete rather than the critical crack opening displacement. Fracture energy results in better accuracy due to the difficulty of capturing a crack behavior after being initiated for a brittle material [26], in addition to, the availability of detailed data on concrete fracture energy in the literature. A general static step was chosen for the analysis, and 4-noded bilinear plane strain quadrilateral element (CPE4R) with reduced integration is selected for meshing. The model is meshed using structured meshing control.
6. Specimens Preparation and Testing
General use Portland cement similar to ASTM (2012b) C150 Type I [27] was used to produce mixtures. Natural crushed stone was used as a coarse aggregate with a maximum size of 10 mm while natural sand was used as a fine aggregate. Both types of aggregates had a specific gravity of 2.6 and absorption of 1%. Several trial mixtures were cast to determine the water/cement ratio and coarse/fine aggregate ratio. Based on the trial mixtures stage, the optimized flowability of mixtures and enhanced mechanical properties were achieved by using a coarse to fine aggregate weight ratio of 0.7 and 0.9. Meanwhile, the water-to-cement ratios of 0.4, 0.45, and 0.5 were used. Table 1 provides the details of the mix designs used to build test specimens.
The material mechanical properties, namely, compressive strength, flexural strength, splitting tensile strength, and the modulus of elasticity for each mixture, were measured from standard testing. First, the compressive strength was measured according to ASTM C39 standard testing [28]. Second, the flexural strength was measured according to ASTM C78 [25]. Third, splitting tensile strength was measured according to ASTM C496 [29]. Finally, the modulus of elasticity was measured according to ASTM C469 [30]. The load cell accuracy of the testing machine was ±0.062% according to the last calibration report of the load-frame. The measured mechanical properties associated with each mixture are recorded in Table 2.
All mechanical properties were measured except for Poisson’s ratio which was assumed to be 0.2 [18] and fracture energy which is dependent on two main parameters, namely, size of aggregate used and compressive strength of the mixture. Given the knowledge of the maximum aggregate size and measuring the compressive strength, the fracture energy of concrete was calculated using (12) and (13) provided by Comité Euro-International du Béton (CEB) [31].where is the calculated fracture energy, is the compressive strength of mixture, and is a factor accounting for the maximum aggregate size used in the preparation of concrete mixtures.
7. Results and Discussion
This section is devoted to discussing the results from testing unnotched concrete prisms. The corresponding results predicted by the proposed algorithm are compared to those of testing for the validation purposes. Assessments are introduced delineating two significant aspects, namely, the critical load causing crack onset and corresponding crack location. Also, the proposed approach is compared to the conventional XFEM regarding computational efficiency.
The failure loads measured from testing of the six concrete prisms are reported in Table 3. Correspondingly, predicted values using the automation algorithm in addition to the percentage error are provided. As can be seen from the results, predictions compared to measurements indicated percentage error from 1.91% to 1.96%. The algorithm was able to precisely predicted failure loads with a relatively minimal error.
Table 4 provides a comparison between observed and predicted crack locations. The first column showing simulation predictions for crack onset compared to the failure from testing. Each figure was scaled down to one-quarter of its original length while preserving its aspect ratio. The second column is devoted to calculated percentage offset, which is based on beam’s mid-segment length. The distances were measured with a resolution of 1 mm. As can be observed, the minimum error value recorded was slightly more than 1% while the maximum was of 5%. The average offset distance value for the total number of specimens is less than 3%. It is noteworthy to mention that the predicted crack offset direction relative to the center of the sample is in the same direction as the offset compared to testing results.
Finally, comparisons were conducted regarding computational efficiency for proposed algorithm. Correlations assumed the absence of an expert user; hence conventional XFEM required enriching the entire domain. Since the proposed algorithms automatically arrives at convergent mesh, the same exact mesh was used for conventional XFEM assuming an expert user has arrived to this convergent mesh. Therefore, identical mesh size (i.e., the same number of elements) was utilized for the comparison purposes as shown in Table 5. The proposed algorithm enabled optimal enrichment of nodes in the critical region only. The total number of enriched nodes by the algorithm is 53% less than those of conventional XFEM. Consequently, the total number of unknowns in the system is 33% less than conventional XFEM requirements. Table 5 shows the detailed computational effort comparison. As can be seen from the comparison, the conventional implementation of XFEM results in almost doubling the total number of enriched nodes. Meanwhile, a minimal increase in the enrichment unknowns was enabled through the proposed approach resulting in runtime reduction by 39% and enhancing overall computational efficiency. Through this simple example, the reader can conclude the high potential of the current approach when applied to problems of complex geometry and combined loading scenarios. Furthermore, minimizing computational effort is significant in the case of dynamic analysis and 3D problems.
The following section is dedicated for comparisons with relevant test data (plain concrete testing) from the literature. The same aspects of comparisons outlined in Section 7 were used.
8. Algorithm Validation with Test Data from Literature
The proposed algorithm was tested for validation purposes against relevant test data from the literature. Only test data for unnotched specimens were compared with the algorithm predictions in the current section. Comparisons were based on damage initiation loads and/or displacements (if available). The first case adopted from literature for validation is from the work of Unger and Eckardt [18]. Their work presented adaptive multiscale modeling for concrete combining two different scale models, namely, mesoscale and macroscale. Nonlinear analysis of an L-shaped large-scale panel was provided using multiscale modeling. Their results were validated with specimens prepared and tested by Winkler et al. [32]. The mechanical properties as reported in [18] are reproduced in Table 6. The mechanical properties required for the current analysis are minimal compared to multiscale modeling techniques.
The problem geometry together with the loading conditions is shown in Figure 13(a). The proposed algorithm was used to analyze the L-shaped problem using the methodology illustrated in Section 4. The third stress invariant contour plot identifying the potential region for crack onset is shown in Figure 13(b).
(a)
(b)
Based on the maximum value of the third stress invariant, the algorithm automatically identified the potential region for XFEM enrichment. The critical region mesh is refined to ensure accurate failure predictions. Finally, the optimal enriched zone size is determined. The optimized final mesh along with the optimal enrichment zone constructed automatically by the proposed algorithm is shown in Figure 14(a). The crack onset location predicted using the proposed algorithm is shown in Figure 14(b).
(a)
(b)
The onset location is in excellent agreement with the test data reported in [18].
Table 7 shows the predicted failure load and displacement using the mesoscale-macroscale coupled model from [18] in addition to the predicted values using the proposed algorithm. The proposed algorithm recorded excellent agreement with measurements from testing in predicting both failure load and failure displacement. As can be observed from Table 7, the percentage error of the proposed algorithm shows enhanced predictions on both failure loads and displacements.
For completeness, the proposed algorithm was also compared to the conventional XFEM regarding computational efficiency as shown in Table 8. It can be observed that same meshes were used for the sake of comparison assuming that the user of conventional XFEM has performed mesh convergence analysis which was automatically performed using the current algorithm.
The algorithm enabled enriching less than 47% of the nodes in the entire domain for the current problem, in turn reducing the total number of unknowns in the system. As a consequence, the total number of required increments and iterations to solve using the proposed algorithm is significantly less than that required by conventional XFEM. It can be observed that the total runtime of the proposed algorithm is 42% in comparison to that of conventional XFEM. It can be concluded, from Tables 7 and 8, that the proposed algorithm enhanced the overall predictions while minimizing computational effort.
The second case for comparison is a large-scale T-section adopted from [33]. The specimen geometry and the loading conditions are presented in Figure 15(a). The section has a total height of 1 m and a hanging span of 0.6 m with a uniform thickness of 0.25 m. A concentrated tip load is applied to the free end of the hanging span. The proposed algorithm was utilized for analyzing the T-section problem following the same methodology discussed in Section 4. The contour plot of the third stress invariant is shown in Figure 15(b); based on its maximum values, the potential region for crack onset is automatically identified. The mechanical properties for the plain concrete mix as reported in [33] are reproduced in Table 9.
(a)
(b)
As per the current algorithm, the initial enrichment zone identification is followed by an automatic mesh refinement for the critical region. The final optimized mesh is shown in Figure 16(a); the automatically optimized enrichment zone is highlighted in red. The crack onset location as predicted by the proposed algorithm is shown in Figure 16(b).
(a)
(b)
The crack onset location is in excellent agreement with testing results from [33]. The failure load causing the first crack from testing was reported to be 30.9 KN, whereas the proposed algorithm predicted the first crack to occur at 30.6 KN. It can be concluded that the proposed algorithm is in excellent agreement with testing for both predicting the crack location and the failure load causing the first crack. Regarding computational efficiency, the proposed algorithm is compared to conventional XFEM as shown in Table 10.
The proposed algorithm was able to enrich approximately 50% of the entire domain resulting in a 50% reduction in the total number of unknowns. The required runtime to analyze this problem using the proposed algorithm was 36% of the runtime required by conventional XFEM.
9. Conclusions
In the presented work, a novel technique is proposed for automating XFEM implementation in user-built and commercial FE codes with minimal user intervention. The proposed technique was developed to overcome XFEM existing challenges related to prediction accuracy and to enhance computational efficiency. The novel technique was developed into an algorithm which was implemented using Python scripting in mainstream FE code ABAQUS performing three major tasks: (1) automatic identification of critical region(s) based on material-specific failure criterion, (2) automatic mesh refinement based on stress convergence for accurate predictions, and (3) optimal enrichment for identified critical region(s).
Three different problems were used for validation of prediction accuracy, namely, beam under four-point bending, L-shaped, and T-section problems. Furthermore, regarding computational efficiency, the performance of the implemented algorithm was compared to that of conventional XFEM. The following observations were recorded in these comparisons: (1) predicted failure loads/displacements corresponding to damage initiation were in excellent agreement with measurements of unnotched specimens obtained from literature and in-house testing; (2) crack onset location predictions showed excellent agreement with observations from all tested specimens of different geometries/loading conditions, and (3) significant reduction was observed in all model parameters in terms of total number of enriched nodes, number of unknowns, increments, and iterations when compared to conventional XFEM. It is noteworthy to mention that the runtime for three different problems under diverse loading conditions was in the order of seconds. Therefore and based on these observations it can be concluded that the proposed algorithm(i)alleviated the user-dependency for identifying the critical region(s);(ii)enabled accurate predictions due to rigorous and automatic mesh refinement;(iii)enhanced the overall computational efficiency.
In conclusion, applying the proposed approach was proven to have significant effects on providing accurate and computationally efficient analysis. Therefore this approach is recommended for analysis of complex structures and combined loading scenarios where critical region(s) identification can be challenging even for an expert user. Moreover, the approach possesses high potential of minimizing computational effort in the case of dynamic and 3D analyses.
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This work is sponsored by the President’s Doctoral Student Investment Fund (PDSIF) at Memorial University of Newfoundland and Natural Sciences and Engineering Research Council of Canada (NSERC) and Discovery Grant Program (NSERC DG no. 210415). These supports are gratefully acknowledged.