(FEM) 1D time-dependent heat equation convergence problem

Question

I'm simulating a simple 3-node bar with convection BCs at the edges to validate my FEM code. The following data was used:

Initial temperature = 25 ºC

Temperature surrounding the rod = 10 ºC

Thermal conductivity = 0

Specific heat capacity = Density = Convection coefficient = 1

For time -> infinity the two nodes at the edge converged to 10 ºC (which was expected). However, at the middle node the temperature should had maintained its initial value of 25 ºC since there's no thermal conductivity but it raised up to a higher temperature (32.5 ºC) and I have no clue why is this happening.

Edit1: I followed these steps:

1) Assuming the standard transient heat-equation:

$$\rho.c_p\ \frac{\partial T}{\partial t} - k\frac{\partial^2 T}{\partial x^2} = 0$$

2) After the formal FEM procedure and using Euler-Backward time-integration I arrived to the following formulation:

$$ \rho.c_p \int N_i\ N_j\ d\Omega\ \frac{T_j^{t+1}-T_j^t}{\Delta t} +\ k\int N'_i\ N'_j\ d\Omega\ \ T_j^{t+1}\ + h\int N_i\ N_j\ d\Gamma\ \ T_j^{t+1}\\= h\ T_{surrounding} \int N_i\ d\Gamma$$

$\rho:\ $Density

$c_p:\ $Specific heat capacity

$k:\ $ Thermal conductivity

$h:\ $Convective heat transfer coefficient

If $k\rightarrow \ 0$ then the above related problem appears.

Edit 2: The step-by-step procedure I followed to arrive at the equation presented at point 2 of Edit 1 was the following:

Step 1) Multiply the PDE at point 1 by the weight function (w) and integrate over the domain:

$$\rho c_p \int \frac{\partial T}{\partial t} w\ d\Omega + \int \overrightarrow{\nabla}.(-k\overrightarrow{\nabla} T) w\ d\Omega = 0$$

Step 2) Apply divergence theorem on the second term at the Left Hand-Side of the equation:

$$\rho c_p \int \frac{\partial T}{\partial t} w\ d\Omega - k\oint (\overrightarrow{\nabla} T .\overrightarrow{n}) w\ d\Gamma + k\int \overrightarrow{\nabla} T . \overrightarrow{\nabla} w\ d\Omega = 0$$

or

$$\rho c_p \int \frac{\partial T}{\partial t} w\ d\Omega + k\int \overrightarrow{\nabla} T . \overrightarrow{\nabla} w\ d\Omega = \oint k(\overrightarrow{\nabla} T .\overrightarrow{n}) w\ d\Gamma$$

where $\overrightarrow{n}$ and $\overrightarrow{\nabla}$ denotes the outward normal vector at the boundary and the gradient vector, respectively.

Step 3) Using the following relation

$$k\overrightarrow{\nabla}T.\overrightarrow{n}\ + h(T-T_{\infty}) = 0$$

at the Right Hand-Side of latter equation at Step 2 one arrives to:

$$\rho c_p \int \frac{\partial T}{\partial t} w\ d\Omega + k\int \overrightarrow{\nabla} T . \overrightarrow{\nabla} w\ d\Omega = \oint -h(T-T_{\infty}) w\ d\Gamma$$

Where $$T_{\infty} = T_{surrounding}$$ is the surrounding temperature.

Step 4) Rearranging the Right Hand-Side

$$\rho c_p \int \frac{\partial T}{\partial t} w\ d\Omega + k\int \overrightarrow{\nabla} T . \overrightarrow{\nabla} w\ d\Omega + \oint h\ T\ w\ d\Gamma = \oint h\ T_{\infty}\ w\ d\Gamma$$

and using $T = N_jT_j$, $w = N_i$ and assuming a Backward-Euler time integration scheme we'll have

$$\rho c_p \int N_iN_j\ d\Omega\ \frac{T_j^{t+1}-T_j^t}{\Delta t} + k\int \overrightarrow{\nabla}N_i.\overrightarrow{\nabla} N_j\ d\Omega\ T_j^{t+1} + h\oint N_i N_j\ d\Gamma\ T_j^{t+1} = \\h\oint T_{\infty} N_i\ d\Gamma$$

