Test Case : Feelpp CFPDEs - Bulk Cylinder in Axisymmetric coordinates

1. Introduction

This is the test case of High-Temperature Superconductors with the A-V Formulation and Gauge Condition on a cylinder geometry surrounded by air in axisymmetric coordinates. The formulation used here is the Formulation 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 --cfpdes.solver=Picard --cfpdes.verbose_solvertimer=1

This simulation uses the last version 109 of Feelpp.

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 A Formulation in axisymmetric coordinates is :

A Formulation in axisymmetric coordinates

\[ -\frac{1}{\mu}\Delta A_\theta + \frac{1}{\mu r^2}A_\theta + \sigma \frac{\partial A_\theta}{\partial t} =0 \text{on } \Omega^{axis} \text{(A Axis)}\]

With :

\(A_{\theta}\) : \(\theta\) component of potential magnetic field
\(\sigma\) : electric conductivity \(S/m\) described by the e-j power law : \(\sigma=\frac{J_c}{E_c}\left(\frac{\mid\mid e\mid\mid}{E_c}\right)^{(1-n)/n} = \frac{J_c}{E_c}\left(\frac{\mid\mid -\partial_t A_\theta\mid\mid}{E_c}\right)^{(1-n)/n}\)
\(\mu\) : electric permeability \(kg/A^2/S^2\)

5. Geometry

The geometry is a rectangle in axisymmetric coordinates \((r,z)\) representing a bulk cylinder, surrounded by air.

Geometry in Axisymmetrical cut

The geometrical domains are :

Cylinder : the cylinder, composed by a conductor
Air & Spherical Shell : the air surrounding Conductor
- Symmetry Line : Air 's bound, correspond to \(Oz\) axis (\(\{(z,r), \, z=0 \}\))
- Exterior Boundary : the rest of the Air 's bound

Symbol

Description

value

unit

\(R\)

radius of cylinder

\(0.00125\)

\(H_{cylinder}\)

height of cylinder

\(0.001\)

\(R_{inf}\)

radius of infty border

\(0.1\)

6. Parameters

The parameters of the problem are :

On Cylinder :

Symbol

Description

Value

Unit

\(\mu=\mu_0\)

magnetic permeability of vacuum

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

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

\(t_f\)

final time

\(15\)

\(s\)

\(b_{max}\)

maximum applied field

\(1\)

\(T\)

\(rate\)

rate of the applied field raise

\(\frac{3}{tf}b_{max}\)

\(T/s\)

\(hsVal\)

time ramp for the applied field

\(\begin{cases}rate*t &\quad\text{if }t<\frac{t_f}{3}\\b_{max} &\quad\text{if }t<\frac{2t_f}{3}\\b_{max} - (t-\frac{2t_f}{3})*rate &\quad\text{if }t>\frac{2t_f}{3}\end{cases}\)

\(T\)

\(J_c\)

critical current density

\(3.10^8\)

\(A/m^2\)

\(E_c\)

threshold electric field

\(10^{-4}\)

\(V/m\)

\(n\)

material dependent exponent

\(20\)

\(\sigma\)

electrical conductivity (described by the \(E-J\) power law)

\(\frac{J_c}{E_c}\left(\frac{\mid\mid e\mid\mid}{E_c}\right)^{(1-n)/n}\)

