Multiscale thermoelastic modelling of composite strands using the “fundamental solutions” method

A novel computationally effective approach to multiscale thermoelastic modelling of composite structures and its application to a thermomechanical analysis of two ITER superconducting strands is presented. Homogenisation and recovering problems are solved by means of the “fundamental solutions” method, which was expanded to the case of thermoelastic analysis. We describe a general procedure of multiscale analysis on the basis of this method and apply it to recover stresses at the microscopic scale of a composite strands using a two-level procedure. The recovered micro-stresses are found to be in good correlation with the stresses obtained in the reference problem where the entire composite structure was modelled with a fine mesh.


Introduction
Practical impossibility of direct accounting for the entire microstructure when performing analysis of whole composite structures gave birth to the development of various multiscale analysis techniques. Such methods allow for consideration of the composite structure on different scales, and various kinds of techniques are used to ensure coupling between those scales [1,2]. One of the widespread approaches is based on two procedures: homogenisation and recovering. The first one, homogenisation, implies calculation of effective properties of a composite medium, which eventually allows to replace the heterogeneous medium with a homogenous one in the framework of a given analysis. First attempts to solve this problem go back to more than a century ago with the pioneering studies of Voigt and Reuss [3,4], and a number of different techniques have been introduced Dedicated to the memory of Prof. V.A. Palmov (1934Palmov ( -2018 1 since then. Among the most popular ones for periodic structures are those based on the asymptotic expansion method for differential equation with oscillating coefficients [5,6]. When finite element analysis is applied, solving boundaryvalue problems with periodic boundary conditions for a unit cell [7][8][9][10][11][12][13] has become arguably the most popular method for identifying effective mechanical properties of a composite medium. The ability to obtain effective properties of a composite medium and to perform the analysis on macroscopic scale (a homogenised solution) is, however, not always sufficient. The second procedure, recovering, which is also often referred to as unsmearing or heterogenization, is aimed at calculation of fields of interest on the microscopic scale starting from a homogenised solution. Typical applications of the recovering techniques are associated with situations where stress and strain states of individual components of a composite structure are required to account for non-linear behaviour. There are, however, problems where microscopic states affect the overall performance of the structure also in different ways. One of such examples is the composite superconducting strand of the ITER magnet coil considered below where superconductivity itself depends on the residual strain state.
ITER (International Thermonuclear Experimental Reactor) [14] is a scientific and engineering project unique in its complexity, which has raised many challenges. Despite the fact that its design is primarily completed and assembly of the machine is ongoing, looking for new ways to tackle problems raised by the ITER project is still relevant, especially in light of the development of forthcoming new fusion devices. Among the most challenging ITER components is its magnet system, which has significantly pushed the limits of current superconducting technology. It is known [15][16][17] that the critical current of Nb 3 Sn conductors, used for the generation of the strong toroidal magnetic field in ITER, substantially depends on the accumulated mechanical strains. The ability to predict strains in individual Nb 3 Sn filaments in a superconducting cable becomes, therefore, crucial in the design of superconducting coils. In view of the multiscale structure of the strand [18] (Fig. 1), it is necessary to obtain effective properties of such a structure, perform a macroscopic analysis for the magnetic system, step back down to microscopic scale and evaluate microscopic fields at this scale. The knowledge of residual strains after cool-down operation, needed for the operation of the coils is crucial for evaluating their performance. Various homogenisation and recovering techniques have been used by researchers to solve these problems under different assumptions and with different computational cost [19][20][21][22].
A novel method called "fundamental solutions and regular expansions" was proposed in 2006 by Palmov and Borovkov [23] for analysing stress and strain fields in periodic linear elastic heterogeneous media. The method is based on the idea that stress, strain, and displacement fields in a periodic medium can be represented as a linear combination of solutions of six problems for a unit cell. These solutions are called "fundamental solutions", and the linear combinations-"regular expansions". Beyond the inherent homogenisation abilities, the method features computationally effective recovering for the fields at microscopic scale. This paper presents a generalization of the "fundamental solutions and regular expansions" method to thermoelastic problems by adding the seventh "fundamental problem" and acquiring thereby the ability to recover stress and strain fields in a composite medium subjected to mechanical and thermal loads simultaneously. The superconducting strand of the ITER magnetic system described above presents a good example of a complex composite structure in which both types of loads are essential, and hence the application of the presented method to a multiscale analysis of the strand was used to test the method and to estimate its accuracy.

"Fundamental solutions" in thermoelasticity
We assume that a periodic composite structure is characterized by a three-dimensional unit cell, which has three planes of symmetry. We use local coordinates, with the coordinate planes parallel to the sides of the unit cell, and the origin located in the geometric centre: where h 1 , h 2 and h 3 are dimensions of the unit cell. Let us formulate seven fundamental boundary-value problems with simple mixed displacement-traction boundary conditions.

«Regular expansions» in thermoelasticity
The seven boundary-value problems (2)-(8) above have an important property-their solutions are periodic in terms of stresses. It can be proven that values of all stress tensor components at the corresponding points of opposite faces are equal. That can be done, for example, using proof by contradiction in the manner of [24] as it is demonstrated below for one of the normal components of the stress tensor.
Let us consider the first fundamental problem (2) for a symmetric unit cell, and let us assume that the normal stress σ 11 has a part that takes opposite values at the corresponding points of the opposite faces. Stress vectors produced by this part of σ 11 on the opposite faces are denoted with σ A 11 in Fig. 2a. Performing a reflection of the unit cell (together with all boundary conditions and internal forces) with respect to the vertical axis x 2 yields the system presented in Fig. 2b. The symmetry of the unit cell and of the boundary conditions (2) with respect to the vertical axis x 2 implies that the boundary value problems presented in Fig. 2a and b (after the reflection) are exactly the same. Opposite directions of σ A 11 in these two identical boundary value problems yields σ A 11 0, which means that σ 11 in the first fundamental problem (2) takes equal values at the two considered opposite faces of the unit cell.
Using similar arguments, it can be rigorously proven that all components of the stress tensor except for those prescribed as zeroes in the boundary conditions of (2)-(8), take equal values at the corresponding points of all opposite faces. The periodicity of stresses implies similar periodicity of strains and quasi-periodicity of displacements. Due to the linear nature of the equations and boundary conditions, any linear combination of the "fundamental solutions" can also be considered as a solution of elasticity equations for a periodic media- (9). Using the mentioned periodicity of the "fundamental solutions", we can consider their linear combination (9) an approximation for any periodic solution.
are the fields of stress tensor, strain tensor and displacement vector in the (k)-th fundamental problem, and λ (k) are constant multipliers (weight coefficients, which determine contribution of the (k)-th fundamental problem). We will refer to expressions (9) as "regular expansions" in thermoelasticity.
Let us introduce an operator of averaging across the unit cell volume < . . . > 1 V V dV and apply it to Eqs. (9). It can be easily shown that for some of the stress and strain tensor components in the "fundamental solutions", averaging results in zero values. For instance, for ε x x in the 4th fundamental problem (5): where only the values of u 1 at the faces of the unit cell are present. Using symmetry of the unit cell with respect to x 2 axis, it can be proven in similar manner as above that these values are bound to be equal. Indeed, a reflection with respect to x 2 axis results in the boundary value problem for the same unit cell but with opposite boundary conditions. This "reflected" boundary value problem, due to its linearity, must result in opposite displacements with respect to the initial one. That yields the conclusion that u 1 takes equal values at corresponding points of the faces x 1 0. Repeating the same reasoning for fundamental problems (6) and (7), we come to the conclusions that in all these problems associated with shear deformation of the unit cell, averaged values of the normal components of the strain tensor take zero values. For the first three fundamental problems, which are associated with stretching of the unit cell, similar reasoning results in zero values of the averaged non-diagonal components of the strain tensor. Finally, for the 7th fundamental problem (8), all averaged components of the strain tensor are zero. That allows to represent the result of averaging (9) in the matrix form as follows (γ ij 2ε ij are engineering shearing strains).
The averaged normal strains in the first six fundamental problems can be calculated analytically using the corresponding boundary conditions in the same way as demonstrated above for ε (4) 11 , which reveals that E is a diagonal matrix. The column of multipliers λ (k) can be calculated from the second expression in (10) as: Substitution of (14) into the first expression in (10) yields the following relation between averaged stresses and averaged strains in a periodic solution: which can be written in the following form: Equation (16) in fact represents a constitutive equation for averaged stresses and averaged strains in a thermoelastic medium: where C * · E −1 is the matrix of effective elastic moduli, which can be calculated on the basis of averaged solutions for the first six fundamental problems, and α * 1 , α * 2 , α * 3 are the effective thermal expansion coefficients, and T is the averaged change in temperature with respect to the reference one. (17) implies that the following substitution was made: On the other hand, the column of effective thermal expansion coefficients can be found by applying averaged constitutive Eqs. (17) to the (7)-th fundamental problem and taking into account zero values of all averaged strains in the solution: which yields: Knowing C * · E −1 and comparing Eqs. (18) and (20), we conclude that λ (7) T δ . Finally, the homogenisation problem in thermoelasticity can be solved by means of "fundamental solutions" method as:

Recovering micro stresses and micro strains by use of the «fundamental solutions» method
An important advantage of the "fundamental solutions" method is that not only it allows calculation of effective moduli according to (21), but it can also be used for recovering of micro fields in the homogenised solution, i.e., the solution of a macroscopic problem, in which a homogeneous material model with effective properties was used for the composite media. Using the column notations for non-averaged stress and stress tensors analogous to the one introduced in (11), we write present expansions (9) in the matrix form: Matrix notations and E in (22) are analogous to the ones introduced by (12)- (13), but are constructed for nonaveraged stresses and strains in the "fundamental solutions": (6) -the matrix of stresses in the first six fundamental problems; (6) -the matrix of strains in the first six fundamental problems.
It was shown above that the first six multipliers λ (k) can be calculated using the values of the strains averaged across the unit cell in the first six "fundamental solutions" (14). The 7th multiplier λ (7) can be identified using the averaged change in temperature T and the value of temperature applied in the 7 th fundamental problem as λ (7) T δ . By substitution of these expressions in (22) and using notations σ (7) and ε (7) for the columns of stress and strains components in the 7th fundamental problem, we obtain: ε E · E −1 · ε + T δ ε (7) It is necessary to note that Eq. (23) were obtained for perfectly periodic stresses and strains, which implies that strains averaged through the unit cell are constant across the whole composite medium. That is obviously not the case in most real applications, but we can expect that regular expansions (23) also provide a good estimate for stress and strain field in a composite if ε is variable, but smooth across the composite medium. Under that assumption, we can use (23) for recovering microscopic fields on the basis of a homogenised solution by substituting strains from that solution for the averaged across the unit cell strains: ε rec E · E −1 · ε homog + T δ ε (7) Note that in (24) we assume that strains from the homogenised solution ε homog are no longer constant across the cell of periodicity, but depend on coordinates. That obviously violates the condition of periodicity of strains and stress  (24). However, in view of the large ratio of macroscopic to microscopic dimensions in real composite structures, we can expect that even a non-homogenous distribution of strains in the homogenised solution leads to a small deviation of this field from a periodic one within a single unit cell. That provides grounds for using (24) for recovering microscopic stresses and stresses in real applications. With the purpose of estimating inaccuracy caused by the absence of periodicity in such applications, in the numerical examples presented below we use this method to recover microscopic stresses in regions of a composite structures that are relatively close to the boundary, and we compare results with those of the reference solutions obtained with a very fine mesh.

Recovering microscopic fields in three-scale thermoelastic problems
In this section we consider two thermomechanical problems for ITER superconducting strands. In the first problem, the ITER superconducting strand described briefly in the introduction is considered. As shown in Fig. 1, it is natural to identify several scales in the composite structure of the strand. Here we deal with three of them: the macroscopic scale M and two meso-scales m 1 and m 2 . Unit cells corresponding to these meso-scales and their finite element (FE) models (we use conventional finite element method and ANSYS code to obtain "fundamental solutions" as well as to solve the homogenised problem) are presented in Fig. 3. Our procedure is applied here in series as explained below.
With the aim to illustrate the application of the "fundamental solutions" method to thermomechanical analysis, the following problem was investigated. A partial cool down process is considered-the temperature of the strand is  [18]. Obviously, the strand is never used alone, but always as a part of a superconducting cable obtained with a twisting procedure to avoid eddy currents. To simulate the situation when a strand is in contact with three other strands, the strand is subjected to a distributed pressure at three regions- Fig. 4. The pressure is applied according to Hertz distribution: where s is a local coordinate for each "contact" area, with the origin at its centre, and a defines the size of the area (Fig. 4). In the considered numerical example, P max was taken 200 MPa. In real conditions, due to the complex multiscale twisting of strands inside the cable, contact areas are usually located in a non-symmetrical way [26], but they are supposed to be symmetric within the problem under consideration with the goal to have a self-balanced system. Material properties of individual components of the strand were assumed according to [22]. The problem was treated under plane strain assumption, which corresponds to the situation when the strand is straight and its ends are fixed. This is another assumption that is to some extent violated in reality, but radii of curvature of real strands in the cable are very large in comparison to its dimensions in a cross section [26], so we can consider that assumption as a reasonable approximation.
The key benefit of the plane strain assumption is the possibility to solve the same problem with a reference model where the entire composite structure is reproduced. This model, developed for testing the accuracy of the recovering method, is presented in Fig. 5. The model includes 4675 Nb 3 Sn filaments, consists of~450,000 s-order elements and has2 ,680,000 degrees of freedom. We will refer to the solution obtained with this model as the "reference solution" while the solution obtained for a model in which the whole composite area was represented with a homogenous material with effective characteristics, will be referred to as the homogenised solution.
The distribution of von Mises stresses in the cross section of the strand is presented in Fig. 6 for the reference solution. The stress field is due to both the mechanical loading and the cool down process.
The distribution of von Mises stresses in the reference solution (Fig. 6) reveals that whereas in the central zone of the cross section the stress field is nearly periodic, there are areas in the composite structure close to its boundaries where the periodicity condition is apparently violated. For the purpose to select reasonably difficult conditions for testing the recovering algorithm, we choose one of the unit cells located in such high-gradient area shown in Fig. 7.
The procedure of three-scale thermoelastic analysis for the composite strand using "fundamental solutions" method consists of the following steps.
Step 1 Homogenisation at m 1 scale.
We develop a finite element model of the unit cell associated with m 1 scale of Fig. 3a. We solve seven fundamental problems (2)-(8) for this model and, after calculating averaged values of stress and strain components in each problem, apply Eqs. (21) to obtain effective thermoelastic parameters at the meso-scale m 1 : We also store the distributions of stresses and strains in the seven fundamental problems for further use in the recovering of microscopic fields. Using the matrix notation introduced before, we can form the following matrices: Step 2 Transformation of the effective moduli between two mesoscopic scales. Boundary conditions in fundamental problems (2)- (8) are formulated in local coordinates associated with the unit cell. Therefore, due to the fact that unit cells of m 1 and m 2 meso-scales are oriented differently with respect to the strand, effective thermoelastic parameters C * m 1 , α * m 1 obtained at step 1 need to be translated into the local coordinates of the mese-scale m 2 . For the components of the fourthorder tensor C * m 1 the usual transformation takes the following form: C * (m 1 )ijkm l ip l jq l kr l ms C * (m 1 )pqrs , (26) where l ip are direction cosines of the angles between corresponding coordinate axes (ξ 1 ,η 1 ) of m 1 and (ξ 2 ,η 2 ) of m 2 scales, and the Einstein summation convention is used. A reduced form of the same transformation is applied for the diagonal secondorder tensor composed of the effective coefficients of thermal expansion α * m 1 . In the considered example, axes ξ 1 and ξ 2 have an angle of 73°.
Step 3 Homogenisation at m 2 scale.
We develop a finite element model of the unit cell associated with m 2 scale of Fig. 3b. For the region corresponding to the structure considered at m 1 scale we use the effective moduli obtained at step 2. We solve seven fundamental problems (2)-(8) for this model and, after calculating averaged values of stress and strain components, apply Eqs. (21) to obtain effective thermoelastic parameters at the meso-scale m 2 : C * m 2 , α * m 2 . We also store the distributions of stresses and strains in the seven fundamental problems for this scale: Step 4 Solving the homogenised problem at macroscopic scale.  We develop a finite element model of the structure at macroscopic scale, in which the entire composite structure is modelled with a homogeneous material with effective characteristics C * m 2 , α * m 2 (Fig. 8a). For the considered example, the macroscopic problem may be solved in the coordinates (x,y) which are identical to (ξ 2 ,η 2 ) used for the homogenisation at m 2 scale, so there is no need for a transformation of effective moduli obtained at step 3. We solve the macroscopic problem and obtain the distribution of stresses and strains within the domain σ hom , ε hom . The distribution of von Mises stresses in the cross section of the strand is presented in Fig. 8b to illustrate the "homogenised" solution.
Step 5 Recovering fields at m 2 scale.
Using Eqs. (24), we recover the distribution of stresses and strains at m 2 scale on the basis the homogenised solution and solutions of fundamental problems at m 2 scale: Step 6 Transformation of the recovered strains between two mesoscopic scales. The recovered strains according to (28) at m 2 scale are further used for recovering microscopic fields at the lowest scale m 1 . Due to the mentioned difference in orientation of the unit cells corresponding to these two mesoscopic scales, an appropriate transformation must be applied for ε rec(m 2 ) : ε rec (m 2 ) ij l ik l jm ε rec (m 2 ) km .
where the same notations as in (26) The presented procedure of three-scale analysis on the basis of "fundamental solutions" was implemented in APDL (Ansys Parametric Design Language) and used with the ANSYS FE code. A comparison of the recovered homogenised solutions and the reference solution for the axial stress σ z across the m 1 -scale unit cell shown in Fig. 7, is presented in Fig. 9. Large positive values of σ z are mostly caused by the cooldown from 80 to 4.2 K, but because of the difference in thermal expansion coefficients and in elastic moduli of Nb 3 Sn filaments and bronze matrix, there are "jumps" in σ z at the interface, which were successfully captured by the reference and recovered solutions. Contrary to that, the homogenised solution features approximately constant value of σ z across the interface, which is between the The difference between the recovered and reference solutions in this example does not exceed 2%, which is acceptable given that the considered unit cell is located relatively close to the boundary of the composite area in the strand (Fig. 7). However, both the recovered and reference solutions feature rather simple distribution of the axial stress, in which temperatureinduced stresses clearly dominate.
A significantly more difficult situation for the recovering algorithm is observed for some other components of the stress tensor that are strongly influenced by the pressure imitating contact with other strands. Figure 10 presents results of recovering σ x σ y σ xy component of the stress tensor in local coordinates of m 1 -scale as well as the von Mises stress along the same line in the same unit cell as in Fig. 9. The recovered solutions for all normal components of the stress tensor (Figs. 9, 10a and b) are fairly close to the reference solution (maximum deviation is observed for σ x (Fig. 10a), but is still less than 10%). The shear stress σ xy , however, appears to be recovered with much poorer accuracy in terms of relative deviation in % (Fig. 10c), which can be explained by small absolute values of this component-maximum values are roughly by an order of magnitude smaller than those for the normal components. That is illustrated by the comparison of recovered and reference values of the von Mises stress, which is important for predicting plastic deformations in such superconducting cables. In can be seen from Fig. 10d that in the considered example, the accuracy in recovering the von Mises stress is roughly as good as it is for the normal components of the stress tensor.
Overall, for the considered example the solution obtained using the two-scale recovering procedure provides a reasonably good approximation given that periodicity conditions are apparently poorly satisfied. Furthermore, it needs to be mentioned that the homogenised solution in case of σ x is clearly outside the range of reference and recovered values along the line. This can be explained by the three-scale nature of the composite structure of the strand combined with the choice of m 1 -scale unit cell located close to the boundary of the homogenised area within the m 2 -scale cell, which in turn is located close to the boundary of the homogenised area of the entire strand (Fig. 7). Nevertheless, the two-level recovering process based on "fundamental solutions" successfully captured that inhomogeneity, and the recovered σ x is fairly close to the reference one.
To make sure that the method works well for other composite structures and other loading options, we considered another superconducting ITER strand with significantly different microstructure. Figure 11 presents a photo of the cross section of the Luvata (LUV) strand [21,22] as well as the finite element models used for performing multiscale simulations and obtaining a reference solution. The LUV strand shown in Fig. 11 features a composite structure with less clear separation of scales-unit cells at both scales are not so small in comparison with the dimensions of the entire composite structure of the strand. Because of that reason, asymptotic homogenisation method is usually considered non applicable to such structure [21]. Nevertheless, with a purpose to test the proposed "fundamental solutions" approach in difficult conditions, the following problem was considered for the LUV strand. A partial cool down process when the temperature of the strand is decreasing from the cryostat temperature 80 K down to the operating temperature of the Nb 3 Sn superconductor 4.2 K is considered, and the strand is subjected to external pressure simulating contact with other strands. To make the pressure boundary conditions different from those used in the previous example, only two areas of contact were introduced- Fig. 12. The pressure is applied according to Hertz distribution with the maximum value P max 200 MPa. Because of the symmetrical nature of the problem, one quarter of the cross section region is considered. Figure 13 shows on the reference model two paths across the unit cells where the recovered solution is studied. Path 1 represents a unit cell that is located relatively far from the boundaries of the composite structure, whereas Path 2 represents a cell located in one of the areas of high stress gradient and almost at the boundary of the composite structure.
Comparison of the recovered and reference solutions for all components of the stress tensor and for von Mises stress along Path 1 is presented in Fig. 14. A fairly good correlation can be observed for all components, including the shear ones, and for the von Mises stress. A better accuracy in the recovering shear stresses in comparison with the previous example can be explained by a "better" location of Path 1, which is located relatively far from the boundaries of the composite structure and from the area where external load is applied.
To test the procedure in the worst possible within the considered problem scenario, a recovering along Path 2 was also studied. This path (Fig. 13) is located very close to the boundary of the composite structure and right under the area of pressure application-that makes gradients of stresses in this zone particularly large, and the periodicity assumption, Fig. 12 A schematic of the pressure application used in the simulation for the LUV strand, also called OUK strand in [21] which is the foundation of the method, is therefore poorly satisfied.
Comparison of the recovered and reference solutions for the stresses along Path 2 is presented in Fig. 15. The accuracy of the recovering decreased noticeably in comparison to what was observed on Path 1 for all components of the stress tensor except for σ z , which is mostly determined by thermal stresses and not so much affected by the applied pressure. However, the accuracy in recovering the von Mises stress is still about 5%, which may be considered a good correlation and indicates applicability of the method even in situations when the solutions are not exactly periodic. To highlight the gradients in the stress field in the considered area, a contour plot of the von Mises stress in the reference solution is shown in Fig. 16.

Conclusions
An extension of the "fundamental solutions" method in the mechanics of periodic composites to thermoelastic problems is presented and its accuracy is studied using two different composite structures. The method allows for solving homogenisation and recovering problems in multiscale thermomechanical analysis of composite structures characterized by linear behaviour of their components.
A three-scale analysis was performed using this method for two different superconducting strands of the ITER magnet system, both of which feature a complex composite structure. A problem of the strand cool down from the cryostat temperature 80 K down to the working conditions 4.2 K with simultaneous side compression from other strands was simulated. Models of the strands with a very fine mesh, in which the entire composite structure was represented, were used to obtain reference solutions for controlling the accuracy of  the recovered micro fields. The comparison revealed that the accuracy of the method depends on the extent, to which the distribution of stresses and strains in the composite structure can be considered periodic. In the considered examples, the recovered major stress tensor components were within 10% from the ones provided by the reference solution even in the most difficult conditions. The von Mises stress was found to be recovered with even better accuracy of about 5% even in the areas located close to the boundary of the composite structure. Since in the considered linear case stresses and strains are related by the generalised Hooke's lay, one may expect that recovering the various strain functions used to describe the dependence of critical current superconducting cables on strains [27] may be carried out with good accuracy using the proposed method.
The method only requires seven fundamental problems to be solved once for unit cells representing each mesoscopic scale of a composite structure, and the recovery of microscopic fields is performed on the basis of algebraic calculations involving these solutions of fundamental problems and a homogenised solution on the macroscopic scale. Thereby, the method can be considered as computationally not expensive, whereas the main limitation of its application follows from the underlying idea of "regular expansions", which requires linear behaviour of the composite components. However, an explicit incremental solution as in [20] may allow to include some nonlinearity if the nature of the problems allows considering the materials behaviour linear elastic within each increment. In that case the fundamental problems have to be solved at each step making the procedure more expensive. This will be further investigated.