A Modified Mechanical Model for the Optimization of High Field Solenoid Magnets

A modified mechanical model for the stress analysis of solenoid pulsed magnets is presented, which includes the influence of bending on the plane stress analysis. The bending effect in the axial direction is included by introducing the discrete axial strain with the assumption of the uniform distribution of the axial force load. The modified model improves the accuracy remarkably compared to the classical Generalized Plane Strain (GPS) analysis. The results are verified with the FEA software. The algorithm of identifying the optimum positions of the separations between windings is also proposed. This model achieves both high efficiency and high precision, which are of great significance to optimize high field solenoid pulsed magnets and explore the existing potential for higher magnetic fields.


I. INTRODUCTION
High magnetic field is a basic research tool that is extensively employed in various scientific areas [1]. It has contributed to many important scientific discoveries, such as the Nobel Prizes of the integer and fraction hall effects [2]. Pulsed magnets can provide higher fields (up to 100 T) than DC magnets (45 T), and they are preferred by some frontier scientific researches [3]- [10]. However, the field strength of pulsed magnets is severely constrained by the massive Lorentz force exerted on the conductors [11]. The magnetic stress generated in a 100 T pulsed magnet is 4 GPa, which exceeds the ultimate tensile stress (UTS) of any metal or metal alloy that are commercially available and approaches closely to the UTS of high strength fiber materials such as Carbon or Zylon [12], [13]. Any minor design error will result in the malfunction of the magnet [14], [15]. Therefore, an accurate mechanical analysis is very important for the optimization of pulsed magnets to tap the existing potential for higher fields.
In a pulsed magnet, the Lorentz force in the radial direction expands the magnet outwards, resulting in hoop and radial stress. Meanwhile, the axial Lorentz force compresses the magnet axially. Due to the massive force, the The associate editor coordinating the review of this manuscript and approving it for publication was Jenny Mahoney. conductors (copper or copper alloy) are always stretched into the plastic region [16]- [18]. In order to sustain the Lorentz force, both internal and external reinforcement made from orthotropic fibers of Zylon or Glass with high modulus and high strength are introduced to support the conductor layers [19], [20]. Besides, the free-standing technique is usually applied between the inner layers to eliminate the stress concentration [21], [22].
Under high magnetic field, the magnet tends to bend axially due to the uneven distributed Lorentz force. The non-linear plastic property of the conductors, the orthotropy of reinforcing fiber, the separations between the inner layers and the bending caused by the shear stress altogether complicate the calculation of the stress in a high field pulsed magnet. Commercially available finite element software can solve the mechanical model, but they are time-consuming and likely to fail to converge because of the nonlinearity and contact analysis in pulsed magnets [23]. Meanwhile, the design and optimization of pulsed magnets is an elaborate process that requires a substantial repeated calculation to acquire an optimum structure. Therefore, an accurate and efficient model for the mechanical analysis of pulsed magnets is preferred.
A lot of work has already been done on this subject [22], [24]- [29]. A simplified model of a solenoid with infinite length was firstly presented [24], but the axial stress is not included and the accuracy is not enough. The Generalized Plane Strain (GPS) model was proposed by neglecting the shear stress and focusing on the stress distribution in the mid-pane of the solenoid magnets that has the maximum stress [25], [26]. Owing to the high efficiency and accuracy, the GPS model is widely adopted in the stress analysis of superconducting and pulsed magnets. Analytic solution of the GPS model was extended to include plastic deformation of the conductors and the shear stress in more complicated two-dimensional coordinates by a power series expansion [27], [28]. In order to solve the non-linear property of the conductors in the plastic region, numerical method is presented based on the GPS model, and a proprietary and professional mechanical analysis software PMDS was developed [22], which calculates the stress in the mid-plane in a finite difference method and shows great efficiency in pulsed magnet design.
However, in the classical GPS model, the magnet is considered as a right cylinder with no bending. This assumption is valid in thin layers where the effect of the bending on the mid-plane is negligible. When many layers are glued together to form a thick ''layer'', the radial uniformity of the axial Lorentz force deteriorates as the thickness of the ''layer'' increases and the degree of the bending becomes more significant. This is especially true for the pulsed magnets with external reinforcement, where the conductors in the plastic region undergo large deformation and become more susceptible to bend under the uneven force load. The bending effect can alter the stress distribution on the mid-plane dramatically, and results in different separation conditions.
In this paper, the Modified Plane Strain (MPS) model for the stress analysis of high field pulsed magnets is proposed. The influence of the bending is included by introducing the approximate discrete axial plane strain. The separations between the layers are analyzed with the MPS model considering the bending effect. An algorithm of identifying the optimum locations of the separations between the layers is presented for the optimization of the stress distribution. The proposed model shows a great improvement in accuracy without sacrificing computational efficiency.

