Transient Heat Transfer (Mixed Form)

Summary

Solves a transient heat conduction problem using a mixed weak form, representing temperature on piecewise constant $var element = document.getElementById("moose-equation-942e12ad-1ce6-4628-bdd4-90f0a992971b");katex.render("L^2", element, {displayMode:false,throwOnError:false});$ conforming finite elements, and heat fluxes on $var element = document.getElementById("moose-equation-5e3f57c3-5920-4c86-b760-d39e40aae989");katex.render("H(\\mathrm{div})", element, {displayMode:false,throwOnError:false});$ conforming Raviart-Thomas elements.

Description

This problem solves the transient heat equation with strong form:

$var element = document.getElementById("moose-equation-a032b72a-676b-419a-b94c-d969580048e8");katex.render("\\begin{split} \\rho c_p \\frac{\\partial T}{\\partial t} - \\vec \\nabla \\cdot \\left(k\\vec\\nabla T\\right) = 0 \\,\\,\\,&\\mathrm{on}\\,\\, \\Omega \\\\ k \\vec \\nabla T \\cdot \\hat n = 0 \\,\\,\\, &\\mathrm{on}\\,\\, \\Gamma_\\mathrm{h} \\\\ T = T_0 \\,\\,\\, &\\mathrm{on}\\,\\, \\Gamma_\\mathrm{T} \\end{split}", element, {displayMode:true,throwOnError:false});$

