Add time-averaged joule heating auxkernel

README.md

@@ -1,6 +1,6 @@
 # HIVE: Heating by Induction to Verify Extremes
 
-_Nuno Nobre and Karthikeyan Chockalingam_
+_Nuno Nobre, Josh Williams and Karthikeyan Chockalingam_
 
 ###
 
@@ -30,12 +30,14 @@ Installation is as usual for any MOOSE app:
 4) Run it with a no. of `processes` of your choice,
    `mpirun [-n processes] ./hive-opt -w -i input/THeat.i`.
 
-## App design
+## A brief conceptual introduction
 
-For simulation purposes, we discretize HIVE using a tetrahedral mesh, see
-[mesh/](mesh/) for both coarse and fine options, and segment it into three
-components: a cuboid vacuum chamber, an electromagnetic coil, and a target
-prototype component.
+For simulation purposes, we discretize HIVE using a tetrahedral mesh and
+segment it into three components: a cuboid vacuum chamber, an electromagnetic
+coil, and a target prototype component.
+See [mesh/](mesh/) for both a coarse and a fine option, but note the latter is
+_not_ a uniform refinement of the former and that, in particular, the coil is
+discretized using the same number of elements in both.
 Both the coil and the target sit within the vacuum chamber and are assumed to
 be made of electrically conductive materials.
 Here we use copper for both, but other conductors can also be used.
@@ -47,37 +49,46 @@ a time-varying electric current flows through the coil, creating a time-varying
 magnetic field which induces a time-varying electric current in the target
 that gradually warms up via Joule heating.
 