II. MATHEMATICAL FOUNDATION
The picture and the schematic diagram of the cutaway view of a typical pulsed magnet with both internal and external reinforcement are given in Fig. 1 [19]. The magnet consists of N layers of conductor. Each layer of conductor is supported by the adjacent layer of reinforcement, and together they are regarded as a couple (CP). The fourth to the N th CPs are glued and form a Mechanical Independent Unit (MIU), which has free boundary conditions for stress calculation in the radial direction. The inner three CPs separate due to the magnetic force and become MIUs respectively [21], [22].

A. MAGNETIC FIELD AND LORENTZ FORCE
In cylindrical coordinates, the magnetic field generated by a solenoid layer with uniform hoop current density J θ can be calculated by equations presented in [30]. For a solenoid layer with inner radius a 1 , outer radius a 2 , height H , and the origin of cylindrical coordinates set at the center of the solenoid, the axial and radial components of the magnetic field B z and B r at any position (r, z) can be obtained as The parameters in (2) can be found in Appendix A. The total magnetic field can be obtained by superimposing the field generated by each conductor layer. The interaction of the magnetic field and the current density gives rise to a distributed Lorentz force density as given in (3), where f r and f z are the Lorentz force per unit volume in r and z direction respectively.

B. STRESS BALANCE EQUATION
The cylindrical form of the general stress equilibrium equations on the r-z plane for an axisymmetric solenoid magnet with a distributed force is well known as [26] ∂σ r ∂r + and the constitutive relations can be written as  where σ is the normal stress, τ the shear stress, D and C the components of the stiffness matrix. In the GPS model, the constant axial strain is introduced on the assumption that the magnet remains upright during the deformation with no bending and no shear stress. By neglecting the shear stress on the mid-plane, the stress equation can be reduced to Considering that the solenoid magnets are symmetric to the mid-plane, the axial component of Lorentz force f z is zero in the mid-plane. This leads to a constant ε z in the z-direction The stress balance equation can be formulated with the radial displacement u by introducing the strain-displacement relationship

C. AXIAL FORCE BALANCE
The axial force balance on the mid-plane requires b a σ z 2πrdr = F ztotal (11) where a, b are the inner and outer radius of a mechanical independent unit (MIU), F ztotal the total axial body load on the mid-plane of the MIU.
In the classical GPS model, the axial strain is assumed to be constant throughout an MIU. This assumption is valid in MIUs with small thickness, which remains almost upright during deformation as the inner three MIUs in Fig. 2(a). However, in the thick glued MIU, such as the 4∼10th CPs in Fig. 2(a), severer bending occurs in the axial direction, resulting in a radially varying axial strain on the mid-plane.
The assumption of constant axial stain is no longer reasonable.
The axial bending of the thick MIU is due to the uneven axial Lorentz force, which is generated by the radial magnetic field B r . Since B r peaks at the axial edge and declines dramatically towards the mid-plane as shown in Fig. 2(b), the axial compressive force mainly concentrates at both ends of the conductors. The distributed axial Lorentz force gives rise to the shear stress between the layers, which reaches its maximum near the ends of the windings and declines rapidly toward the mid-plane as shown in Fig. 2(a). It tends to redistribute the axial force over the whole MIU more evenly on the mid-plane, resulting in a macroscopical uniform distribution of axial stress among each CP.
In order to include the effect of the bending on the stress distribution, the MPS model is proposed by assuming that the equivalent average axial stress σ z on the mid-plane is uniform in all CPs except for the CP with the external reinforcement layer and calculating the corresponding axial strain in each CP. The distribution of the σ z in an MIU consisting of N CPs is given in Fig. 3. Since the outermost reinforcement layer of most pulsed magnets is designed to be thick with a  large mechanical strength margin for stress, the axial strain at the outer diameter is negligibly small and can be taken as zero ( r out = 0). The axial force load on the conductor couples through the shear stress to the external reinforcement. The coupled axial load peaks at the inner radius and declines to 0 at the outer radius of the external reinforcement. Here, a linear decline of the σ z in the outermost reinforcement layer is introduced for simple approximation.
By integrating the axial force f z over the whole volume of the MIU, and then averaging the total axial force load on the mid-plane according to the stress distribution in Fig. 3, the equivalent average axial stress σ z (r) can be obtained as In the MPS model, the thick MIU is divided into several CPs and the axial strain is considered to be constant in each CP, neglecting the shear stress on the mid-plane. Each CP can be seen as a small unit with the composite properties of conductor and reinforcement. Notice that the axial strain of the MIU is discontinuous in r direction. It is an equivalent value in each CP corresponding to the redistributed axial force load, rather than the true axial strain. The redistributed axial force load F zi in each CP can be calculated by integrating the σ z over the area S i of each CP on the mid-plane With the redistributed axial force load, the axial force balance equation in each CP can be obtained respectively as (14) and the axial strain of each CP ε zi can be calculated by solving the axial force balance equation (14) in each CP together with the stress equilibrium in (7). In the MPS model, the axial strain in a thick MIU is a radially discrete variable, which is an approximation of the continuously varying axial strain caused by the bending.