where $var element = document.getElementById("moose-equation-d2dbd007-8c16-4315-b13f-4ee86791ec43");katex.render("\\rho c_p = 1.0", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-a0aff31b-bbaa-46e6-8c8b-a84b6039185e");katex.render("k^{-1} = 1.0", element, {displayMode:false,throwOnError:false});$ , subject to the initial condition $var element = document.getElementById("moose-equation-52fc2ee7-46ff-444a-b84d-415f26f7820f");katex.render("T(t=0)=0.0", element, {displayMode:false,throwOnError:false});$ .

In this example, we solve this by introducing the heat flux $var element = document.getElementById("moose-equation-1d322103-3544-4563-a3dd-60c086c65907");katex.render("\\vec h = -k \\vec \\nabla T", element, {displayMode:false,throwOnError:false});$ to generate the mixed weak form of the transient heat equation

$var element = document.getElementById("moose-equation-395842ff-a62f-4d15-952c-42768dffa45a");katex.render("\\begin{split} \\left(\\rho c_p \\frac{\\partial T}{\\partial t}, T'\\right)_\\Omega + \\left(\\vec \\nabla \\cdot \\vec h, T'\\right)_\\Omega = 0 \\,\\,\\, \\forall T' \\in V_T \\\\ \\left(k^{-1} \\vec h, \\vec h'\\right)_\\Omega - (T, \\vec \\nabla \\cdot \\vec h')_\\Omega = (T , \\vec h' \\cdot \\hat n)_{\\partial \\Omega} \\,\\,\\, \\forall \\vec h' \\in V_h \\end{split}", element, {displayMode:true,throwOnError:false});$

where

$var element = document.getElementById("moose-equation-d45364fa-b7f8-4906-aedf-e02c1348195e");katex.render("\\begin{split} T, T' &\\in L_2(\\Omega)\\\\ \\vec h, \\vec h' &\\in H(\\mathrm{div})(\\Omega) : \\vec h \\cdot \\hat n = 0 \\,\\,\\, \\mathrm{on}\\,\\, \\Gamma_\\mathrm{h} \\end{split}", element, {displayMode:true,throwOnError:false});$

and the polynomial order of the FEs used to represent $var element = document.getElementById("moose-equation-6e9db4be-8712-454d-9783-6ecaea67ddda");katex.render("\\vec h", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-e61412ae-ab7c-4657-aa66-963399ad05a2");katex.render("\\vec h'", element, {displayMode:false,throwOnError:false});$ is one order higher than $var element = document.getElementById("moose-equation-68437a84-e1a1-424a-9524-3156cf8133cc");katex.render("T", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-d352d677-b36e-4924-961b-533f6edb3bd6");katex.render("T'", element, {displayMode:false,throwOnError:false});$ .

Notably, in contrast to the primal weak form for the transient heat equation solved here, the heat flux on $var element = document.getElementById("moose-equation-2703d4a8-2b0b-4cdb-9b21-278017ccfbd5");katex.render("\\Gamma_\\mathrm{h}", element, {displayMode:false,throwOnError:false});$ is strongly imposed via a Dirichlet condition, and temperatures on boundaries $var element = document.getElementById("moose-equation-cf3640b4-2a5f-4e02-986d-b3f418783ff8");katex.render("\\Gamma_\\mathrm{T}", element, {displayMode:false,throwOnError:false});$ are imposed weakly via boundary integrals.

Example File

[Mesh<<<{"href": "../Mesh/index.html"}>>>]
  type = MFEMMesh
  file = ../mesh/square.e
[]

[Problem<<<{"href": "../Problem/index.html"}>>>]
  type = MFEMProblem
[]

[FESpaces]
  [HDivFESpace]
    type = MFEMVectorFESpace
    fec_type = RT
    fec_order = FIRST
  []
  [L2FESpace]
    type = MFEMScalarFESpace
    fec_type = L2
    fec_order = CONSTANT
  []
[]

[Variables<<<{"href": "../Variables/index.html"}>>>]
  [time_integrated_heat_flux]
    type = MFEMVariable
    fespace = HDivFESpace
    time_derivative = heat_flux
  []
  [temperature]
    type = MFEMVariable
    fespace = L2FESpace
  []
[]

[Kernels<<<{"href": "../Kernels/index.html"}>>>]
  [dT_dt,T']
    type = MFEMTimeDerivativeMassKernel
    variable = temperature
  []
  [divh,T']
    type = MFEMVectorFEDivergenceKernel
    trial_variable = heat_flux
    variable = temperature
  []
  [h,h']
    type = MFEMTimeDerivativeVectorFEMassKernel
    variable = time_integrated_heat_flux
  []
  [-T,div.h']
    type = MFEMVectorFEDivergenceKernel
    trial_variable = temperature
    variable = time_integrated_heat_flux
    coefficient = -1.0
    transpose = true
  []
[]

[BCs<<<{"href": "../BCs/index.html"}>>>]
  [right]
    type = MFEMVectorFEBoundaryFluxIntegratedBC
    variable = time_integrated_heat_flux
    coefficient = 0.0
    boundary = 2
  []
  [left]
    type = MFEMVectorFEBoundaryFluxIntegratedBC
    variable = time_integrated_heat_flux
    coefficient = -1.0
    boundary = 4
  []
  [topbottom]
    type = MFEMVectorDirichletBC
    variable = time_integrated_heat_flux
    vector_coefficient = '0.0 0.0'
    boundary = '1 3'
  []
[]

[Solver]
  type = MFEMSuperLU
[]

[Executioner<<<{"href": "../Executioner/index.html"}>>>]
  type = MFEMTransient
  device = cpu
  assembly_level = legacy
  dt = 0.03
  start_time = 0.0
  end_time = 0.09
[]

[Outputs<<<{"href": "../Outputs/index.html"}>>>]
  [ParaViewDataCollection]
    type = MFEMParaViewDataCollection
    file_base = OutputData/MixedHeatTransfer
    vtk_format = ASCII
  []
[]

(moose/test/tests/mfem/kernels/mixed_heattransfer.i)

[Mesh]
  type = MFEMMesh
  file = ../mesh/square.e
[]

[Problem]
  type = MFEMProblem
[]

[FESpaces]
  [HDivFESpace]
    type = MFEMVectorFESpace
    fec_type = RT
    fec_order = FIRST
  []
  [L2FESpace]
    type = MFEMScalarFESpace
    fec_type = L2
    fec_order = CONSTANT
  []
[]

[Variables]
  [time_integrated_heat_flux]
    type = MFEMVariable
    fespace = HDivFESpace
    time_derivative = heat_flux
  []
  [temperature]
    type = MFEMVariable
    fespace = L2FESpace
  []
[]

[Kernels]
  [dT_dt,T']
    type = MFEMTimeDerivativeMassKernel
    variable = temperature
  []
  [divh,T']
    type = MFEMVectorFEDivergenceKernel
    trial_variable = heat_flux
    variable = temperature
  []
  [h,h']
    type = MFEMTimeDerivativeVectorFEMassKernel
    variable = time_integrated_heat_flux
  []
  [-T,div.h']
    type = MFEMVectorFEDivergenceKernel
    trial_variable = temperature
    variable = time_integrated_heat_flux
    coefficient = -1.0
    transpose = true
  []
[]

[BCs]
  [right]
    type = MFEMVectorFEBoundaryFluxIntegratedBC
    variable = time_integrated_heat_flux
    coefficient = 0.0
    boundary = 2
  []
  [left]
    type = MFEMVectorFEBoundaryFluxIntegratedBC
    variable = time_integrated_heat_flux
    coefficient = -1.0
    boundary = 4
  []
  [topbottom]
    type = MFEMVectorDirichletBC
    variable = time_integrated_heat_flux
    vector_coefficient = '0.0 0.0'
    boundary = '1 3'
  []
[]

[Solver]
  type = MFEMSuperLU
[]

[Executioner]
  type = MFEMTransient
  device = cpu
  assembly_level = legacy
  dt = 0.03
  start_time = 0.0
  end_time = 0.09
[]

[Outputs]
  [ParaViewDataCollection]
    type = MFEMParaViewDataCollection
    file_base = OutputData/MixedHeatTransfer
    vtk_format = ASCII
  []
[]

Summary
Description
Example File