Mechanical design of a nested 4-layer Canted Cosine Theta (CCT) dipole

Nested CCT dipoles could be used in particle accelerators and proton therapy machines to bend a beam of charged particles in any direction. In the present study, a mechanical design of a 4 layer nested CCT dipole is evaluated with Finite Element Method (FEM). A full parametric 3D model of a 2.5 Tm 4-layer CCT dipole has been developed using the APDL scripting in the ANSYS software. The, so called, bottom up approach with direct generation of nodes and elements has been utilized, optimizing for the speed of the model generation as the time consuming meshing was bypassed. The properties of the Nb-Ti strands with the surrounding CTD-101K epoxy were obtained with a dedicated homogenization model. Resulting orthotropic properties were fully accounted for in the simulation. The shear stresses in the bonding composite layer: Kapton+S2-glass+CTD-101K were computed. With the results above the limit of 10 MPa, the necessity to provide additional rigidity against the torque was confirmed – castellated design. The shear stresses were mostly caused by thermal effects, and only ~25 % was caused by the Lorentz forces. The influence of the boundary conditions was analysed, leading to their optimal choice limiting the deformation due to Lorentz forces to 91 μm.


Introduction
The CCT magnets (originally named skewed solenoids [1]) are being strongly developed as offering many potential advantages over the cosine-theta, one of the most important being the lower stresses in the superconducting coils due to the support of each individual turn [2]. CCT magnets have flexibility of producing any multipole field, via control of the cable path. Up to now, 2 prototypes of CCT magnets wound with Rutherford Nb-Ti cables were build and tested [3] as well as three CCT dipole prototypes with Nb3Sn Rutherford cables [4]. The MCBRB (2.65 T bore field, 5 Tm integrated field) orbit correctors for the recombination dipole D2 for the High Luminosity upgrade of the Large Hadron Collider -HL-LHC [5] are made in CCT technology [6][7]. The development of the CCT dipoles for the HL-LHC brought the idea of nested CCT dipoles, with possible applications in particle accelerators (e.g. orbit correctors) and proton therapy machines [8]. Nested dipoles allow control of the magnetic field vector in the full range 0-360⁰. Currently, the PSI Institute is developing CCT technology with the goal of advancing the technology for 16 T magnets to satisfy the FCC requirements [9][10]. Recently a study showed validity of building CCT magnets from REBCO CORC ® HTS wires [11].
Several works have been devoted to 3D periodic FEM models of CCT magnets [12][13][14][15][16]. Whereas only one paper was devoted to the full 3D model of a CCT magnet [17]. In that work the experimental IOP Publishing doi: 10.1088/1742-6596/1559/1/012073 2 results from the strain gauges placed on the CCT4 prototype have been compared to the 3D simulations with good agreement.
In the present paper, a parametric 3D model of the nested CCT dipole is presented. The model was implemented in the APDL script in the ANSYS software to provide a full parametrization to the magnet's geometry. In addition, the meshing was bypassed via direct generation of nodes and element to maximally speed-up the pre-processing phase. This work is the direct development of the previous work [18] which was focused on developing a 3D periodic model and 2D models capable of modeling the nested CCT dipoles. Two types of boundary conditions were considered in this study according to the 1 st and the 2 nd design iterations. The interlayer composite was modeled according to the 1 st design, having thickness of 0.3 mm The results of the simulations allowed concluding on the final design choice of the castellated design for the engineering solution of the nested CCT dipole.

