Problem with Analytical Solution

Simulation Setup

The following instructions assumes that there is a working version of MOOSE installed on the system, for information on getting it setup please refer to the install instructions.

The simulation to be performed relies on the heat transfer module, thus it is necessary to compile that specific module. It is assumed that MOOSE was cloned into a directory named "projects" within the home directory. With this assumption, the following commands will compile and test the heat transfer module executable.

cd ~/projects/moose/modules/heat_conduction
make -j8
./run_tests -j8

The number used after the "-j" should be the number of processes available on your system. If tests execute without failure, then the following sections should be followed to setup the simulation.

Simulation Input

MOOSE operates using input files, as such, an input file is detailed below that is capable of simulating this problem and its use demonstrated. This section will describe the various sections of the input file that should be created for the current problem. The location of this input file is arbitrary, but it is assumed that the content below will be in a file named ~/projects/problems/verification/1d_analytical.i. This file can be created and edited using any text editor program.

Within this file, the domain of the problem is defined using the Mesh System using the [Mesh] block of the input file, as follows, which creates a one-dimensional mesh in the -direction from 0 to 0.03 using 200 elements.

[Mesh<<<{"href": "../../../syntax/Mesh/index.html"}>>>]
  [gen]
    type = GeneratedMeshGenerator<<<{"description": "Create a line, square, or cube mesh with uniformly spaced or biased elements.", "href": "../../../source/meshgenerators/GeneratedMeshGenerator.html"}>>>
    dim<<<{"description": "The dimension of the mesh to be generated"}>>> = 1
    xmax<<<{"description": "Upper X Coordinate of the generated mesh"}>>> = 0.03
    nx<<<{"description": "Number of elements in the X direction"}>>> = 200
  []
[]

There is a single unknown, temperature (), to compute. This unknown is declared using the Variables System in the [Variables] block and used the default configuration of a first-order Lagrange finite element variable. Note, within this block the name of the variable is arbitrary; "T" was selected here to match the equation.

[Variables<<<{"href": "../../../syntax/Variables/index.html"}>>>]
  [T]
  []
[]

The "volumetric" portions of the weak form are defined using the Kernel System in the [Kernels] block, for this example this can be done with the use of two Kernel objects as follows.

[Kernels<<<{"href": "../../../syntax/Kernels/index.html"}>>>]
  [T_time]
    type = HeatConductionTimeDerivative<<<{"description": "Time derivative term $\\rho c_p \\frac{\\partial T}{\\partial t}$ of the thermal energy conservation equation.", "href": "../../../source/kernels/HeatConductionTimeDerivative.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = T
    density_name<<<{"description": "Property name of the density material property"}>>> = 7800
    specific_heat<<<{"description": "Name of the volumetric isobaric specific heat material property"}>>> = 450
  []
  [T_cond]
    type = HeatConduction<<<{"description": "Diffusive heat conduction term $-\\nabla\\cdot(k\\nabla T)$ of the thermal energy conservation equation", "href": "../../../source/kernels/HeatConduction.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = T
    diffusion_coefficient<<<{"description": "Property name of the diffusion coefficient"}>>> = 80.2
  []
[]

The sub-block "T_time" defines the time derivative and "T_cond" defines the conduction portion of the equation, please refer to Heat Conduction for details regarding the weak form of the heat equation. Both blocks include the "variable" parameter, which is set equal to "T". The remaining parameters define the values for , , and as the constant values defined in Table 1.

The boundary portions of the weak form are defined using the Boundary Condition System in the [BCs] block. At a Neumann condition is applied with a value of 7E5. The generated mesh, by default, labels the boundary at as "left".

[BCs<<<{"href": "../../../syntax/BCs/index.html"}>>>]
  [left]
    type = NeumannBC<<<{"description": "Imposes the integrated boundary condition $\\frac{\\partial u}{\\partial n}=h$, where $h$ is a constant, controllable value.", "href": "../../../source/bcs/NeumannBC.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = T
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = left
    value<<<{"description": "For a Laplacian problem, the value of the gradient dotted with the normals on the boundary."}>>> = 7e5
  []
[]

The natural boundary condition does not require definition, as the name suggests, this condition is applied naturally based on the weak form derivation.

The problem is solved using Newton's method with a second-order backward difference formula (BDF2) with a timestep of 0.01 seconds up to a simulation time of one second. These settings are applied within the [Executioner] block using the Executioner System.

[Executioner<<<{"href": "../../../syntax/Executioner/index.html"}>>>]
  type = Transient
  scheme = bdf2
  solve_type = NEWTON
  dt = 0.01
  end_time = 1
[]

Finally, the output method is defined. In this case the ExodusII and CSV format are enabled within the [Outputs] block using the Outputs System.

[Outputs<<<{"href": "../../../syntax/Outputs/index.html"}>>>]
  exodus<<<{"description": "Output the results using the default settings for Exodus output."}>>> = true
  csv<<<{"description": "Output the scalar variable and postprocessors to a *.csv file using the default CSV output."}>>> = true
[]

Simulation Execution

Executing the simulation is straightforward, simply execute the heat transfer module executable with the input file included using the "-i" option as follows.

cd ~/projects/moose/tutorials/tutorial03_verification/step03_analytical
~/projects/moose/modules/heat_transfer/heat_conduction-opt -i 1d_analytical.i

When complete an output file will be produced with the name "1d_analytical_out.e", this file can be viewed using Paraview or Peacock.

Spatial Convergence

To perform a spatial convergence the input file created needs to be executed with decreasing element size, which in MOOSE can be done by adding the Mesh/uniform_refine=x option to the command line with "x" being an integer representing the number of refinements to perform.

The mms package includes a function (run_spatial) for performing and automatically modifying the supplied input file to perform increasing levels of uniform refinement. After each run the results are collected into a pandas.DataFrame object, which can be used in the ConvergencePlot object to create a plot.

For example, for Listing 1 results of the spatial convergence study are shown as a plot in Figure 2. In this case, the input file created above (1d_analytical.i) is executed with six levels of uniform refinement (0–5). Notice that an extra argument Mesh/nx=10 is supplied to reduce the number of elements to 10 for the initial mesh. The console flag suppress the screen output from the simulations.

import mms
df = mms.run_spatial('1d_analytical.i', 6, 'Mesh/gen/nx=10', console=False)
fig = mms.ConvergencePlot(xlabel='Element Size ($h$)', ylabel='$L_2$ Error')
fig.plot(df, label='1st Order', marker='o', markersize=8)
fig.save('1d_analytical_spatial.png')

The resulting convergence data indicate that the simulation is converging to the known solution at the correct rate for first-order elements, thus verifying that the finite element solution is computing the results as expected.

Figure 2: Results of the spatial convergence study comparing the finite element simulation with the known one-dimensional analytical solution.

Temporal Convergence

The temporal convergence study is nearly identical to the spatial, as shown in Listing 2. For a temporal study, the run_temporal is used and the initial time step is set using dt; the number of elements is not altered. The resulting converge graph is shown in Figure 3.

import mms
df = mms.run_temporal('1d_analytical.i', 6, dt=0.5, console=False)
fig = mms.ConvergencePlot(xlabel='$\Delta t$', ylabel='$L_2$ Error')
fig.plot(df, label='2nd Order (Backward Difference)', marker='o', markersize=8)
fig.save('1d_analytical_temporal.png')

The resulting convergence data indicate that the simulation is converging to the known solution at the correct rate for second-order time integration, thus verifying that the numerical integration scheme is operating as expected.

Figure 3: Results of the spatial convergence study comparing the finite element simulation with the known one-dimensional analytical solution.