Test case of Thermo-Magnetism with Linear Coefficients in Axisymmetrical case

1. Introduction

This page presents the simulation of electromagnetism in A-V Formulation coupled with thermic problem on geometry of torus in transient case in axisymmetrical case with CFPDEs application.

2. Run the calculation

The command line to run this case is :

    mpirun -np 4 feelpp_toolbox_coefficientformpdes --config-file=thermo-mag.cfg --cfpdes.gmsh.hsize=1e-3

To run the case with adaptative time-step :

    sh feelpp_adaptative.sh

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
SH script for adaptative time-step - Edit the file

4. Equation

Differential Formulation in Axisymmetric

\[\text{(AV+Heat Axis)} \left\{ \begin{matrix} \sigma \frac{\partial \Phi}{\partial t} - \frac{1}{\mu} \Delta \Phi + \frac{2}{\mu \, r} \frac{\partial \Phi}{\partial r} + \sigma \frac{\partial V}{\partial \theta} = 0 & \text{ on } \Omega & \text{(AV-1 Axis)} \\ V = \frac{U}{2\pi} \, \theta & \text{ on } \Omega_c^{axis} & \text{(AV-2 Axis)} \\ A_{\theta} = 0 & \text{ on } \Gamma_D^{axis} & \text{(D Axis)} \\ \frac{\partial A_{\theta}}{\partial \mathbf{n}^{axis}} = 0 & \text{ on } \Gamma_N^{axis} & \text{(N Axis)} \\ \rho C_p \, \frac{\partial T}{\partial t} - k \Delta T = \sigma \, \left( \frac{1}{r} \, \frac{\partial V}{\partial \theta} + \frac{\partial A_{\theta}}{\partial t} \right)^2 & \text{ on } \Omega_c^{axis} & \text{(Heat Axis)} \\ \frac{\partial T}{\partial \mathbf{n}^{axis}} = 0 & \text{ on } \Gamma_{Nc}^{axis} & \text{(Nc Axis)} \\ - k \, \frac{\partial T}{\partial \mathbf{n}^{axis}} = h \, \left( T - T_c \right) & \text{ on } \Gamma_{Rc}^{axis} & \text{(Rc Axis)} \end{matrix} \right.\]

With \(\Phi = r \, A_{\theta}\)

5. Geometry

The geometry is a torus of the conductor in cartesian coordinates \((x,y,z)\) or rectangle in axisymmetric coordinates \((r,z)\), surrounded by air.

Geometry in axisymmetrical cut loop on Conductor

Geometry in axisymmetrical cut

Geometry in three dimensions

Geometry in three dimensions cut loop on Conductor

The geometrical domains are :

Conductor : the torus is composed by conductor material
- Interior : interior surface of conductor ring
- Exterior : exterior surface of conductor ring
- Upper : upper of conductor ring
- Bottom : bottom of conductor ring
Air : the air surround Conductor
- zAxis : a bound of Air, correspond to \(Oz\) axis (\(\{(z,r), \, z=0 \}\))
- infty : the rest of bound of Air

Symbol

Description

value

unit

\(r_{int}\)

interior radius of torus

\(75e-3\)

\(r_{ext}\)

exterior radius of torus

\(100.2e-3\)

\(z_1\)

half-height of torus

\(25e-3\)

6. Initial/Boundary Conditions

We impose the boundary conditions :

Dirichlet :
- On zAxis : \(\Phi = 0\) (\(A_{\theta} = 0\) by symetric argument)
- On infty : \(\Phi = 0\) (\(A_{\theta} = 0\) we consider the bound of resolution like infty for magnetic field)
Neumann : \(\frac{\partial T}{\partial \mathbf{n}} = 0\) on Interior and Exterior
Robin : \(-k \, \frac{\partial T}{\partial \mathbf{n}} = h \, \left( T - T_c \right)\) on Upper and Bottom

We initialize :

On Conductor \(\Phi(t=0,r,z) = 0\) (\(A_{\theta}(t=0,r,z) = 0\)) and \(T(t=0,r,z) = T_i\)
On Air \(\Phi(t=0,r,z) = 0\) (\(A_{\theta}(t=0,r,z) = 0\))

On JSON file, the boundary conditions are writed :

Boundary conditions on JSON file

    "BoundaryConditions":
    {
        "magnetic":
        {
            "Dirichlet":
            {
                "ZAxis":
                {
                    "expr":"0"
                },
                "Infty":
                {
                    "expr":"0"
                }
            }
        },
        "heat":
        {
            "Robin":
            {
                "Interior":
                {
                    "expr1":"h*x:h:x",
                    "expr2":"h*T_c*x:h:T_c:x"
                },
                "Exterior":
                {
                    "expr1":"h*x:h:x",
                    "expr2":"h*T_c*x:h:T_c:x"
                }
            }
        }
    }

On JSON file, the intial conditions are writed : .Initial conditions on JSON file

    "InitialConditions":
    {
        "temperature":
        {
            "Expression":
            {
                "myic":
                {
                    "markers":"Conductor",
                    "expr":"T_i:T_i"
                }
            }
        }
    }

7. Weak Formulation

We obtain :

Weak Formulation

\[\scriptsize{ \text{(Weak Thermo-Mag Axis)} \begin{eqnarray*} \int_{\Omega^{axis}}{ \sigma \frac{\partial \Phi}{\partial t} \mathbf{ϕ} \ r\ drdz} + \int_{\Omega^{axis}}{ \frac{1}{\mu} \tilde{\nabla} \Phi \cdot \tilde{\nabla} \mathbf{ϕ} \ r \ drdz } + \int_{\Omega^{axis}}{\frac{2}{\mu \, r} \frac{\partial \Phi}{\partial r} \, \mathbf{ϕ} \ r \ drdz} &=& - \int_{\Omega_c^{axis}}{ \sigma \frac{\partial V}{\partial \theta} \mathbf{ϕ} \ r \ drdz } & \text{(Weak AV Axis)} \\ \int_{\Omega_c^{axis}}{ \rho C_p \, \frac{\partial T}{\partial t} \, \mathbf{ϕ} \, r \, drdz } + \int_{\Omega_c^{axis}}{ k \, \tilde{\nabla} T \cdot \tilde{\nabla} \mathbf{ϕ} \, r \, drdz } + \int_{\Gamma_R^{axis}}{ h \, T \, \mathbf{ϕ} \, r \, d\Gamma } &=& \int_{\Omega_c^{axis}}{ Q \, \mathbf{ϕ} \, r \, drdz } + \int_{\Gamma_R^{axis}}{ h \, T_c \, \mathbf{ϕ} \, r \, d\Gamma } & \text{(Weak Heat Axis)} \end{eqnarray*} }\]

With \(\tilde{\nabla} = \begin{pmatrix} \frac{\partial}{\partial r} \\ \frac{\partial}{\partial z} \end{pmatrix}\)

8. Parameters

The parameters of problem are :

On Conductor :

Symbol

Description

Value

Unit

\(V\)

scalar electrical potential

\( U \, \frac{\theta}{2\pi}\)

\(Volt\)

\(U\)

electrical potential

\(\left\{ \begin{matrix} t \text{ on } 0s \leq t \leq 1s \\ 1 \text{ on } 1s \leq t \leq 20s \\ 0 \text{ on } 20s \leq t \leq 22s \end{matrix} \right.\)

\(Volt\)

\(\sigma\)

electrical conductivity

\(4.8e7\)

\(S/m\)

\(\mu=\mu_0\)

magnetic permeability of vacuum

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

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

\(k\)

thermal conductivity

\(380\)

\(W/m/K\)

\(C_p\)

thermal capacity

\(380\)

\(J/K/kg\)

\(\rho\)

density

\(10000\)

\(kg/m^3\)

\(h\)

convective coefficient

\(80000\)

\(W/m^2/K\)

\(T_c\)

cooling temperature

\(293\)

\(K\)

\(T_i\)

initial temperature

\(293\)

\(K\)

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 parmeters are writed :

Parameters on JSON file

    "Parameters":
    {
        "h":80000,  // W*m-2*K-1
        "T_c":293,  // K
        "T_i":293,  // K
        "U":"t/1.*(t<1.)+(t<20.)*(t>1.):t", // Volt

        // Constants of analytical solve
        "a":1933.10,    // K
        "b":0.40041,    // K
        "rmax":0.0861910719118454,  // m
        "Tmax":364.446  // K
    }

9. Coefficient Form PDEs

We use the application Coefficient Form PDEs. The coefficient associate to Weak Formulation are :

For MQS equation (Weak MQS Axis) :

On Conductor :

Coefficient

Description

Expression

\(d\)

damping or mass coefficient

\(\sigma\)

\(c\)

diffusion coefficient

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

\(\beta\)

convection coefficient

\(\begin{pmatrix} \frac{2}{\mu} \\ 0 \end{pmatrix}\)

\(f\)

source term

\(- \sigma \frac{U}{2\pi} \, r\)

On Air :

Coefficient

Description

Expression

\(c\)

diffusion coefficient

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

For heat equation (Weak Heat Axis), on Conductor (the temperature isn’t computed on Air) :

Coefficient

Description

Expression

\(d\)

damping or mass coefficient

\(\rho \, C_p\)

\(c\)

diffusion coefficient

\(k\)

\(f\)

source term

\(\sigma \left( \frac{U}{2\pi} \right)^2 \frac{1}{r}\)

On JSON file, the coefficients are writed :

CFPDEs coefficients on JSON file

    "Materials":
    {
        "Conductor":
        {
            "sigma":58e+6,  // S.m-1
            "mu":"4*pi*1e-7",   // kg.m/A2/S2
	        "magnetic_c":"x/mu:x:mu",
            "magnetic_beta":"{2/mu,0}:mu",
            "magnetic_f":"-sigma*U/2/pi*x:sigma:U:x",
            "magnetic_d":"sigma*x:x:sigma",

            "k":380,        // W/m/K
            "rho":10000,    // kg/m3
            "Cp":380,       // J/K/kg
            "heat_c":"k*x:k:x",
            "heat_f":"sigma*(U/(2*pi)+magnetic_dphi_dt)*(U/(2*pi)+magnetic_dphi_dt)/x:sigma:U:magnetic_dphi_dt:x",
            "heat_d":"rho*Cp*x:rho:Cp:x"
        },
        "Air":
        {
            "physics":"magnetic",
            "mu":"4*pi*1e-7",   // kg.m/A2/s2
            "magnetic_c":"x/mu:x:mu",
            "magnetic_beta":"{2/mu,0}:mu"
        }
    }

10. Numeric Parameters

This section show the parameters used to compute the simulation.

Size of mesh : \(1 \, mm\)
Time Parameters :
- Time step : I use adaptative time step. Near of brutal changement of current, the time step is little, in other case, the time step is big. I put adaptative time step to have good precision near of changement (cut of current at \(t=22s\)) and to not explose the time of compute. The shell script feelpp_adaptative.sh implements that.
  
  \(\Delta_t = \left\{ \begin{matrix} 0.1s \, \text{for} \, 0s<t<0.9s \\ 0.008s \, \text{for} \, 0.9s<t<1.1s \\ 0.1s \, \text{for} \, 1.1s<t<19.9s \\ 0.008s \, \text{for} \, 19.9s<t<20.1s \\ 0.1 \, \text{for} \, 20.1s<t<22s \end{matrix} \right.\)
- Initial Time : \(0 \, s\)
- Final Time : \(22 \, s\)
Element type : \(P2\) for magnetic equation (\(\mathbf{A}\) and \(P1\) heat equation (\(T\))

Mesh of Geometry

11. Exact solutions

11.1. Intensity

In this subsection, we see the analytical compute of intensity.

We have :

\[ U = R \, I + L \, \frac{dI}{dt} \hspace{2cm} \text{(Intensity)}\]

In stationary case we have \(U = R \, I\).

And :

\[ U = \left\{ \begin{matrix} t & \text{ for } 0s \leq t \leq 1s \\ 1 & \text{ for } 1s \leq t \leq 20s \\ 0 & \text{ for } 20 \leq t \\ \end{matrix} \right.\]

Resolution of first part (\(0s \leq t \leq 1s\)) :

\[ R \, I + L \, \frac{dI}{dt} = t \hspace{2cm} \text{(Intensity-1)}\]

The homogenous equation \(R \, I_0 + L \, \frac{dI_0}{dt} = 0\) have a solve : \(I_0(t) = e^{-\frac{R}{L}t}\)

We pose \(C(t)\), such that \(I(t) = C(t) \, I_0(t) = C(t) \, e^{-\frac{R}{L}t}\) is solution of equation (Intensity-1).

\[\begin{eqnarray*} R \, C(t) \, e^{-\frac{R}{L}t} + L \, \left( \frac{dC}{dt}(t) e^{-\frac{R}{L}t} - \frac{R}{L} \, C(t) \, e^{-\frac{R}{L}t} \right) &=& t \\ 0 + L \frac{dC}{dt}(t) e^{-\frac{R}{L}t} &=& t \\ \frac{dC}{dt}(t) &=& \frac{1}{L} \, t \, e^{\frac{R}{L}t} \\ \end{eqnarray*}\]

Thus :

\[\begin{eqnarray*} C(t) &=& \int_{\tau = 0}^t{ \frac{1}{L} \, \tau \, e^{\frac{R}{L} \tau}} \\ &=& \left[ \frac{1}{L} \, \tau \, \frac{L}{R} \, e^{\frac{R}{L} \tau} \right]_{\tau=0}^t - \int_{\tau=0}^t{ \frac{1}{L} \frac{L}{R} e^{\frac{R}{L} \tau} } \text{ by integration by part} \\ &=& \frac{1}{R} \, t \, e^{\frac{R}{L} t} - \left[ \frac{1}{L} \, \left( \frac{L}{R} \right)^2 \, e^{\frac{R}{L}\tau} \right]_{\tau=0}^t \\ &=& \frac{1}{R} \, t \, e^{\frac{R}{L} t} - \frac{1}{L} \, \left( \frac{L}{R} \right)^2 \, \left(e^{\frac{R}{L} t} - 1 \right) \end{eqnarray*}\]

Thus :

\[\begin{eqnarray*} I(t) &=& C(t) \, e^{\frac{R}{L} t} \\ &=& \frac{1}{R} \, t - \frac{1}{L} \, \left( \frac{L}{R} \right)^2 \, \left( 1 - e^{-\frac{R}{L} t} \right) \\ &=& \frac{1}{L} \left( \frac{L}{R} \, t - \left( \frac{L}{R} \right)^2 \, \left( 1 - e^{-\frac{R}{L} t} \right) \right) \end{eqnarray*}\]

Resolution of second part (\(1s \leq t \leq 20s\)) :

\[ R \, I + L \, \frac{dI}{dt} = 1 \hspace{2cm} \text{(Intensity-2)}\]

An evident solution is :

\[ I(t) = C_1 \, e^{-\frac{R}{L} t} + 1 \text{ with } C_1 \text{ a constant}\]

Resolution of third part (\(20s \leq t\)) :

\[ R \, I + L \, \frac{dI}{dt} = 0 \hspace{2cm} \text{(Intensity-3)}\]

An evident solution is :

\[ I(t) = C_2 \, e^{-\frac{R}{L} t} \text{ with } C_2 \text{ a constant}\]

To conclude, we have :

\[ I(t) = \left\{ \begin{matrix} \frac{1}{L} \, \left( \frac{L}{R} t - \left( \frac{L}{R} \right)^2 \left( 1 -e^{-\frac{R}{L}t} \right) \right) & \text{ for } 0s \leq t \leq 1s \\ C_1 \, e^{-\frac{R}{L} t} + 1 & \text{ for } 1s \leq t \leq 20s \\ C_2 \, e^{-\frac{R}{L} t} & \text{ for } 20s \leq t \end{matrix} \right.\]

This solution is compute with python script.

12. Result

The analytical solve of potential magnetic field and magnetic field (the cyan curve on plots) is computed by python3’s module MagnetTools.MagnetTools. The module based on paper : spire.

The analytical solve of intensity (and magnetic field) is computed on python3’s script Edit the file.

12.1. Intensity

Feelpp compute in post-process, the intensity of circuit :

\[\begin{eqnarray*} I(t) &=& \int_{\Omega_c \cap OrOz}{ \mathbf{J} \cdot \mathbf{n} \ ds } \\ &=& \int_{\Omega_c^{axis}}{ \sigma \left( \frac{1}{r} \frac{U}{2\pi} + \frac{1}{r} \frac{\partial \Phi}{\partial t} \right) \, drdz } \end{eqnarray*}\]

The intensity verify the equation :

\[ U(t) = R \, I(t) + L \, \frac{d I}{d t} \hspace{1cm} \text{(intensity)}\]

With \(R\) the resistance of torus of conductor and \(L\) the inductance of conductor.

In case of a rectangular cross section \((r_{int},r_{ext}) \times (z_1,z_2)\) torus, we can show that :

\[ R = \frac{2\pi \, r_{int} \, \rho}{r_{int} \, ln(\frac{r_{ext}}{r_{int}}) \, (z_1 - z_2)} = \frac{2\pi \, \rho}{ln(\frac{r_{ext}}{r_{int}}) \, (z_1 - z_2)}\]

We have the value :

Symbol

Description

Value

Unit

resistance of torus

\(7.5313 \times 10^{-6}\)

Ohm

inductance of torus

\(1.9204 \times 10^{-7}\)

Henry

The exact solution is (see the subsection Intensity for more detail) :

Intensity of current \(I(A)\)

In stationary case, we have the Ohm law :

\[ U(t \approx \infty) = R \, I(t \approx \infty)\]

With this relation, we have theorycal intensity : \(I_{th} = \frac{U}{R} = -132779 \, A\). And with numeric solve, we have \(I_{num} = -1.3376228966881469e+05 \, A\), so an error of \(7.4e-3 \, A\).

12.2. Magnetic field \(B_z\)

Recall, with axisymmetric condition :

\[ \mathbf{B} = \begin{pmatrix} B_r \\ 0 \\ B_z \end{pmatrix}_{cyl} = \nabla \times \mathbf{A} = \begin{pmatrix} - \frac{\partial A_{\theta}}{\partial z} \\ 0 \\ \frac{1}{r} \frac{\partial \left( r \, A_{\theta} \right)}{\partial z} \end{pmatrix}\]

Thus \(B_z = \frac{1}{r} \frac{\partial \Phi}{\partial r} \).