Geometry of a 3D CCT layer
The framework of describing geometry of a CCT layer was introduced in [18]. In the simplest case capable of modelling all the necessary components of the nested CCT dipole magnet, each of the layer contains 3 sub-layers: the coil, the spar (bottom part of the former) and the interlayer insulation. For such a geometry ( Figure 1) the parametric equations describing position of the corners (positions of nodes) in cylindrical coordinate system are [18]: where j={A, B, D, E, F, G, H, I, J, K, L, M, O, P}, i -layer number, -is the radius of the CCT axis (passing through the geometric center of the coil) of the i-th layer, is the angle in the xy plane of the center C, r -is the distance from the j point from the center C along the r-axis, b -is the distance of the point j from the center C along the b-axis, -is approximately the tilt angle of i-th layer, 0is the angle corresponding to the beginning of the i-th layer, -is the pitch of i-th layer, 0 -is the zcoordinate corresponding to the beginning of a turn in the i-th layer  After inserting the formulas for the unit vectors in cylindrical coordinate system ̂=̂cos( ) + sin( ), ̂= −̂sin( ) +̂cos( ), ̂=̂ into Equation 1 one obtains Cartesian equivalent: The geometrical parameters of the four layer nested CCT dipole magnet are listed in Table 1. In order to have a control over the mesh in the radial direction each of the 3 sub-layers ( Figure 1  In order to find the necessary values of the angle so all the points {A, B, D÷P} lie on the plane = the following equations needed to be solved:

Introducing Equation 2 and Equation 3 into Equation 7 leads to
Equation 8 has no close-form analytical solution. It should be noted that for the distance = 0, namely if the considered point is the center of the coil C the equation simplifies to = with the obvious solution θ= . The solution algorithm was implemented in the ANSYS APDL script, and it's described in the next section.

Finding intersections of the CCT edges with the plane θ=const
The plot of the Equation 8 with for the 1 st layer of the nested CCT dipole (Table 1) with the angle of intersection θk = θ01 is shown in Figure 3. The equation is always solved in the angle range θ0i..θ0i+360⁰, as the situation repeats for subsequent turns. The Equation 8 contains asymptotes for the denominator equal to zero: For θk=π/2+nπ the tangent tg(θk) in Equation 8 is infinite and the corresponding solution for the angle θ was found from Equation 9. For the 3 rd and 4 th layer as the angle θ0 is 270⁰ the equation 9 has 3 asymptotes, as shown in Figure 4. It should be noted that the solution of the Equation 8 and Equation 9 was always close to the angle θk. For the range of considered geometrical parameters (Table 1)     Depending on the number of angular divisions nθ, the equation 8 had to be solved 360•nP/nθ times per layer. The key aspect in the solution was setting correct angle limits that capture the solution around θk, generally ± 0.7⁰ (for the given geometrical parameters), and close to the asymptotes only on one side of them. The error tolerance was set to 2e-9, for lower values convergence could not be obtained close to the asymptotes. For this error tolerance the solution required less than 30 iterations.

The 3D FEM model of the nested CCT dipole
The model was meshed with n1=n2=n3=3 divisions, giving nP=14 (Equation 5). The angular step was 2.8125⁰ (giving exactly 128 elements around the circumference), so the Equation 8 needed to be solved 14*128=1792 times per layer. Base on the computed angles the nodes were generated from Equation 1 and then the elements. First 5 turns of the 1 st layer are shown in Figure 5a÷c. The end parts were generated base on the coordinates of the existing nodes (Figure 6a), 15 divisions were assumed for the end parts in the z-direction. Apart from 2nP (28), 6-node wedge elements (Figure 6b) per layer, all the other elements were hexagonal in shape. The side view on the 1 st layer without the interlayer insulation elements is shown in Figure 7. The full model with the outer cylinder is plotted in Figure 8a, the coils in Figure 8b, the interlayer insulation in Figure 8c, and the former channels in Figure 8d.
Solid 185 elements were used for the structure and CONTA174/TARGE170 elements for the contacts between the layers ( Figure 9). The option Initial penetration/gap was set to Exclude for all the contacts. There was in total 2323088 Solid185 elements with 2594000 nodes, plus 253200 contact elements. In total 3*2594000=7.782 MDOF.
The direct generation of the meshed model, show in Figure 8, took 19 minutes. With generation of keypoints and volumes and mesh subsequent meshing the generation time was of the order of 12 hours due to the very large number of volumes ~750 000. Writing a direct script made the meshing at least 20x faster allowing more rapid model development and faster parametric runs.

Material properties
Three material models were employed in the simulation. The formers and the outer cylinder were made out of high strength Aluminum 6082-T6 (Table 2) . The interlayer insulation was a composite made of Kapton and S2-glass impregnated with the CTD-101K epoxy. The interlayer composite was described with orthotropic material model with the thermal properties assuming 50 % of the S-2 glass fibers [22]. The mechanical and thermal properties of the coils (Figure 10a, b) were obtained via homogenization based on a separate 3D model (Figure 10c). Due to orthotropic properties of Nb-Ti strands plus the isotropic epoxy, the resulting properties were orthotropic as well. In order to account for the nonisotropic properties of the coils, element coordinate systems were rotated from the default Cartesian to the local frame rbt [18] (Figure 11).

Electromagnetic forces
The electromagnetic (EM) nodal force densities were mapped from the electromagnetic model (having denser mesh: 180 divisions circumferentially, 3 radially and 2 axially) to the element centroids in the mechanical model via ANSYS subroutine *MOPER,,MAP. The mapped force densities were multiplied by element volumes to obtain the forces. In order to visualize the distribution of the EM forces along one turn they were transferred to the local frame rbt (Figure 1 The distribution of the radial force is shown in Figure 12a. With maximum positive values near the angle θ=45⁰ for the layers 1 and 2 and maximum negative values for the angle θ=135⁰. The radial force was responsible for the oval shape of the magnet during powering (Figure 15). The binormal force (Figure 12b) was the largest out of r and t , whereas t the smallest (Figure 12c). The values of the forces r , b , t over the turns 100÷101 are plotted in Figure 13. The same values of forces were repeating along the entire straight section of the magnet. Forces in the ends were smaller due to lower magnetic field there. The period of variation of the force b was twice smaller compared to the forces r and t . Close to sinusoidal shape was observed for the forces r and t , whereas for the binormal force b considerable deviation from sine-like shape occurred.    (Figure 14b). The clearance of the yoke grooves is large enough, so the more contracting Aluminum cylinder can shrink independently from the yoke (smaller αT) in the axial and radial directions. In order to model such boundary condition the displacement of the node corresponding to the center of the keys was blocked in the circumferential direction along the whole length of the cylinder (Figure 14). The top surface of the cylinder was fully constrained in the z-direction, to represent the test in a vertical cryostat.
The castellations in the formers, visible in Figure 14 were not included in the geometry in this study. Here, we wanted to check if the structure can satisfy the design constraints only relying on bonding between the layers. The key issue was checking if the maximum shear stresses in the interlayer insulation were below the assumed design limit of 10 MPa. The answers are given in the next section.

Results of the simulations for the nested CCT dipole
Direct solver was used for the solution with the Distributed Memory Parallel mode. The model required ~550 GB of RAM memory. As the system was limited to 192 GB, the computation run in the out of core mode on fast SSD M.2 PCIe discs. The analysis consisted of 2 steps: cool-down from 293 K to 1.9 K and loading with EM forces. The results related to EM forces only were obtained in post-processing by subtracting the results of the 1 st load step from the results of the 2 nd load step. The 1 st step required 5 iterations of the Direct Solver (factorizations) and the 2 nd step six, with the average amount of written data to SSD discs of ~3.5 TB per iteration. Giving ~38.5 TB for the entire simulation, so an important number for considerations of SSD wear. Solution time was ~24 hours (with 16x3.5 GHz cores).
The maximum total deformation due to EM forces only for the diagonal BCs (as in Figure 14b) was 161 µm (Figure 15), so considerably more than for the case with symmetric BCs (as in Figure 14a) with 92.3 µm (Figure 16). For the field quality reasons, the deformation of the coils should be ideally below 50 µm. The symmetric BCs are therefore beneficial in lowering the deformations due to EM forces. The maximum deformation due to EM forces only for the symmetric BCs did not lie in the center of the magnet, but rather was shifter towards the free end (Fig. 16), an effect not observed for the diagonal BCs.
The shear stresses in the interlayer composite for the model with diagonal BCs are shown in Figure  17 and Błąd! Nie można odnaleźć źródła odwołania.. The maximum value after cool-down was 32.7 MPa and 35.4 MPa after powering, whereas due to EM forces only 11.8 MPa. These results clearly showed that the thermal stresses were 3x larger than Lorentz force induced stresses, so actually the layers would debond already due to the cool-down to 1.9 K. For the symmetric BCs, maximum shear stresses due to cool-down were slightly smaller: after cool-down 32.7 MPa (Figure 19), after powering 34.2 MPa and due to powering 8.0 MPa (Figure 20).   MPa. All shear stresses caused by powering only were below 7.0 MPa, for both boundary condition types, indicating that the debonding from the Aluminum formers is not likely to occur during powering -assuming that the shear limit of 10 MPa holds. The occurrence of debonding is of a large importance as it could become a source of a, very undesirable, training in the magnet. The compressive stresses in the interlayer insulation were low w.r.t. the compressive limit of 1360 MPa at 4 K [22], with maximum value in the radial direction for the symmetric BCs during powering (2 nd step) of 48 MPa, in the circumferential direction of 156 MPa (local effect) and 162 MPa in the axial direction.

Conclusions
The analysis of the 3D model confirmed the necessity of including castellations in the design, as shear stresses were much above the shear strength limit of 10 MPa. The maximum shear stress , after cooldown to 1.9 K, was 32.7 MPa for both diagonal and symmetric boundary conditions, whereas the maximum shear stress due to electromagnetic forces only was 11.8 MPa and 8.0 MPa respectively. So the maximum shear stresses caused by temperature were 2.8 and 4.1 times larger than stresses caused by the Lorentz forces. Large thermal stresses are caused by the incompatibility of the coefficients of thermal contraction between the Aluminum (formers) -α293-1.9=14.2e-6 and the fiberglass impregnated with CTD-101K resin -αx_293-1.9=20.9e-6 (through thickness). In the simulation the volume fraction of the resin in the interlayer composite was assumed as 50 %, if the content of S-2 glass fibers would be larger the thermal contraction coefficient would decrease. However, increasing the fiberglass content can be difficult, as it can hinder the resin flow during impregnation. Therefore, the assumption of 50 % content of fiberglass looks realistic. Another possibility of decreasing the thermal stresses in the interlayer would result from employing other type of epoxy resin with lower coefficient of thermal contraction. However, this would require extensive R&D work in verifying the mechanical behavior of such a resin and its performance with complex coils -such as for the nested CCT dipole.
The symmetric boundary conditions show considerably smaller deformations of the coils -92.3 µm versus 161 µm, as well as smaller shear stresses, therefore they should be considered in the final design due to field quality reasons. The FEM modeling of the castellated nested CCT dipole will be the subject of the subsequent work.