-More specifically, this app leverages a linear, but time-dependent,
-finite element formulation split into three sub-apps. Each
-sub-app solves a different equation for a different field, and feeds
-the solution to the next sub-app. This one-way coupling pipeline
-proceeds as follows:
+Recall that the time-varying electric current flowing through the target will
+self-induce a current, causing the (net) current to flow mostly along the outer perimeter, or skin, of the target. This is known as the skin effect and is more
+pronounced the higher the frequency. Furthermore, since induced currents lag
+their sources by one-quarter cycle, i.e. $\pi/2$, the resulting current on the
+target will, with increasing frequency, tend to lag the current on the coil by
+one-half cycle, i.e. $\pi$.
+
+## App design
+
+This app leverages a linear, but time-dependent, finite element formulation
+split into three sub-apps. Each sub-app solves a different equation for a
+different field, and feeds the solution to the next sub-app. This one-way
+coupling pipeline proceeds as follows:
 
 1) Laplace's equation: $∇^2 V = 0$.
    See [input/VLaplace.i](https://github.com/farscape-project/HIVE/blob/main/input/VLaplace.i).
 
     Solved for the electric potential $V \in \mathcal{P}^1$
     [<sup>(*)</sup>](https://defelement.com/elements/examples/tetrahedron-lagrange-equispaced-1.html),
     only on the coil and only once, with Dirichlet boundary conditions on both
-    its terminals, $V_\mathrm{in} = 1\mathrm{V}$ and
-    $V_\mathrm{out} = 0\mathrm{V}$, and Neumann boundary conditions on its
+    its terminals, $V_\mathrm{in} = V_\mathrm{max}$ and
+    $V_\mathrm{out} = 0$, and Neumann boundary conditions on its
     $\mathbf{n}$-oriented lateral surface, $\mathbf{∇}V \cdot \mathbf{n} = 0$.
     This sole solution can then be scaled appropriately for any time step if
     the time-dependent Dirichlet boundary condition is assumed uniform, i.e.
     $V_\mathrm{in}(\mathbf{r},t) \equiv V_\mathrm{in}(t)$. Here, we take
-    $V_\mathrm{in}(t)=\mathrm{sin}(\omega t)\mathrm{V}$ for some angular
+    $V_\mathrm{in}(t)=V_\mathrm{max}\mathrm{sin}(\omega t)$ for some angular
     frequency $\omega$.
 
 2) The $\mathbf{A}$ formulation: $\mathbf{∇}× \left(ν \mathbf{∇}× \mathbf{A}\right) +σ \partial_t \mathbf{A} = -σ \mathbf{∇}V$.
    See [input/AForm.i](https://github.com/farscape-project/HIVE/blob/main/input/AForm.i).
 
     Solved for the magnetic vector potential $\mathbf{A} \in \mathcal{N}^0_I$
     [<sup>(*)</sup>](https://defelement.com/elements/examples/tetrahedron-nedelec1-lagrange-1.html),
-    everywhere in space and for each time step, with Dirichlet boundary
-    conditions on the $\mathbf{n}$-oriented plane boundary of the vacuum
-    chamber where the coil terminals sit, $\mathbf{A} × \mathbf{n} = 0$, and
-    Neumann boundary conditions on its remaining $\mathbf{n}$-oriented outer
-    surfaces, $\mathbf{∇} × \mathbf{A} × \mathbf{n} = 0$.
+    everywhere in space and for each time step $\Delta t_\mathrm{EM}$ of only a
+    reduced selection of cycles of the voltage source in (1), with Dirichlet
+    boundary conditions on the $\mathbf{n}$-oriented plane boundary of the
+    vacuum chamber where the coil terminals sit, $\mathbf{A} × \mathbf{n} = 0$,
+    and Neumann boundary conditions on its remaining $\mathbf{n}$-oriented
+    outer surfaces, $\mathbf{∇} × \mathbf{A} × \mathbf{n} = 0$.
     $ν$ is the magnetic reluctivity (the reciprocal of the magnetic
     permeability) and $σ$ is the electrical conductivity.
     The right-hand side is non-zero only within the coil, see (1), and is
@@ -88,21 +99,25 @@ proceeds as follows:
 
    Solved for the temperature $T \in \mathcal{P}^1$
    [<sup>(*)</sup>](https://defelement.com/elements/examples/tetrahedron-lagrange-equispaced-1.html),
-   everywhere in space and for each time step, with Neumann boundary conditions
-   on the $\mathbf{n}$-oriented outer surface of the vacuum chamber,
-   $\mathbf{∇}T \cdot \mathbf{n} = 0$, and initial conditions everywhere in
-   space, $T = T_\mathrm{room}$.
+   everywhere in space and for each time step $\Delta t_\mathrm{HT}$, with
+   Neumann boundary conditions on the $\mathbf{n}$-oriented outer surface of
+   the vacuum chamber, $\mathbf{∇}T \cdot \mathbf{n} = 0$, and initial
+   conditions everywhere in space, $T = T_\mathrm{room}$.
    $ρ$ is the density, $c$ is the specific heat capacity, and $k$ is the
    thermal conductivity.
-   The right-hand side is the Joule heating term which, as of this writing, we
-   compute only on the target.
+   The right-hand side is the Joule heating term which we compute only on the
+   target and is time-averaged over the simulated time interval in (2). This
+   enables quicker simulations by solving for the temperature $T$ on a larger
+   time scale than the magnetic vector potential $\mathbf{A}$, i.e.
+   $\Delta t_\mathrm{HT} >> \Delta t_\mathrm{EM}$.
 
 See [input/Parameters.i](input/Parameters.i) for the set of parameters
 influencing the simulation.
-This file is included at the top of both
-[input/AForm.i](input/AForm.i) and [input/THeat.i](input/THeat.i).
-Since neither uses the entire parameter set, MOOSE issues a few warnings that
-can be safely ignored.
+This file is included at the top of the input files for each of the three
+sub-apps, [input/VLaplace.i](input/VLaplace.i), [input/AForm.i](input/AForm.i)
+and [input/THeat.i](input/THeat.i).
+Since none uses the entire parameter set, MOOSE issues a few warnings that
+can be safely ignored (thus the need for the `-w` flag above).
 All material properties are assumed uniform within each of the three components
 and, as of this writing, also temperature-independent.
 
@@ -124,11 +139,6 @@ accuracy, time-to-solution and general usability.
 
 ### Time-to-solution
 
-* Switch to a time-averaged Joule heating source in the heat equation sub-app
-  to keep the problem tractable when simulating over a long physical time span.
-  The $\mathbf{A}$ formulation sub-app will then sub-cycle, i.e. perform
-  multiple time steps for each time step of the heat equation sub-app.
-
 * Study the potential gains of solving all, but most importantly the
   $\mathbf{A}$ formulation sub-app, on the GPU simply via PETSc/hypre flags.
 

include/auxkernels/JouleHeatingAux.h

@@ -16,5 +16,11 @@ class JouleHeatingAux : public AuxKernel
   const VectorVariableValue & _electric_field;
 
   /// The electrical conductivity
-  Real _sigma;
+  const Real _sigma;
+
+  /// Time interval after which the kernel starts computing
+  const Real _skip;
+
+  /// Whether to take the time average
+  const bool _avg;
 };

input/AForm.i

@@ -61,6 +61,7 @@
     variable = P
     vector_potential = A
     sigma = ${copper_econductivity}
+    skip = ${fparse end_t_em/2}
     block = target
     execute_on = timestep_end
   []
@@ -80,7 +81,9 @@
   solve_type = LINEAR
   petsc_options_iname = -pc_type
   petsc_options_value = lu
-  num_steps = 1
+  start_time = 0.0
+  end_time = ${end_t_em}
+  dt = ${delta_t_em}
 []
 
 [MultiApps]

input/Parameters.i

@@ -14,10 +14,12 @@ copper_capacity      = 385                                    # J/(kg*K)
 
 room_temperature     = 293.15                                 # K
 
-voltage_amplitude    = 1                                      # V
+voltage_amplitude    = 0.3052310653                           # V
 voltage_frequency    = 1e5                                    # Hz
 voltage_wfrequency   = ${fparse 2*pi*voltage_frequency}       # rad/s
 voltage_period       = ${fparse 1/voltage_frequency}          # s
 
-end_t                = ${fparse voltage_period}               # s
-delta_t              = ${fparse voltage_period/10}            # s
+end_t_em             = ${fparse voltage_period*10}            # s
+end_t_ht             = 60                                     # s
+delta_t_em           = ${fparse voltage_period/10}            # s
+delta_t_ht           = 5                                      # s

input/THeat.i

@@ -85,8 +85,8 @@
   petsc_options_iname = -pc_type
   petsc_options_value = hypre
   start_time = 0.0
-  end_time = ${end_t}
-  dt = ${delta_t}
+  end_time = ${end_t_ht}
+  dt = ${delta_t_ht}
 []
 
 [Outputs]
@@ -96,9 +96,9 @@
 
 [MultiApps]
   [AForm]
-    type = TransientMultiApp
+    type = FullSolveMultiApp
     input_files = AForm.i
-    execute_on = timestep_begin
+    execute_on = initial
     clone_parent_mesh = true
   []
 []

src/auxkernels/JouleHeatingAux.C

@@ -6,22 +6,31 @@ InputParameters
 JouleHeatingAux::validParams()
 {
   InputParameters params = AuxKernel::validParams();
-  params.addClassDescription(
-      "Computes the differential form of the Joule heating equation (power per unit volume).");
+  params.addClassDescription("Computes (optionally, the time average of) the differential form of "
+                             "the Joule heating equation (power per unit volume). "
+                             "The user may specify a time interval only after which the kernel "
+                             "starts computing. If computing the time average, the right endpoint "
+                             "rectangle rule is used for integration.");
   params.addCoupledVar("vector_potential", "The vector potential variable");
   params.addParam<Real>("sigma", 1, "The electrical conductivity");
+  params.addParam<Real>("skip", 0, "Time interval after which the kernel starts computing");
+  params.addParam<bool>("average", true, "Whether to take the time average");
   return params;
 }
 
 JouleHeatingAux::JouleHeatingAux(const InputParameters & parameters)
   : AuxKernel(parameters),
     _electric_field(coupledVectorDot("vector_potential")),
-    _sigma(getParam<Real>("sigma"))
+    _sigma(getParam<Real>("sigma")),
+    _skip(getParam<Real>("skip")),
+    _avg(getParam<bool>("average"))
 {
 }
 
 Real
 JouleHeatingAux::computeValue()
 {
-  return _sigma * _electric_field[_qp] * _electric_field[_qp];
+  Real p = _sigma * _electric_field[_qp] * _electric_field[_qp];
+  Real w = _t > _skip ? _avg ? _dt / (_t - _skip) : 1 : 0;
+  return (1 - w) * _u[_qp] + w * p;
 }