Coupled FEM-DBEM method to assess crack growth in magnet system of Wendelstein 7-X

A BSTRACT . The fivefold symmetric modular stellarator Wendelstein 7-X (W7-X) is currently under construction in Greifswald, Germany. The superconducting coils of the magnet system are bolted onto a central support ring and interconnected with five so-called lateral support elements (LSEs) per half module. After welding of the LSE hollow boxes to the coil cases, cracks were found in the vicinity of the welds that could potentially limit the allowed number N of electromagnetic (EM) load cycles of the machine. In response to the appearance of first cracks during assembly, the Stress Intensity Factors (SIFs) were calculated and corresponding crack growth rates of theoretical semi-circular cracks of measured sizes in potentially critical position and orientation were predicted using Paris’ law, whose parameters were calibrated in fatigue tests at cryogenic temperature. In this paper the Dual Boundary Element Method (DBEM) is applied in a coupled FEM-DBEM approach to analyze the propagation of multiple cracks with different shapes. For this purpose, the crack path is assessed with the Minimum Strain Energy density criterion and SIFs are calculated by the J-integral approach. The Finite Element Method (FEM) is adopted to model, using the commercial codes Ansys or Abaqus;, the overall component whereas the submodel analysis, in the volume surrounding the cracked area, is performed by FEM (“FEM-FEM approach”) or alternatively by DBEM (“FEM-DBEM


INTRODUCTION
he fivefold symmetric modular stellarator Wendelstein 7-X (W7-X) is currently under construction in Greifswald, Germany [1]. Each of the five modules are build up of two flip symmetric half modules (ten in total) which consist of 5 non-planar (NP) and 2 planar superconducting coils each to be operated at 4K (Fig. 1). The NP coil cases are made of cast stainless steel (SS) EN 1.3960. They are bolted onto a SS central support ring and welded together at the outboard side of the torus via the LSE's consisting of 100 -150 mm long hollow tubes of 30 -35 mm thick forged SS EN 1.4429 (Fig. 1). After welding the LSEs (weld depth 15-30 mm), many surface cracks substantially larger than 8 mm (typical acceptance limit of EN 23277 [2]) were found with dye penetration tests at the accessible surfaces, particularly at the coil side of the weld within the cast steel, oriented parallel to the weld seam. In a first paper [3], the averaged SIF along the crack front of these cracks were predicted using an analytical method and finite element method (FEM) models which were benchmarked against each other and against a single run with a Dual Boundary Element Method (DBEM) model. Assuming proportionality between the SIF and square root of the crack size ( K a   ) and an arbitrary maximum crack size, it was possible to derive the number of load cycles until failure from the SIF of the initial crack. The maximum radial crack size increase was limited to 10 mm to avoid coalescence of the T observed crack with adjacent cracks and to ensure that the crack remains small compared to the cross section, i.e. to ensure the proportionality rule K a   . Based on Paris' law parameters derived from 10 fatigue crack growth rate tests (FCGR) carried out at cryogenic temperature at Karlsruhe Institute of Technology (KIT), the predicted SIF could be related to cyclic crack growth [3]. In the current paper, the assessment is extended, for the most critical crack, using the DBEM [4][5][6][7][8][9][10][11]. This most critical crack was located in LSE-D05 in the cast steel (Fig. 1). In the DBEM method, the stress state and SIFs along the crack front are updated at each step of crack advance; as a consequence, deviation from the proportionality rule is automatically simulated. Moreover, an adjacent crack was included in the model which allows for the simulation of crack coalescence. So, no maximum crack size needs to be defined and crack growth can be continued until the critical SIF is reached, i.e. the SIF where unstable crack growth starts. In addition, the method calculates the crack growth rate along the crack front rather than assuming a pure radial growth as done in [3] with a pure FEM approach. So, it enables the prediction of crack growth of non-circular cracks, differently from what imposed by the FEM sub-modelling approach, where only semicircular crack fronts are modelled [3]. The depth of the crack is fundamentally unknown but from repair experience it was found that the crack depth was typically smaller than half the crack length. To verify the influence of the assumed crack shape, a semi-elliptical crack has been compared with the original semi-circular crack. In the next section, the DBEM and FEM models are presented and the results of the FEM and DBEM submodels are compared. In the DBEM submodel also the effect of an elliptical crack, adjacent to the main crack is evaluated. Finally the conclusions are given.

MODELS AND RESULTS
ince the mechanical behaviour of the magnet system is of crucial importance for the operation of W7-X, two independent global FEM models of the magnet system have been developed, respectively using the commercial codes Ansys and Abaqus, which are successfully benchmarked against each other. The displacements and generalised sectional forces and moments typically deviate less than 10% between both global models when loaded with a magnetic field on the plasma axis of 3 T (Tesla). An FEM submodel analysis around the most critical crack has been made in Abaqus: a crack is introduced in such submodel (loaded on the cutting boundaries with displacements from the Abaqus global model) and the SIF's are calculated along the front. Two different DBEM submodels are considered: one extracted from Ansys (loaded on the cutting boundaries with displacements from the Ansys global model) and another one extracted from the Abaqus submodel, with a further closeup around the cracked area (again the boundary conditions on the cutting face are taken for the Abaqus submodel).

FEM model and submodel
Preparing the FEM submodel, it was found that it suffices to include in the submodel only the LSE without the adjacent coils. So only LSE D05 was modelled with the semi-circular crack of 14 mm diameter at the outside of the cast steel (Fig.  2). The crack is modelled as a seam in the mesh with elements following the contour of the crack. At the same time, it was found that the results are sensitive to the weld modelling. Notably, the weld does not penetrate through the entire thickness of the forged steel tube, so part of the cross section remains un-welded (see right picture of Fig. 2). If the unwelded area is at the outside of the cross section, it shields the crack from stresses perpendicular to the crack. Vice versa, the perpendicular stresses around the crack are increased when the un-welded area is inside of the cross section. For the critical crack under consideration, the stresses perpendicular to the crack are increased. Since the mesh is adapted to the crack shape and size, the model does not allow an automatic crack growth. In Fig. 3 the Von Mises stresses for the cast steel around the crack, under an EM load equal to 3T, are presented for both the Ansys and Abaqus models (and submodel): it is to be pointed out that the LSE geometry is less accurate in the Ansys global model than in the Abaqus global model so that a rigorous comparison between the stress scenarios is prevented.

DBEM submodel
The DBEM sub-modelling is aimed at analyzing the propagation of cracks with elliptical or semicircular shapes and, eventually, affected by nearby cracks. For this purpose, the Dual Boundary Element Method (DBEM) is applied in a coupled FEM-DBEM approach [12][13][14][15] with the crack path assessed by the Minimum Strain Energy density criterion and the Stress Intensity Factors (SIFs) calculated by the J-integral approach [4][5][6]. The Finite Element Method (FEM) is S adopted to model the overall component, whereas the submodel analysis in the volume surrounding the cracked area is performed by DBEM. The boundary conditions, applied on the DBEM local model, are the displacements calculated from the ANSYS global model (Fig. 4), with related material properties listed in Tab. 1, or from the Abaqus submodel ( Fig. 5) with related material properties listed in Tab. 2.       A DBEM thermal-stress analysis [9] is performed on the submodel extracted from Ansys (Fig. 4) in order to allow for the coil cooling from an ambient temperature equal to 293 K to the superconductive temperature equal to 4 K, and subsequent electromagnetic loading from 0 to 3 T. In Fig. 6 the DBEM submodels, with maximum principal stresses (MPa) are shown, in the configuration before and after the crack insertion (the initial crack is semicircular with a radius equal to 7 mm). A DBEM stress analysis (this time the thermal part is considered with negligible effects, namely the stress coming from cooling at 4 K are neglected) is performed on the submodel extracted from Abaqus ( Fig. 5) with an electromagnetic loading equal to 3 T. In Fig. 7 the DBEM submodels, with maximum principal stresses (MPa) are shown for comparison with the results illustrated in Fig. 6: some limited differences pop up as expected due to the huge difference in submodel and global model mesh refinement, for not mentioning the geometry differences; there is also a small contribution to the difference coming from the different loading process (in one case there is no allowance for cooling effects).

FEM-FEM vs. FEM-DBEM approaches
n Fig. 8 it is possible to compare the J-integral values calculated along the initial crack front with the 3 different approaches available: the FEM-DBEM approach provides two sets of results depending upon the software (Ansys or Abaqus) used for the global modelling, whereas the remaining set of results is related to the crack assessment done on the Abaqus submodel. As expected there is a non negligible discrepancy when considering the results coming from the Ansys global model but it is easy to find the reasons in the rough mesh and in the lack of a detailed geometric reproduction of the LSE; on the contrary a satisfactory correspondence is obtained between the results obtained by FEM and DBEM submodels when considering the Abaqus global model, apart from a non negligible discrepancy localised at the break through points that are more affected by extrapolation effects and mesh inaccuracy. With reference to the latter, with reference to DBEM analysis, when the J-path is built on a triangular element there is a loss of accuracy in the J calculation, but this can be overcome by further refining the mesh on the crack so as to minimise the impact of the last unavoidable triangular elements at break through point (see in Fig. 7 the close up of the mesh along the crack front).

DBEM crack propagation analysis
The simulation of crack propagation, in this work, is only performed with reference to the DBEM submodel extracted from the Ansys global model. The crack propagation law adopted is the Paris law (Eq.1) calibrated under cryogenic conditions (T=4K) [3]: with C=7.95·10 -10 (with da/dN measured in mm/cycle and K in MPa√m) and m=3.23; the plane strain fracture toughness is Ic K =106 MPa√m. The fatigue cycle considered is based on the coil cooling, from ambient temperature to 4 I K, and subsequent application of the electromagnetic load corresponding to a magnetic field cyclically varying from 0 to 3 T. After 6 steps of crack propagation, the crack front break through the LSE superior surface as shown in Fig. 9, where again the maximum principal stresses (MPa) and the final crack shape are shown.  The mode I SIFs (model II and mode III SIFs are much lower) along the crack front for each step of crack advance are shown in Fig. 10: it is possible to see the sudden increase of KI, in the final step, along that part of the crack front that breaks through the upper LSE surface (see Fig. 9). One of the aims of this work was to check the influence on the fatigue life of a different initial crack front geometry, so it was also considered an alternative case with an initial semi-elliptical crack with the same value of a=7 mm and with a depth reduced from c=7 to c=5 mm. The fatigue cycles for each crack advance for the two aforementioned initial configurations are shown in Fig. 11: it is clear the impact of the initial crack geometry. The simulation is ended when the crack breaks through the upper LSE surface, because from now on the SIFs will rapidly increase approaching the fracture toughness.  The second aim of this work was to analyse the reciprocal influence of adjacent cracks and consequently a multiple crack propagation analysis has been performed. The two initial cracks considered are respectively semicircular, with radius equal to 7 mm (the N.1 in Fig. 12), and semielliptical (N. 2 in Fig. 12) with semi-axis a and c (depth size) respectively equal to 7 mm and 5 mm; the distance between the centres of the two aforementioned cracks is equal to 28 mm; the corresponding stress scenario is illustrated in Fig. 12.
The stress scenario and the crack configuration after four steps of crack propagation are shown in Fig. 13: it is possible to see that from this stage on the crack coalescence can be considered likely to happen.  Differently for the previous case, this time the SIFs are calculated by the crack opening displacement method (COD), [4] because there were difficulties in enforcing on the secondary crack a mesh along the crack front made of mainly quadrilateral elements and the J-integral accuracy, as previously said, is affected if the J-path is built on a triangular element. The calculated SIFs are shown for the two cracks at each crack advance in Figs. 14-16. It is to be pointed out the formula used (Eq. 2) for calculating the equivalent SIF used in the Paris formula is: (2)   The fatigue cycles for each crack propagation step for the two aforementioned cracks are shown in Fig. 17: comparing the double crack results with the single crack results (Fig. 11) it is possible to observe a clear detrimental effect of the second crack on the main semicircular crack (e.g. in case of multiple cracks, the main crack depth reaches a size of about 2 cm after 4.5·10 4 cycles rather than after 7·10 4 cycles as in the single crack case). Again the simulation is ended when the main crack breaks through the upper LSE surface, because from now on the SIFs will rapidly increase towards the fracture toughness.

CONCLUSIONS
n this paper some partial results of a benchmarking activity between two different approaches for the crack growth assessment in a component of the magnet system of Wendelstein 7-X are presented: both approaches are based on the sub-modelling of the critical cracked domain, as extracted from an FEM overall model, but in one case the submodel is still FEM based whereas in the second case the submodel is DBEM based. The latter approach realises a synergic coupled approach between the two numerical procedures, providing a more flexible tool for analysing complex crack scenarios. The adopted codes are ANSYS and ABAQUS for the FEM modelling and BEASY for the DBEM modelling. The SIFs calculated with the FEM-FEM approach and with the FEM-DBEM approach can be successfully compared on the initial crack front. The effect of the initial crack shape and the presence of an adjacent crack have been evaluated: assuming a non-circular crack of 5 mm depth and 14 mm width, instead of the initial circular crack with 7 mm depth, leads to a reduction of crack growth and increases the number of cycles up to reaching K Ic from 7·10⁴ to 8.5·10⁴. On the other side, an adjacent semicircular crack of equal size of 7 mm radius at initially 28 mm centre to centre spacing accelerates crack growth and reduces the number of cycles to 4.5·10⁴. Therefore, when assessing individual cracks neglecting eventual adjacent cracks, it is necessary to prevent coalescence of the cracks by limiting crack growth to half the crack spacing on top of limiting Keff to K Ic .