Test Case : Feelpp CFPDES - Stacked tapes

1. Introduction

This is the test case of High-Temperature Superconductors with transport current using the T-A Homogeneous Formulation on stacked tapes geometry surrounded by air in axisymmetric coordinates.

2. Run the Calculation

The command line to run this case is :

    mpirun -np 16 feelpp_toolbox_coefficientformpdes --config-file=mqs_axis.cfg

3. Data Files

The case data files are available in Github here :

CFG file - Edit the file
JSON file - Edit the file
GEO file - Edit the file

4. Equation

The T-A Homogeneous Formulation in axysimmetric coordinates is :

T-A Formulation in axisymmetric coordinates

\[\text{(AV 2D)} \left\{ \begin{matrix} -\frac{1}{\mu}\Delta A_\theta + \frac{1}{\mu r^2}A_\theta &=& \frac{\delta}{\Lambda}\partial_z T_r & \text{on } \Omega^{axis} & \text{(A Axis)} \\ \begin{pmatrix} 0 & 0\\ 0 & -\rho \end{pmatrix}\Delta T_r &=& \frac{\partial(\partial_z A_\theta)}{\partial t} & \text{on } \Omega_c^{axis} & \text{(T Axis)} \end{matrix} \right.\]

With :

\(A_\theta\) : \(\theta\) component of potential magnetic field
\(T_r\) : \(r\) component of potential current density
\(\rho\) : electric resistivity \(\Omega \cdot m\)defined by the E-J power law on HTS : \(\rho=\frac{E_c}{J_c}\left(\frac{\mid\mid J \mid\mid}{J_c}\right)^{(n)}=\frac{E_c}{J_c}\left(\frac{\mid\mid \partial_z T_r \mid\mid}{J_c}\right)^{(n)}\)
\(\mu\) : electric permeability \(kg/A^2/S^2\)
\(\delta\) : thickness of the tapes \(m\)
\(\Lambda\) : space between the Tapes \(m\)

5. Geometry

The geometry is a set of stacked tapes in axisymmetric coordinates \((r,z)\), surrounded by air.

Geometry

Zoom on Tapes

The geometrical domains are :

Bulks : the tapes represented by bulks for the Homogeneous approach
- Top : the top boundary of the tapes
- Bottom : the bottom boundary of the tapes
- Sides : the side boundaries of the tapes
Air : the air surrounding Conductor
- Symmetry Line : the boundary for the symmetry in axisymmetric coordinates
Spherical Shell : the air surrounding Conductor
- Exterior boundary : the Spherical Shell 's boundary

Symbol

Description

value

unit

\(W_{tapes} (\delta)\)

tapes width

\(1e-6\)

\(H_{tapes}\)

tapes height

\(12e-3\)

\(cel_{tapes} (\Lambda)\)

space between the tapes

\(250e-6\)

\(R_{inf}\)

radius of infty border

\(0.4\)

6. Parameters

The parameters of the problem are :

On Bulks :

Symbol

Description

Value

Unit

\(\mu=\mu_0\)

magnetic permeability of vacuum

\(4\pi.10^{-7}\)

\(kg \, m / A^2 / S^2\)

\(thickness_{tape}\)

tapes width

\(1e-6\)

\(m\)

\(thickness_{cell}\)

space between the tapes

\(250e-6\)

\(m\)

\(height\)

tapes height

\(12e-3\)

\(m\)

\(f\)

frequency

\(50\)

\(Hz\)

\(Ic\)

critical current

\(300\)

\(A\)

\(Iapp\)

applied current

\(0.5*Ic*sin(2*\pi*f*t)\)

\(A\)

\(j_c\)

critical current density

\(2.5.10^10\)

\(A/m^2\)

\(e_c\)

threshold electric field

\(10^{-4}\)

\(V/m\)

\(n\)

material dependent exponent

\(25\)

\(\rho\)

electrical resistivity (described by the \(e-j\) power law)

\(\frac{e_c}{j_c}\left(\frac{\mid\mid \partial_z T_r \mid\mid}{j_c}\right)^{(n)}\)

\(\Omega\cdot m\)

On Air :

Symbol

Description

Value

Unit

\(\mu=\mu_0\)

magnetic permeability of vacuum

\(4\pi.10^{-7}\)

\(kg \, m / A^2 / S^2\)

On JSON file, the parameters are written :

Parameters on JSON file

    "Parameters":
    {
        "timestep":2e-4,

        "thickness_tape":1e-6,   //delta
        "thickness_cell":250e-6, //Lambda
        "height":12e-3,

        "mu":"4*pi*1e-7",
        "f":50,
        "Ic":300.0,
        "Iapp":"0.5*Ic*sin(2*pi*f*t):Ic:f:t",

        "jc":"Ic/thickness_tape/height:Ic:height:thickness_tape", //=2.5E10
        "ec":1e-4,
        "n":25
    },
    "Materials":
    {
        "Conductor":
        {
            "markers":["Bulks"],

            "rhoHTS":"ec/jc*((abs(current_grad_T_rt_1)/jc)^(n)):ec:jc:n:current_grad_T_rt_1",

            "J":"current_grad_T_1*(thickness_tape/thickness_cell):thickness_tape:thickness_cell:current_grad_T_1"
        },
        "Air":
        {
	        "markers":["Air","Spherical_shell"]
        }
    },

7. Boundary Conditions

For the Dirichlet boundary conditions, we want to impose the transport current :

The transport current is the difference of current potential between the top and the bottom of the tapes divided by the thickness of the tapes, so we can impose \(0\) at the bottom and \(Iapp/thickness\) at the top of the tapes

Finally we have :

On Top : \(T_r = Iapp/\delta\)
On Bottom : \(T_r = 0\)

For the Neumann boundary conditions, the homogeneous approach of the T-A formulation impose that it is \(0\) on both side of the bulk representation.

Finally we have :

On Sides : \(T_r \cdot \mathbf{n} = 0\)

On JSON file, the boundary conditions are written :

Boundary conditions on JSON file

    "BoundaryConditions":
    {
        "current":
        {
            "Dirichlet":
            {
                "Top":
                {
                    "expr":"Iapp/thickness_tape:thickness_tape:Iapp"
                },
                "Bottom":
                {
                    "expr":"0"
                }
            },
            "Neumann":
            {
                "Sides":
                {
                    "expr":0
                }
            }
        }
    },

8. Weak Formulation

Weak formulation of T-A Formulation in axisymmetric coordinates

\[\scriptsize{ \text{(Weak TA Hom Axis)} \left\{ \begin{eqnarray*} \int_{\Omega^{axis}}{ \frac{r}{\mu} \tilde\nabla A_\theta \cdot \tilde\nabla \phi \ drdz } + \int_{\Omega^{axis}}{\frac{1}{\mu r} A_\theta \cdot \phi \ drdz} = \int_{\Omega^{axis}_c}{ \frac{\delta}{\Lambda}\frac{\partial T_r}{\partial z} \cdot\phi \ rdrdz } \\ \int_{\Omega^{axis}_c}{ \begin{pmatrix} 0 & 0 \\ 0 & \rho \end{pmatrix} \tilde\nabla T_r \cdot \tilde\nabla \varphi \ drdz } = \int_{\Omega^{axis}_c}{ \frac{\partial (\partial_z A_\theta)}{\partial t} \cdot\varphi \ rdrdz }\\ \text{for } \phi \in H^1(\Omega) \text{ and } \varphi \in H^1(\Omega_c) \end{eqnarray*} \right. }\]

9. Coefficient Form PDEs

The Feelpp toolboxe Coefficient Form PDEs is used here. The coefficients associated to the Weak Formulation are :

On Bulks :
- A Formulation

Coefficient

Description

Expression

\(c\)

diffusion coefficient

\(\frac{r}{\mu}\)

\(a\)

absorption or reaction coefficient

\(\frac{1}{\mu r}\)

\(f\)

source term

\(r\frac{\delta}{Lambda}\partial_z T_r\)

T Formulation

Coefficient

Description

Expression

\(c\)

diffusion coefficient

\(\begin{pmatrix} 0 & 0\\ 0 & -r\rho \end{pmatrix}\)

\(f\)

source term

\(r\frac{\partial(\partial_z \textbf{A})}{\partial t}\)

On Air :
- A Formulation

Coefficient

Description

Expression

\(c\)

diffusion coefficient

\(\frac{r}{\mu}\)

\(a\)

absorption or reaction coefficient

\(\frac{1}{\mu r}\)

On JSON file, the coefficients are written :

CFPDEs coefficients on JSON file

    "Models":
    {
        "cfpdes":{
            "equations":["magnetic","current"]
        },
        "magnetic":{
            "common":{
                "setup":{
                    "unknown":
                    {
                        "basis":"Pch1",
                        "name":"A",
                        "symbol":"A"
                    }
                }
            },
            "models":[{
                "name":"magnetic_Conductor",
                "materials":"Conductor",
                "setup":{
                    "coefficients":{
                        "c":"x/mu:x:mu",
                        "a":"1/mu/x:mu:x",
                        "f":"x*materials_Conductor_J:x:materials_Conductor_J"
                    }
                }
            },{
                "name":"magnetic_Air",
                "materials":"Air",
                "setup":{
                    "coefficients":{
                        "c":"x/mu:x:mu",
                        "a":"1/mu/x:mu:x"
                    }
                }
            }]
        },
        "current":{
            "common":{
                "setup":{
                    "unknown":
                    {
                        "basis":"Pch1",
                        "name":"T",
                        "symbol":"T"
                    }
                }
            },
            "models":[{
                "name":"current_Conductor",
                "materials":"Conductor",
                "setup":{
                    "coefficients":{
                        "c":"{0,0,0,x*materials_Conductor_rhoHTS}:x:materials_Conductor_rhoHTS",
                        "f":"(magnetic_grad_A_1-magnetic_grad_A_previous_1)*x/timestep:x:magnetic_grad_A_1:magnetic_grad_A_previous_1:timestep"
                    }
                }
            }]
        }
    },

10. Numeric Parameters & Mesh

Time
- Initial Time : \(0s\)
- Final Time : \(0.02s\)
- Time Step : \(2e-4s\)
Mesh :

Mesh of Geometry

11. Solving the Model

In order to solve the non-linearity in the model, several parameters are used in the json and the cfg file :

The non-linear solver used is Newton, more efficient when \(\rho\) is the parameter that contain the e-j power law.

in CFG file

solver=Newton

The type of preconditioner chose is LU :

in CFG file

pc-type=lu
ksp-type=preonly

The maximum number of iteration for the non-linear solver is fixed at 200

in CFG file

snes-maxit=200

The number of iteration is usually low except when the applied current change of direction of variation :

Number of iteration compared to the evolution of the applied current

12. Results

The results that we obtain with this formulation with Feelpp are compared to the results of the article Real-time simulation of large-scale HTS systems : multi-scale and homogeneous models using the T-A formulation where the software Comsol is used.

The time evolution of the applied current is :

Time evolution of the applied current

12.1. Electric current density

The electric current density \(j_\theta\) is defined by :

\[ j_\theta=\frac{\partial T}{\partial z}\]

Electric current density \(j_\theta (A/m^2)\)

We compare the distribution of the electric current density on the Oz axis between the tapes at the instant \(t=0.005s\) with Feelpp and Comsol.

L2 Relative Error Norm : \(9.98 \%\)

12.2. Magnetic flux density

The magnetic flux density \(B\) is defined by:

\[ B=\nabla\times A =\begin{pmatrix}-\partial_z A_\theta\\ 0\\ \frac{1}{r}\partial_r (rA_\theta)\end{pmatrix}\]

Magnitude of the Magnetic flux density \(B (T)\)

As such, \(B_r=-\partial_z A_\theta\) and \(B_z=\frac{1}{r}\partial_r (rA_\theta)\) :

r_component of the magnetic flux density \(B_r (T)\)

z_component of the magnetic flux density \(B_z (T)\)

We compare the distribution of the r-component of the magnetic flux density on the Oz axis between the tapes at the instant \(t=0.005s\) with Feelpp and Comsol.

L2 Relative Error Norm : \(2.86 \%\)

We also compare the distribution of the z-component of the magnetic flux density on the Or axis across the tapes at the instants \(t=0.005s\) with Feelpp and Comsol.

L2 Relative Error Norm : \(3.02 \%\)

13. References

Real-time simulation of large-scale HTS systems: multi-scale and homogeneous models using the T–A formulation, Edgar Berrospe-Juarez et al 2019 Supercond. Sci. Technol. 32 065003, PDF