\(S/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":
    {
        "mu":"4*pi*1e-7",
        "tf":15,
        "bmax":1.0,
        "rate":"3.0/tf*bmax:tf:bmax",
        "hsVal":"bmax - (t-2.0*tf/3.0)*rate + (t<2.0*tf/3.0)*(t-2.0*tf/3.0)*rate + (t<tf/3.0)*(t*rate - bmax):t:tf:rate:bmax",

        "jc":3e8,
        "ec":1e-4,
        "n":20,
        "epsSigma":1e-8
    },
    "Materials":
    {
        "Conductor":
        {
            "sigma":"jc / ec * 1.0 / ( epsSigma + ( abs(-magnetic_dAtheta_dt)/ec )^((n-1.0)/n) ):jc:ec:n:epsSigma:magnetic_dAtheta_dt",
            "sigma_mod":"x*jc / ec * 1.0 / ( epsSigma + ( abs(-magnetic_dAtheta_dt)/ec )^((n-1.0)/n) ):x:jc:ec:n:epsSigma:magnetic_dAtheta_dt",
            \\ [...]
        },
    },

7. Boundary Conditions

For the Dirichlet boundary conditions, we want to impose the applied magnetic field :

So we have \(B=hsVal\), therefore \(\nabla\times A = hsVal\) and so \(\frac{1}{r}\partial_r (rA_\theta)=hsVal\).

Finally we have :

On Symmetry Line & Exterior Boundary : \(A_{\theta} = \frac{r}{2}hsVal\)

On JSON file, the boundary conditions are written :

Boundary conditions on JSON file

    "BoundaryConditions":
    {
        "magnetic":
        {
            "Dirichlet":
            {
                "magdir":
                {
                    "markers":["Symmetry_line","Exterior_boundary"],
                    "expr":"x/2 *hsVal:x:hsVal"
                }
            }
        }
    },

8. Weak Formulation

Weak formulation of A Formulation in two dimensions approximation

\[ \int_{\Omega_c^{axis}}{\sigma \frac{\partial A_\theta}{\partial t} \, \phi \ rdrdz } + \int_{\Omega^{axis}}{ \frac{1}{\mu} \nabla A_\theta \cdot \nabla \phi \ rdrdz }+\int_{\Omega^{axis}}{ \frac{1}{\mu r} A_\theta \cdot \phi \ drdz } = 0 \hspace{1cm} \text{(Weak A Axi)}\]

9. Coefficient Form PDEs

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

On Cylinder :

Coefficient

Description

Expression

\(d\)

damping or mass coefficient

\(\sigma r\)

\(c\)

diffusion coefficient

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

\(a\)

absorption or reaction coefficient

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

\(f\)

source term

\(0\)

On Air :

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

    "Materials":
    {
        "Conductor":
        {
            "markers":["Cylinder"],
            "sigma":"jc / ec * 1.0 / ( epsSigma + ( abs(-magnetic_dAtheta_dt)/ec )^((n-1.0)/n) ):jc:ec:n:epsSigma:magnetic_dAtheta_dt",
            "sigma_mod":"x*jc / ec * 1.0 / ( epsSigma + ( abs(-magnetic_dAtheta_dt)/ec )^((n-1.0)/n) ):x:jc:ec:n:epsSigma:magnetic_dAtheta_dt",

	        "magnetic_c":"x/mu:x:mu",
            "magnetic_a":"1/mu/x:mu:x",
            "magnetic_f":"0.",
            "magnetic_d":"sigma_mod:sigma_mod"
        },
        "Air":
        {
	        "markers":["Air","Spherical_shell"],
            "magnetic_c":"x/mu:x:mu",
            "magnetic_a":"1/mu/x:mu:x"
        }
    },

10. Numeric Parameters

Time
- Initial Time : \(0s\)
- Final Time : \(15s\)
- Time Step : \(1s\)
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 Picard, more efficient when \(\sigma\) is the parameter that contain the e-j power law.

in CFG file

solver=Picard

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

in CFG file

snes-maxit=600

The time-step is high to help the convergence of the model :

in CFG file

time-step=1

The number of iterations ranges from 80 to 120 :

Number of iteration compared to the evolution of the applied magnetic field

12. Results

The results that we obtain with this formulation with Feelpp are compared to the results of the article Finite-Element Formulations for Systems With high-temperature Superconductors where the solver getDP is used.

The time evolution of the applied field is :

Time evolution of the external applied field

With \(t_1=5s\), \(t_2=10s\) and \(t_3=15s\)

12.1. Electric current density

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

\[ J_\theta=\sigma E = -\sigma \frac{\partial A_\theta}{\partial t}\]

With :

\[ \sigma = \frac{J_c}{E_c}\left(\frac{\mid\mid e\mid\mid}{E_c}\right)^{(1-n)/n}= \frac{J_c}{E_c}\left(\frac{\mid\mid -partial_t A_\theta\mid\mid}{E_c}\right)^{(1-n)/n}\]

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

We compare the current density profiles with Feelpp and getDP on the \(O_r\) axis, at the mid-height of the cylinder, at time \(t_3\) for a maximum applied field of 1 T and \(n=20\).

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

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}\]

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

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

We compare the distribution of the z-component of the magnetic flux density 2mm above the cylinder at the instants \(t_1\), \(t_2\) and \(t_3\) with Feelpp and getDP.

t1 \(=5s\)

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

t2 \(=10s\)

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

t3 \(=15s\)

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

13. References

Finite-Element Formulation for Systems with High-Temperature Superconductors, Julien Dular, Christophe Gauzaine, Benoît Vanderheyden, IEEE Transactions on Applied Superconductivity VOL. 30 NO. 3, April 2020, PDF