D. SEPARATION OF LAYERS
The free-standing technique is employed at the inner layers of the windings to improve the stress distribution. To take the separation into consideration in the design and optimization of pulsed magnets, it is crucial to identify the optimum locations of the separations between the CPs.
For a compact pulsed magnet, each layer of the windings is tightly wound and they can be treated as sticking together. It is determined whether the CP tends to separate according to the radial stress. Tensile radial stress means pulling between CPs and the free boundary is applied to eliminate the stress concentration. In practice, thin layers of Teflon are inserted between these CPs to facilitate the separations [22]. However, the introduction of the free boundaries will redistribute the stress and strain, and the positions of possible separations will also change accordingly. The initial radial stress of all the layers stuck together cannot reflect the optimum locations of the separations clearly, which might result in an excessive preset number of the separated CPs and an incorrect stress analysis.
A reliable method is to determine the separations of the CPs by the radial displacement. When two MIUs are set to be independent, the criterion for true separation is that the radial displacement of two MIUs at the interface follows where U i−1 r is the radial displacement of the inner MIU at the outer radius, and U i r the radial displacement of the outer MIU at the inner radius. As depicted in Fig. 4, the MIU i−1 and the MIU i satisfy (15) and separate from each other. Whereas, the MIU i and the MIU i+1 overlap in radial displacement, which means that extrusion occurs between them.
Based on this criterion, an algorithm for automatic identification of the optimum locations of the separations is proposed. Initially, assume all the CPs to be MIUs. Then, combine the MIUs in contact into one MIU according to the radial displacement and repeat this process until all the MIUs satisfy the criterion of separation. Eventually, the optimum locations of the separations are found.
The criterion of separation relies on the calculation of the radial displacement. However, the bending has an impact on the axial strain, which can couple through the Poisson's ratio and changes the radial displacement, resulting in incorrect separation conditions. So a closer approximation of the axial strain in the MPS model will contribute to identifying the precise locations of the optimum separations.

