We consider homogeneous parabolic equation here: $$ \begin{equation} \rho C_\textrm{p}\frac{\partial T}{\partial t} + \nabla\cdot(\lambda \nabla T) = 0 \quad \text{in }\Omega \end{equation}$$ w.r.t boundary conditions $$ \eqalign{ T(x, t=0) = T_0,\cr T(x) = g_D(x) &\quad \text{on }\Gamma_D,\cr \lambda {\partial T(x) \over \partial n} = g_N(x) &\quad \text{on }\Gamma_N, } $$ where $T$ could be temperature, the subscripts $D$ and $N$ denote the Dirichlet and Neumanntype boundary conditions, $n$ is the normal vector pointing outside of $\Omega$, and $\Gamma = \Gamma_D \cup \Gamma_N$ and $\Gamma_D \cap \Gamma_N = \emptyset$.
We solve the parabolic equation on a line domain $[60\times 1]$ with $\lambda = 3.2$ and $\rho C_\textrm{p} = 2.5 \times 10^6$ w.r.t. the specific boundary conditions:
$$ \eqalign{ \lambda {\partial T(x) \over \partial n} = q &\quad \text{on } (x=0) \subset \Gamma_N.\cr \lambda {\partial T(x) \over \partial n} = 0 &\quad \text{on } (x=60) \subset \Gamma_N. } $$
with flux value being $q = 2$.
Approximate solution of this problem is the solution for semiinfinite rod, which is applicable as long as the temperature does not come close to the right side of the domain: $$ \begin{equation} T(x,t) = \frac{2q}{\lambda}\left(\left(\frac{\alpha t}{\pi}\right)^{\frac{1}{2}} e^{x^2/4\alpha t}  \frac{x}{2}\textrm{erfc}\left( \frac{x}{2\sqrt{\alpha t}}\right)\right), \end{equation} $$ where $T_\textrm{b}$ is the boundary temperature, $\textrm{erfc}$ is the complementary error function and $\alpha = \lambda/(C_p \rho)$ is the thermal diffusivity.
The main project file is picard.prj
. It describes the processes to be solved and the related process variables together with their initial and boundary conditions. It also references the mesh and geometrical objects defined on the mesh.
The geometries used to specify the boundary conditions are given in the
line_60_geometry.gml
file. The input mesh is stored in the mesh.vtu
file.
Same project was also setup to be solved with the NewtonRaphson method, see
newton.prj
input file.
Another two project files were tested with a masslumping method used to
reduce oscillations for particularly small ratios of element size to time step
increment. See picard_masslumping.prj
and newton_masslumping.prj
files.
For the simple equation, as expected, the difference between the Picard and Newton nonlinear solvers is reasonable for the given nonlinear solver tolerances.
Using the masslumping method reduces the oscillations for the first time steps on cost of accuracy as the error is significantly larger.
This article was written by Dmitri Naumov, Tianyuan Zheng. If you are missing something or you find an error please let us know.
Generated with Hugo 0.117.0
in CI job 389515

Last revision: November 7, 2023
Commit: [PL/SF] clangformat 2a6e9af
 Edit this page on