Or in matrix form

$$[M]\frac{T^{t+1}-T^t}{\Delta t} + [K]T^{t+1} + [H]T^{t+1} = F$$

Answer 1

You don't need MFnx answer to find out why if $k = 0$, $h$ also MUST be zero ($h = 0$). Look at your initial PDE for heat transfer:

$$\rho C_{p}\frac{\partial T}{\partial t} = k \frac{\partial^{2} T}{\partial x^{2}}$$

If $k = 0$, so:

$$\rho C_{p} \frac{\partial T}{\partial t} = 0$$

Since your density and heat capacity are not zero, then:

$$\frac{\partial T}{\partial t} = 0$$

It means: $T = T_{0}$ every where in your 1D domain, where $T_{0}$ is your temperature at initial time or $t = 0$, cause you MUST satisfy initial condition.

On the other hand, you have this convection BC at the boundaries: $-k \frac{\partial T}{\partial x} = h (T - T_{\infty})$. Since $k = 0$, then: $h (T-T_{\infty}) = 0$. Let's assume that $h$ is not zero, then, $T - T_{\infty}$ MUST be zero: $T = T_{\infty}$. It's CONTRADICTION. why?! Cause, there is no reason whatsoever that $T_{\infty}$ MUST be equal to $T_{0}$. Then, my initial assumption of $h$ non-zero is FALSE. So, if $k= 0$, I MUST have $h = 0$. Otherwise, I would see a clear contradiction as you saw in your numerical scheme that your initial point reached a 32 degree Celsius, despite the fact that it MUST remain at 25 degree Celsius. You chose a really bad example to verify your FEM scheme.

Answer 2

From the comments, I believe you are struggling with the FEM formulations and therefore fail to see what happens if the diffusion coefficient tends to zero. I don't have access to your book, but I will post here an answer which shows how to discretize the heat equation with the standard Galerkin method. You will then see that, if there is no diffusion, there is no change of temperature within the material.

The heat equation

The linear heat equation in absence of any heat source can be written as:

$$ \frac{\partial T}{\partial t} - \alpha \nabla \cdot \nabla T = 0 $$

which becomes in 1D:

$$ \frac{\partial T}{\partial t} - \alpha \frac{\partial^2 T}{\partial x^2} = 0 $$

where the thermal diffusivity $\alpha = \frac{k}{\rho c_p}$.

Note that, if $\alpha = 0$ (or $k = 0$), then we would have:

$$ \frac{\partial T}{\partial t} = 0 $$

which is a different differential equation!

Discretizing the time derivative, we obtain:

$$ \frac{T^{i+1} - T^{i}}{\Delta t} = 0 \Rightarrow T^{i+1} = T^{i} \text{ } \forall \text{ } (t,x) \text{ in } \Omega_t \times \Omega $$

where $\Omega$ and $\Omega_t$ refer to the spatial and time spaces.

This should be sufficient to show that, in absence of thermal diffusivity, temperature does not change over time within $\Omega$. However, in order to show that the same holds after discretizing the differential equation by means of the standard FEM, let's develop the strong, weak and matrix formulations of the initial boundary value problem.

Strong formulation

The strong form of the initial boundary value problem can be stated as follows:

$$ \text{Given } T_0 : \Omega \rightarrow \mathbb{R} \text{, } \hat{T} : \Omega_t \rightarrow \mathbb{R} \text{ and } q : \Omega_t \rightarrow \mathbb{R} \text{, find } T : \Omega_t \times \Omega \rightarrow \mathbb{R} \text{ such that:} $$ $$ \frac{\partial T}{\partial t} - \alpha \frac{\partial^2 T}{\partial x^2} = 0 \text{ in } \Omega_t \times \Omega \\ T = T_0(x) \text{ in } \Omega_{t=0} \times \Omega \\ T = \hat{T}(t) \text{ on } \Omega_t \times \Gamma_D \\\ -\alpha \frac{\partial T}{\partial x} = q(t) \text{ on } \Omega_t \times \Gamma_N $$