III. FINITE ELEMENT MODELING
In order to solve the stress in the mid-plane of a pulsed magnet with discontinuous media structure and nonlinear material properties, the finite element numerical method is chosen for its high accuracy, stability and reliability.
Substituting the constitutive relation (6) and strain displacement relation (10) into the stress equilibrium equations (7) and (14), the stress equations for each MIU can be obtained as 92878 VOLUME 8, 2020 where D ij are the elements in the constitutive matrix, B z the z component of the magnetic field, J θ the current density of the conductor layer, and all of them are varying with radius. m is the number of CPs in the MIU, ε zi the uniform axial strain of the ith CP, and F zi the correspondent axial force load in the ith CP. The Galerkin-weighted residual method [31] is introduced here to get the finite element discrete form of the stress equations. The one-dimension meshing of an MIU on the mid-plane of a pulsed magnet is shown in Fig. 5. The MIU is divided into n-1 elements and n nodes. U i is the displacement of each node to be solved, x i 1,2 the interpolation of the ith element.
Usually, a pulsed magnet has several MIUs, which can be solved one by one. Each MIU may compose of m CPs with different axial strains that correspond to the redistributed axial force loads. Accordingly, a total stiffness matrix can be assembled up from one element to another in each MIU (see Appendix B). By substitution, (16) and (17) can be integrated into where is to be solved as The coefficient matrix K and F are given in Appendix B. If any point in a pulsed magnet yield (plastic deformation), (18) becomes nonlinear and must be solved by iteration methods based on the incremental plastic theory or the total strain plastic theory (Hencky-Theory) [32]. Since pulsed magnets are usually monotonically loaded during each shot, the total strain plastic theory is applied in this paper to obtain accurate results and high computing efficiency. In the total strain plastic theory, the constitutive matrix divides into two parts, elasticity D e and plasticity D p , as Accordingly, the stress balance relation (18) also divides into two parts as The elastic matrix K e is constant and the plastic matrix K p is variable. The iteration scheme given in (21) is applied to avoid repeatedly inverting the matrix and improve the computing efficiency. The iteration cycle stops when the relative error of between two iterations is considerably small (10 −3 in this paper).
The optimum locations of the separations are identified by iterative calculation with a decreasing number of independent CPs. Firstly, the maximum number k of possible separate CPs is derived from the distribution of the tensile radial stress calculated with no separation. As separations only occurs in the inner consecutive CPs, the remaining outer CPs are glued together. Then, calculate the displacement of the kth MIU, namely the outermost independent CP and the (k + 1)th MIU (Glued). If (15) is satisfied at the interface between the two MIUs, the optimum locations of the separations are found. Otherwise, stick the kth CP to the glued MIU and repeat the previous process with k = k − 1 until the criterion is met. Eventually, the optimum locations of the separations can be identified with the minimum calculations.

IV. RESULTS
A typical pulsed magnet with both internal and external reinforcement and the properties of the materials are given in TABLE 1. The magnet has a bore radius of 5.45 mm and a height of 98.7 mm. This magnet is made up of 10 layers of conductor and 13 layers of reinforcement. The windings locate between 5.45 mm and 67.5 mm of the radius. The Glass, stainless steel and carbon layers provide reliable external reinforcement for the magnet [23]. An average current density of 30×10 8 A/m 2 is applied on the conductor layers. The magnetic field distribution of the magnet is given in Fig. 6. It can generate 75.2 T at the center of the bore.

A. RESULTS OF THE MPS MODEL
The results of the MPS model are compared with that of a two-dimension axisymmetric model in the FEA software.  The stress distribution of the magnet is calculated with the classical GPS model, the MPS model and the FEA software respectively, with the first three CPs separated. The von Mises equivalent stress, hoop stress, radial stress and axial stress are shown in Fig. 7. The correspondent axial strain is shown in Fig. 8.
In the 1st∼3rd thin MIUs where the separations occur, each MIU is made up of one CP. The MPS model is equivalent to the GPS model because the axial strain is assumed to be uniform throughout the whole MIU in both cases. For these MIUs, the results of the MPS model and the GPS model meet quite well with the results of the FEA software as shown in Fig. 7. The axial strain is almost constant in the 1st∼3rd MIUs in Fig. 8, which confirms the validity of this assumption.
As to the thick MIU (4th∼10th CPs), the MPS model shows good agreement with the FEA software, while the results of the GPS model deviate substantially as shown in Fig. 7. In Fig. 8, the axial strain varies evidently through the thickness of the thick MIU due to the significant bending in the axial direction. The MPS model approximates the variation of the axial strain quite well with the discrete axial strain. Whereas, the constant axial strain assumption in the GPS model leads to an ''average'' strain over the whole MIU, such that the axial stress and strain is lower in the inner 4th∼9th CPs, and higher in the 10th CP as shown in Fig. 7(d) and Fig. 8. The lower axial strain of the GPS model couples through Poisson's ratio and results in positive radial stress (tensile) in the 4∼6th CPs in Fig. 7(c).
The radial tensile stress obtained with the GPS model suggests that more CPs should be separated to eliminate the stress concentration. However, due to the bending effect, the axial strain in the inner part of the thick MIU (4∼10th CPs) is actually higher than the ''average'' axial strain. The inner part of the thick MIU expands radially inwards under the axial compression. Therefore, some of the separate CPs calculated with the GPS model will contact with the swelled thick MIU and fail to separate. The incorrect judgment of the optimum locations of the separations will lead to much uncertainty. On the one hand, the excessively disparted CP fails to separate and it will squeeze the outer layer, increasing the stress in the outer CP. The actual stress in the outer CP will be larger than the calculation results. On the other hand, another CP is unnecessarily separated and becomes axially independent, which increases the risk of buckling due to the lower stability of thin cylinder windings.
The deviation between the MPS model and the FEA software grows lager as the radius increases in the outer region. This is mainly attributed to the error of redistribution of the axial force load and the constant axial strain assumption within each CP, which is inevitable in the one-dimension stress analysis on the mid-plane of solenoid magnets. Considering that the stress is usually low in the outer section, the simplicity and accuracy brought by the MPS model are highly valuable for efficient stress analysis of high field solenoid pulsed magnets.

B. OPTIMUM SEPARATION
Since orthotropic fiber has low transverse modulus, an effort has been made to enhance the lateral strength by compositing metal with fiber [33]. Nonetheless, isotropic materials are preferred for their uniform material properties in all directions. They have better performance in radial stress transfer and are more susceptible to the locations of the separations. Therefore, to illustrate the significance of the optimum separations, the stress is calculated assuming isotropic properties of the reinforcement, with E θ and ν rz as their equivalent isotropic Young's modulus E and Poisson's ratio ν.
The optimum positions of the separations between the CPs obtained with the proposed algorithm locate between the inner four CPs. The stress distribution calculated with no separation and four separate CPs respectively is shown in Fig. 9. When the CPs are all stuck together, positive radial stress (tensile) occurs in the inner six CPs, resulting in the stress concentration in the inner two CPs. By employing the free-standing technique on the inner four CPs, the peak von Mises stress is reduced from 4.5 GPa to 2.7 GPa.
The stress results of the magnet with different numbers of separate CPs are given in Fig. 10(a). If six CPs are set to separate according to the positive radial stress in Fig. 9, there will be a sharp rise of the stress in the excessive separate CPs (5th and 6th CP), with the peak von Mises stress of 3.1 GPa. This is because the separate boundaries impede the stress to be transferred to the outer CP. Fig. 10(b) shows that the 6th CP overlaps with the 7th CP, which indicates that contact occurs between them and the separate boundary is incorrect. Fewer separate units (three separate CPs) will lead to residue tensile stress and the peak stress of 2.8 GPa. Misjudgment of the locations of the separations will ultimately result in an increase of the peak stress.
Precisely identifying the locations of the optimum separations is critical to improving the stress distribution and seeking for more reliable mechanical structures. With the MPS model and the algorithm, the stress distribution and the optimum locations of the separations can be calculated efficiently and precisely.

V. CONCLUSION
The MPS model, which is based on the redistribution of the axial force load with the uniform axial plane stress assumption, can solve the stress distribution on the mid-plane of a solenoid pulsed magnet more accurately with the inclusion of notable bending effect. It also contributes to the precise identification of the optimum locations of the separations, which can be achieved with the proposed algorithm. The algorithm of identifying the optimum locations of separations not only improves the stress distribution in the magnet, but also avoids the contact analysis in the FEA model and brings higher computational efficiency. The proposed MPS model together with the algorithm achieves both high efficiency and high accuracy compared to the FEA model and the GPS model. It provides a more powerful tool for the mechanical analysis in the design and optimization of pulsed magnets where the bending effect is nonnegligible, such as ultra-high field pulsed magnets (above 50 T). These magnets are widely employed in physics, medical science, biology and many other areas, e.g. to study the superconductivity of FeSe [34] and the phase transitions in graphite [35]. The proposed model can also be embedded in an optimization algorithm to search for better structures and explore the potential of higher magnetic field strength (>100 T) with the existing materials.

APPENDIXES APPENDIX I. CALCULATION OF MAGNETIC FIELD
With the origin of the coordinates set at the center of the solenoid, a 1 and a 2 the inner and outer radius of the solenoid, z 1 and z 2 the minimum and maximum axial location, the boundaries of the solenoid magnet are normalized to the radius r of the field point (r, z) as A 1 = a 1 /r, A 2 = a 2 /r, Z 1 = z 1 /r, Z 2 = z 2 /r. (22) The intermediate function G(A, Z ) in the radial and axial direction are given in (23) and (24), as shown at the bottom of this page.

APPENDIX II. FINITE ELEMENT MODELING
According to the Galerkin-weighted residual method, we haveû whereû is an approximate solution, u j the solution for each node, R(û) the residual error between the approximate solution and the real solution, N j the interpolation basic function, the domain to be solved.