where $\Gamma_D$ and $\Gamma_N$ are the Dirichlet and Neumann boundaries.

Weak formulation

We will use the following notation:

$$ <a,b> = \int_\Omega ab d\Omega \\ <<a,b>> = \int_{\Gamma} ab d\Gamma \\ <<a,b>>_N = \int_{\Gamma_N} ab d\Gamma \\ \dot{T} = \frac{\partial T}{\partial t} $$

where $\Gamma$ refers to the boundary of $\Omega$ ($\Gamma = \Gamma_D \cup \Gamma_N$).

Let's define $w$ to be the test functions which are 0 on $\Gamma_D$. Also, let's assume that $T$ satisfies the boundary conditions defined on $\Gamma_D$ in the strong formulation. Typically in FEM analysis, the Dirichlet boundary conditions are imposed once the matrix formulation is obtained from the weak formulation.

Then, in order to obtain the weak formulation, we multiply the differential equation by the test function and integrate on $\Omega$:

$$ <w,\dot{T}> + <w,-\alpha \frac{\partial^2 T}{\partial x^2}> = 0 $$

By virtue of the chain rule, this can be rewritten as:

$$ <w,\dot{T}> + <\frac{\partial w}{\partial x},\alpha \frac{\partial T}{\partial x}> + \int_\Omega \frac{\partial }{\partial x} \big(-w \alpha \frac{\partial T}{\partial x} \big)d\Omega = 0 $$

Thanks to the divergence theorem, we can finally write:

$$ <w,\dot{T}> + <\frac{\partial w}{\partial x},\alpha \frac{\partial T}{\partial x}> + <<w,-\alpha \frac{\partial T}{\partial x}>> = 0 $$

Since $w = 0 \text { on } \Gamma_D$ and $-\alpha \frac{\partial T}{\partial x} = q(t)$, the weak formulation of the initial boundary value problem reads:

$$ <w,\dot{T}> + <\frac{\partial w}{\partial x},\alpha \frac{\partial T}{\partial x}> + <<w,q(t)>>_N = 0 $$

Note that, if $k=0$, so should $q(t)$, and once more we arrive to the conclusion that $T$ doesn't change over time if there is no thermal diffusivity. If $k=0$ but $q(t) \neq 0$, than you are not solving the problem stated in the strong formulation.

Matrix formulation

The matrix formulation can now be obtained easily from the weak formulation by using the standard Galerking discretization. Using a similar notation as in the question:

$$ \int_\Omega N_i N_j d\Omega \dot{T}_j + \int_\Omega \alpha\frac{dN_i}{dx}\frac{dN_j}{dx}d\Omega T_j + \int_{\Gamma_N} N_i q(t) d\Gamma = 0 $$

Again, note that, $k=0 \Rightarrow q(t) = 0$.

It seems that in the book you mentioned, they used the implicit Euler method. Applying the same here gives the following result:

$$ \int_\Omega N_i N_j d\Omega \frac{T_j^{i+1}-T_j^i}{\Delta t} + \int_\Omega \alpha\frac{dN_i}{dx}\frac{dN_j}{dx}d\Omega T_j^{i+1} + \int_{\Gamma_N} N_i q(t) d\Gamma = 0 $$

If $k=0$, then the second and third term vanish, and you are left with $T^{i+1} = T^{i} \text{ } \forall \text{ } (t,x) \text{ in } \Omega_t \times \Omega$

EDIT 1

Even if you set $k=0$ and $q(t) \neq 0$, the inner nodal values should not change, as the imposition of the natural boundary conditions should only affect the nodal values of those nodes which ly on the Neumann boundary. Indeed, since $k=0$, it is impossible to transfer heat from a surface node to an inner node.

EDIT 2

To take into account convection at the boundary of the material, you can impose the temperature gradient on the Neumann boundary as follows:

$$ q(t) = -\alpha \nabla T = \frac{h}{\rho c_p}(T-T_{sur}) $$

or, equivalently:

$$ -k \nabla T = h(T-T_{sur}) $$