Nonlinear finite elements/Solution of heat equation

Finite element solution for the Heat equation

Let us now try to create a finite element approximation for the variational initial boundary value problem for the heat equation . This problem can be stated as

$\begin{array}{cc}& 𝖵𝖺𝗋𝗂𝖺𝗍𝗂𝗈𝗇𝖺𝗅𝖨𝖡𝖵𝖯𝖿𝗈𝗋𝗍𝗁𝖾𝖧𝖾𝖺𝗍𝖤𝗊𝗎𝖺𝗍𝗂𝗈𝗇\\ & \\ \text{Find a function}& T(t)\in {𝒮}_t,t\in [0,\tau ]\text{such that for all}w\in 𝒱\\ & \\ & \underset{\mathrm{\Omega }}{\int }\left(\rho C_v\frac{\partial T}{\partial t}\right)wdV+\underset{\mathrm{\Omega }}{\int }(\mathrm{\nabla }w)\bullet (\mathit{\kappa }\bullet \mathrm{\nabla }T)dV=\underset{\mathrm{\Omega }}{\int }swdV-\underset{{\mathrm{\Gamma }}_q}{\int }w\bar{q}dA\\ & \underset{\mathrm{\Omega }}{\int }w\rho C_{v}T(0)dV=\underset{\mathrm{\Omega }}{\int }w\rho C_v{T}_{0}dV.\end{array}$

Let ${𝒱}_h$ be a finite dimensional approximation to $𝒱$ and let ${𝒮}_h$ be a finite dimensional approximation to $𝒮$. Let the weighting functions $w_h\in {𝒱}_h$ be such that $w_h=0$ on ${\mathrm{\Gamma }}_T$. Similarly, any other function $v_h\in {𝒱}_h$ also goes to zero on the boundary ${\mathrm{\Gamma }}_T$.

Assume that the trial solutions $T_h\in {𝒮}_h$ can be represented as

$T_h=v_h+{\bar{T}}_h\text{where}{v}_h\in {𝒱}_h.$

The additional quantity ${\bar{T}}_h$ results in the satisfaction of the boundary condition $T=\bar{T}$ on ${\mathrm{\Gamma }}_T$.

If we substitute these quantities into the variational initial boundary value problem, we get the Galerkin formulation. The steps in this process are shown below.

Approximate IBVP

The approximate form of the variational IBVP is

$\underset{\mathrm{\Omega }}{\int }\left(\rho C_v\frac{\partial T_h}{\partial t}\right)w_{h}dV+\underset{\mathrm{\Omega }}{\int }(\mathrm{\nabla }{w}_h)\bullet (\mathit{\kappa }\bullet \mathrm{\nabla }{T}_h)dV=\underset{\mathrm{\Omega }}{\int }s{w}_{h}dV-\underset{{\mathrm{\Gamma }}_q}{\int }{w}_h\bar{q}dA.$

Replacing $T_h$ with $(v_h+{\bar{T}}_h)$ we get

$\underset{\mathrm{\Omega }}{\int }\left(\rho C_v\frac{\partial (v_h+{\bar{T}}_h)}{\partial t}\right)w_{h}dV+\underset{\mathrm{\Omega }}{\int }(\mathrm{\nabla }{w}_h)\bullet (\mathit{\kappa }\bullet \mathrm{\nabla }(v_h+{\bar{T}}_h))dV=\underset{\mathrm{\Omega }}{\int }s{w}_{h}dV-\underset{{\mathrm{\Gamma }}_q}{\int }{w}_h\bar{q}dA.$

After expanding and rearranging terms, we get

$\begin{array}{cc}\underset{\mathrm{\Omega }}{\int }{w}_h\rho C_v\frac{\partial v_h}{\partial t}dV+\underset{\mathrm{\Omega }}{\int }(\mathrm{\nabla }{w}_h)\bullet (\mathit{\kappa }\bullet \mathrm{\nabla }{v}_h)dV=& \underset{\mathrm{\Omega }}{\int }{w}_{h}sdV-\underset{{\mathrm{\Gamma }}_q}{\int }{w}_h\bar{q}dA\\ & -\underset{\mathrm{\Omega }}{\int }{w}_h\rho C_v\frac{\partial {\bar{T}}_h}{\partial t}dV-\underset{\mathrm{\Omega }}{\int }(\mathrm{\nabla }{w}_h)\bullet (\mathit{\kappa }\bullet \mathrm{\nabla }{\bar{T}}_h)dV.\end{array}$

In compact notation

$⟨w_h,\rho C_v\dot{v}{⟩}_h+a(w_h,v_h)=⟨w_h,s⟩-⟨w_h,\bar{q}{⟩}_{\mathrm{\Gamma }}-⟨w_h,\rho C_v\dot{\bar{T}}{⟩}_h-a(w_h,{\bar{T}}_h).$

Similarly, the initial condition can be written as

$\underset{\mathrm{\Omega }}{\int }{w}_h\rho C_v{T}_h(0)dV=\underset{\mathrm{\Omega }}{\int }{w}_h\rho C_v{T}_{0}dV.$

Substituting for $T_h$ and expanding out terms, we have

$\underset{\mathrm{\Omega }}{\int }{w}_h\rho C_v{v}_h(0)dV=\underset{\mathrm{\Omega }}{\int }{w}_h\rho C_v{T}_{0}dV-\underset{\mathrm{\Omega }}{\int }{w}_h\rho C_v{\bar{T}}_h(0)dV.$

In compact form,

$⟨w_h,\rho C_v{v}_h(0)⟩=⟨w_h,\rho C_v{T}_0⟩-⟨w_h,\rho C_v{\bar{T}}_h(0)⟩.$

Therefore the approximate form of the variational BVP can be written as

$\text{(40)}\begin{array}{cc}& 𝖠𝗉𝗉𝗋𝗈𝗑𝗂𝗆𝖺𝗍𝖾𝖵𝖺𝗋𝗂𝖺𝗍𝗂𝗈𝗇𝖺𝗅𝖨𝖡𝖵𝖯𝖿𝗈𝗋𝗍𝗁𝖾𝖧𝖾𝖺𝗍𝖤𝗊𝗎𝖺𝗍𝗂𝗈𝗇\\ & \\ \text{Find a function}& T_h(t)=v_h+{\bar{T}}_h\in {𝒮}_h,t\in [0,\tau ]\text{such that for all}{w}_h\in {𝒱}_h\\ & \\ \underset{\mathrm{\Omega }}{\int }{w}_h\rho C_v\frac{\partial v_h}{\partial t}dV& +\underset{\mathrm{\Omega }}{\int }(\mathrm{\nabla }{w}_h)\bullet (\mathit{\kappa }\bullet \mathrm{\nabla }{v}_h)dV=\underset{\mathrm{\Omega }}{\int }{w}_{h}sdV-\underset{{\mathrm{\Gamma }}_q}{\int }{w}_h\bar{q}dA\\ & -\underset{\mathrm{\Omega }}{\int }{w}_h\rho C_v\frac{\partial {\bar{T}}_h}{\partial t}dV-\underset{\mathrm{\Omega }}{\int }(\mathrm{\nabla }{w}_h)\bullet (\mathit{\kappa }\bullet \mathrm{\nabla }{\bar{T}}_h)dV.\\ \underset{\mathrm{\Omega }}{\int }{w}_h\rho C_v{v}_h(0)dV& =\underset{\mathrm{\Omega }}{\int }{w}_h\rho C_v{T}_{0}dV-\underset{\mathrm{\Omega }}{\int }{w}_h\rho C_v{\bar{T}}_h(0)dV.\end{array}$

Note that the derivatives with respect to time are still continuous in this approximate variational BVP.

Finite element approximation

We now discretize our domain into elements ${\mathrm{\Omega }}_e$ where $1<e<E$ and $E$ is the total number of elements. Nodes may exist anywhere within the element domains. But the most common locations are at the element vertices and along the interelement boundaries.

Let the set of global node numbers be

$\eta =\{1,2,\dots ,N\}.$

Let the set nodes at which a temperature (essential BC) is prescribed be

${\eta }^T\subset \eta .$

Then the set of nodes at which a temperature is to be determined is the complement

$\eta -{\eta }^T.$

Let $n$ be the number of nodes in $\eta -{\eta }^T$. Then the number of equations to be solved is also $n$. See Figure 1 to get an idea about what these sets of nodes refer to.

File:HeatCondProbFEM.png Figure 1. FEM discretization for the heat conduction problem.

As we have seen before, a typical weighting function $w_h\in {𝒱}_h$ is assumed to have the form

$\text{(42)}{w}_h(𝐱)=\sum _{i=1}^n{b}_i{N}_i(𝐱)\equiv \sum _{i\in \eta -{\eta }^T}{b}_i{N}_i(𝐱).$

Note that $w_h=0$ only if $b_i=0$ for every $i\in \eta -{\eta }^T$. On the other hand $w_h=0$ for every node in ${\eta }^T$.

Similarly, a typical trial solution $v_h\in {𝒱}_h$ can be written as

$\text{(43)}{v}_h(𝐱,t)=\sum _{i\in \eta -{\eta }^T}{T}_i(t)N_i(𝐱)$

where $T_i(t)$ is the unknown temperature at node $i$ at time $t$.

We also often represent the approximate form of the essential BCs in a similar way, i.e.,

$\text{(44)}{\bar{T}}_h(𝐱,t)=\sum _{i\in {\eta }^T}{\bar{T}}_i(t)N_i(𝐱),{\bar{T}}_i(t)=\bar{T}({𝐱}_i,t).$

Substitute equations (42), (43), and (44) into the first of equations (40). You will get

$\begin{array}{cc}\underset{\mathrm{\Omega }}{\int }\left(\sum _j{b}_j{N}_j\right)& \rho C_v\left(\sum _i\frac{\partial T_i}{\partial t}{N}_i\right)dV+\underset{\mathrm{\Omega }}{\int }\left(\sum _j{b}_j\mathrm{\nabla }{N}_j\right)\bullet \left[\sum _i{T}_i(\mathit{\kappa }\bullet \mathrm{\nabla }{N}_i)\right]dV=\\ & \underset{\mathrm{\Omega }}{\int }\left(\sum _j{b}_j{N}_j\right)sdV-\underset{{\mathrm{\Gamma }}_q}{\int }\left(\sum _j{b}_j{N}_j\right)\bar{q}dA\\ & -\underset{\mathrm{\Omega }}{\int }\left(\sum _j{b}_j{N}_j\right)\rho C_v\left(\sum _k\frac{\partial {\bar{T}}_k}{\partial t}{N}_k\right)dV\\ & -\underset{\mathrm{\Omega }}{\int }\left(\sum _j{b}_j\mathrm{\nabla }{N}_j\right)\bullet \left[\sum _k{\bar{T}}_k(\mathit{\kappa }\bullet \mathrm{\nabla }{N}_k)\right]dV,i,j\in \eta -{\eta }^T,k\in {\eta }^T.\end{array}$

Assuming that $\mathrm{\Omega }$ does not change with time, we can take the constants and time-dependent parameters outside the integrals to get

$\begin{array}{cc}\sum _j{b}_j& [\sum _i\frac{\partial T_i}{\partial t}\left(\underset{\mathrm{\Omega }}{\int }\rho C_v{N}_i{N}_{j}dV\right)+\sum _i{T}_i(\underset{\mathrm{\Omega }}{\int }\mathrm{\nabla }{N}_j\bullet (\mathit{\kappa }\bullet \mathrm{\nabla }{N}_i)dV)]=\\ & \sum _j{b}_j[\underset{\mathrm{\Omega }}{\int }{N}_{j}sdV-\underset{{\mathrm{\Gamma }}_q}{\int }{N}_j\bar{q}dA-\sum _k\frac{\partial {\bar{T}}_k}{\partial t}\left(\underset{\mathrm{\Omega }}{\int }\rho C_v{N}_j{N}_{k}dV\right)\\ & -\sum _k{\bar{T}}_k(\underset{\mathrm{\Omega }}{\int }\mathrm{\nabla }{N}_j\bullet (\mathit{\kappa }\bullet \mathrm{\nabla }{N}_k)dV)],i,j\in \eta -{\eta }^T,k\in {\eta }^T.\end{array}$

Using the usual argument about the arbitrary nature of $b_j$, we get

$\begin{array}{cc}\sum _i& \frac{\partial T_i}{\partial t}\left(\underset{\mathrm{\Omega }}{\int }\rho C_v{N}_i{N}_{j}dV\right)+\sum _i{T}_i(\underset{\mathrm{\Omega }}{\int }\mathrm{\nabla }{N}_j\bullet (\mathit{\kappa }\bullet \mathrm{\nabla }{N}_i)dV)=\\ & \underset{\mathrm{\Omega }}{\int }{N}_{j}sdV-\underset{{\mathrm{\Gamma }}_q}{\int }{N}_j\bar{q}dA-\sum _k\frac{\partial {\bar{T}}_k}{\partial t}\left(\underset{\mathrm{\Omega }}{\int }\rho C_v{N}_j{N}_{k}dV\right)\\ & -\sum _k{\bar{T}}_k(\underset{\mathrm{\Omega }}{\int }\mathrm{\nabla }{N}_j\bullet (\mathit{\kappa }\bullet \mathrm{\nabla }{N}_k)dV),i,j\in \eta -{\eta }^T,k\in {\eta }^T.\end{array}$

In compact notation,

$\sum _i⟨N_j,\rho C_v{N}_i⟩\dot{T_i}+\sum _{i}a(N_j,N_i)T_i=⟨N_j,s⟩-⟨N_j,\bar{q}{⟩}_{\mathrm{\Gamma }}-\sum _k⟨N_j,\rho C_v{N}_k⟩{\dot{\bar{T}}}_k-\sum _{k}a(N_j,N_k){\bar{T}}_k$

where $i,j\in \eta -{\eta }^T$, and $k\in {\eta }^T$.

Let us rename parts of the above equation as follows:

$\begin{array}{cc}{M}_{ij}& :=⟨N_i,\rho C_v{N}_j⟩\\ K_{ij}& :=a(N_i,N_j)\\ f_j& :=⟨N_j,s⟩-⟨N_j,\bar{q}{⟩}_{\mathrm{\Gamma }}-\sum _k⟨N_j,\rho C_v{N}_k⟩{\dot{\bar{T}}}_k-\sum _{k}a(N_j,N_k){\bar{T}}_k.\end{array}$

Then we get

$\sum _i{M}_{ji}\dot{T_i}+\sum _i{K}_{ji}{T}_i=f_j.$

In matrix form,

$\text{(45)}𝐌\dot{𝐓}+𝐊𝐓=𝐟.$

As before, we can break up the global matrices into sums over elements. Thus,

$\begin{array}{cc}𝐌& =\sum _{e=1}^E{𝐌}^e\\ 𝐊& =\sum _{e=1}^E{𝐊}^e\\ 𝐟& =\sum _{e=1}^E{𝐟}^e\end{array}$

where

$\begin{array}{cc}{M}_{ij}^e& =\underset{{\mathrm{\Omega }}^e}{\int }\rho C_v{N}_i{N}_{j}dV,\\ K_{ij}^e& =\underset{{\mathrm{\Omega }}^e}{\int }\mathrm{\nabla }{N}_i\bullet (\mathit{\kappa }\bullet \mathrm{\nabla }{N}_j)dV\equiv \underset{{\mathrm{\Omega }}^e}{\int }{𝐁}_i^T𝐃{𝐁}_{j}dV,\\ f_j^e& =\underset{{\mathrm{\Omega }}^e}{\int }{N}_{j}sdV-\underset{{\mathrm{\Gamma }}_q^e}{\int }{N}_j\bar{q}dA-\sum _{k=1}^{n^e}[M_{jk}^e{\dot{\bar{T}}}_k-K_{jk}^e{\bar{T}}_k]\end{array}$

and $n^e$ is the number of nodes in an element.

We can do the same for the initial conditions. However, a more common approach is to specify the values of ${𝐓}_0$ directly at the nodes instead of going through an assembly process.

Computing M, K, f

The finite element system of equations at time $t$:

$\begin{array}{cc}{M}_{ij}^e& =\underset{{\mathrm{\Omega }}^e}{\int }\rho C_v{N}_i{N}_{j}dV,\\ K_{ij}^e& =\underset{{\mathrm{\Omega }}^e}{\int }\mathrm{\nabla }{N}_i\bullet (\mathit{\kappa }\bullet \mathrm{\nabla }{N}_j)dV\equiv \underset{{\mathrm{\Omega }}^e}{\int }{𝐁}_i^T𝐃{𝐁}_{j}dV,\\ f_j^e& =\underset{{\mathrm{\Omega }}^e}{\int }{N}_{j}sdV-\underset{{\mathrm{\Gamma }}_q^e}{\int }{N}_j\bar{q}dA-\sum _{k=1}^{n^e}[M_{jk}^e{\dot{\bar{T}}}_k-K_{jk}^e{\bar{T}}_k]\end{array}$

Isoparametric Map

File:Isoparametric.png Isoparametric map.

Coordinate Transformation

$\begin{array}{cc}x& =x_1^e{N}_1^e(\xi ,\eta )+x_2^e{N}_2^e(\xi ,\eta )+x_3^e{N}_3^e(\xi ,\eta )+x_4^e{N}_4^e(\xi ,\eta )\\ & \\ y& =y_1^e{N}_1^e(\xi ,\eta )+y_2^e{N}_2^e(\xi ,\eta )+y_3^e{N}_3^e(\xi ,\eta )+y_4^e{N}_4^e(\xi ,\eta )\end{array}$

Integrating Stiffness Matrix

Stiffness matrix terms:

$K_{ij}^e=\underset{{\mathrm{\Omega }}^e}{\int }\mathrm{\nabla }{N}_i\bullet (\mathit{\kappa }\bullet \mathrm{\nabla }{N}_j)dV$

Index notation:

$K_{ij}^e=\underset{{\mathrm{\Omega }}^e}{\int }{N}_{i,m}{\kappa }_{mn}{N}_{j,n}dV$

Expanded out (in 2D) (assume ${\kappa }_{12}=0$):

$K_{ij}^e=\underset{{\mathrm{\Omega }}^e}{\int }(\frac{\partial N_i}{\partial x}{\kappa }_{11}\frac{\partial N_j}{\partial y}+\frac{\partial N_i}{\partial x}{\kappa }_{22}\frac{\partial N_j}{\partial y})dxdy$

Transformation

Chain Rule:

$\begin{array}{cc}\frac{\partial N_i}{\partial \xi }& =\frac{\partial N_i}{\partial x}\frac{\partial x}{\partial \xi }+\frac{\partial N_i}{\partial y}\frac{\partial y}{\partial \xi }\\ \frac{\partial N_i}{\partial \eta }& =\frac{\partial N_i}{\partial x}\frac{\partial x}{\partial \eta }+\frac{\partial N_i}{\partial y}\frac{\partial y}{\partial \eta }\\ \end{array}$

Matrix notation:

$\left[\begin{array}{c}\frac{\partial N_i}{\partial \xi }\\ \\ \frac{\partial N_i}{\partial \eta }\end{array}\right]=\left[\begin{array}{ccc}\frac{\partial x}{\partial \xi }& & \frac{\partial y}{\partial \xi }\\ \\ \frac{\partial x}{\partial \eta }& & \frac{\partial y}{\partial \eta }\end{array}\right]\left[\begin{array}{c}\frac{\partial N_i}{\partial x}\\ \\ \frac{\partial N_i}{\partial y}\end{array}\right]={𝐉}^e\left[\begin{array}{c}\frac{\partial N_i}{\partial x}\\ \\ \frac{\partial N_i}{\partial y}\end{array}\right]$

Inverse Transformation

$\left[\begin{array}{c}\frac{\partial N_i}{\partial x}\\ \\ \frac{\partial N_i}{\partial y}\end{array}\right]=({𝐉}^e{)}^{-1}\left[\begin{array}{c}\frac{\partial N_i}{\partial \xi }\\ \\ \frac{\partial N_i}{\partial \eta }\end{array}\right];{𝐉}^{*}=({𝐉}^e{)}^{-1}\text{exists if}\det ({𝐉}^e)\ne 0.$

with

$\begin{array}{cc}\frac{\partial x}{\partial \xi }& =x_1^e\frac{\partial N_1^e}{\partial \xi }+x_2^e\frac{\partial N_2^e}{\partial \xi }+x_3^e\frac{\partial N_3^e}{\partial \xi }+x_4^e\frac{\partial N_4^e}{\partial \xi }\\ \frac{\partial y}{\partial \xi }& =y_1^e\frac{\partial N_1^e}{\partial \xi }+y_2^e\frac{\partial N_2^e}{\partial \xi }+y_3^e\frac{\partial N_3^e}{\partial \xi }+y_4^e\frac{\partial N_4^e}{\partial \xi }\\ \frac{\partial x}{\partial \eta }& =x_1^e\frac{\partial N_1^e}{\partial \eta }+x_2^e\frac{\partial N_2^e}{\partial \eta }+x_3^e\frac{\partial N_3^e}{\partial \eta }+x_4^e\frac{\partial N_4^e}{\partial \eta }\\ \frac{\partial y}{\partial \eta }& =y_1^e\frac{\partial N_1^e}{\partial \eta }+y_2^e\frac{\partial N_2^e}{\partial \eta }+y_3^e\frac{\partial N_3^e}{\partial \eta }+y_4^e\frac{\partial N_4^e}{\partial \eta }\end{array}$

Returning to the integration of the stiffness matrix terms:

$K_{ij}^e=\underset{{\mathrm{\Omega }}^e}{\int }(\frac{\partial N_i}{\partial x}{\kappa }_{11}\frac{\partial N_j}{\partial y}+\frac{\partial N_i}{\partial x}{\kappa }_{22}\frac{\partial N_j}{\partial y})dxdy$

In parent element coordinates:

$\begin{array}{cc}{K}_{ij}^e& =\underset{-1}{\overset{1}{\int }}({\kappa }_{11}(J_{11}^{*}\frac{\partial N_i}{\partial \xi }+J_{12}^{*}\frac{\partial N_i}{\partial \eta })+(J_{11}^{*}\frac{\partial N_j}{\partial \xi }+J_{12}^{*}\frac{\partial N_j}{\partial \eta })+\\ & {\kappa }_{22}(J_{21}^{*}\frac{\partial N_i}{\partial \xi }+J_{22}^{*}\frac{\partial N_i}{\partial \eta })+(J_{21}^{*}\frac{\partial N_j}{\partial \xi }+J_{22}^{*}\frac{\partial N_j}{\partial \eta }))d\xi d\eta \end{array}$

Gaussian Quadrature

$\underset{-1}{\overset{1}{\int }}{F}_{ij}(\xi ,\eta )d\xi d\eta =\sum _{I=1}^M\sum _{J=1}^N{F}_{IJ}({\xi }_I,{\eta }_J)W_I{W}_J$

Gaussian Quadrature Example

Find the constants Co, C1, and X so that the quadrature formula

$\underset{0}{\overset{1}{\int }}f(x)dx=Cof(0)+C1f(x_1).$

This has the highest possible degree of precision.

Solution.

Since there are three unknowns, Co, C1 and X1, we will expect the formula to be exact for

$f(x)=1,x,and{x}^2$

Thus

$f(x)=1,\underset{0}{\overset{1}{\int }}f(x)dx=1=C_0+C_{1}Equation1,f(x)=x,\underset{0}{\overset{1}{\int }}f(x)dx=\frac{1}{2}=C_1{x}_{1}Equation2.f(x)=x^2,\underset{0}{\overset{1}{\int }}f(x)dx=\frac{1}{3}=C_1{x}_1^2.$

Equation 2 and 3 will yield.

$\frac{c_1{x}_1}{c_1{x}_1^2}=x_1=\frac{2}{3}.c_1=\frac{3}{4}.c_0=\frac{1}{4}.$

Hence

$\underset{0}{\overset{1}{\int }}f(x)dx=\frac{1}{4}f(0)+\frac{3}{4}f(\frac{2}{3})$

Now,

$f(x)=x^3.\underset{0}{\overset{1}{\int }}{x}^{3}dx=\frac{1}{4}$

And

$\frac{1}{4}(0{)}^3+\frac{3}{4}(\frac{2}{3}{)}^3=\frac{2}{9}.$

Thus the degree of the precision is 2

Robin boundary conditions

A similar approach can be used when w:Robin boundary conditions are needed to be applied. Let us assume that the Robin boundary conditions take the form

$\alpha T+\beta (\mathrm{\nabla }T\cdot 𝐧)=h\text{on}\mathrm{\Gamma }.$

Recall from the previous section that the boundary term in the weak form of the heat equation is

$\underset{\mathrm{\Gamma }}{\int }w\bar{q}dA$

where

$\bar{q}=(\mathit{\kappa }\cdot \mathrm{\nabla }T)\cdot 𝐧.$.

Let us assume, for simplicity, that the material is isotropic. In that case $\mathit{\kappa }$ becomes a scalar quantity $\kappa$ and we can write

$\bar{q}=\kappa (\mathrm{\nabla }T\cdot 𝐧).$

Now, we can write the Robin boundary condition as

$\mathrm{\nabla }T\cdot 𝐧=\frac{h}{\beta }-\frac{\alpha }{\beta }T\equiv \hat{h}-\hat{\alpha }T\text{on}\mathrm{\Gamma }.$

Using this expression we can get of the flux term in $\bar{q}$ to get

$\bar{q}=\kappa (\hat{h}-\hat{\alpha }T).$

Then the boundary term takes the form

$\underset{\mathrm{\Gamma }}{\int }w\bar{q}dA=\underset{\mathrm{\Gamma }}{\int }w\kappa (\hat{h}-\hat{\alpha }T)dA.$

If we ignore the contributions to the trial function from the Dirichlet boundary conditions we can write

$w_h(𝐱)=\sum _{j=1}^n{b}_j{N}_j(𝐱);T_h=v_h(𝐱,t)=\sum _{i=1}^n{T}_i(t)N_i(𝐱)$

Plugging these expressions into the boundary term leads to

$\begin{array}{cc}\sum _j{b}_j\underset{\mathrm{\Gamma }}{\int }{N}_j\bar{q}dA& =\sum _j{b}_j\underset{\mathrm{\Gamma }}{\int }{N}_j\kappa (\hat{h}-\hat{\alpha }\sum _i{T}_i{N}_i)dA\\ & =\sum _j{b}_j(\underset{\mathrm{\Gamma }}{\int }\kappa \hat{h}{N}_{j}dA-\sum _i{T}_i\underset{\mathrm{\Gamma }}{\int }\kappa \hat{\alpha }{N}_i{N}_{j}dA)\end{array}$

Now recall that spatially discretized equation for transient heat conduction can be expressed as (in simplified form; after getting rid of contributions due to Dirichlet boundary conditions and with a scalar $\kappa$)

$\begin{array}{cc}\sum _j{b}_j& [\sum _i\frac{\partial T_i}{\partial t}\left(\underset{\mathrm{\Omega }}{\int }\rho C_v{N}_i{N}_{j}dV\right)+\sum _i\kappa T_i(\underset{\mathrm{\Omega }}{\int }\mathrm{\nabla }{N}_j\bullet \mathrm{\nabla }{N}_{i}dV)]=\\ & \sum _j{b}_j[\underset{\mathrm{\Omega }}{\int }{N}_{j}sdV-\underset{\mathrm{\Gamma }}{\int }{N}_j\bar{q}dA]\end{array}$

Substituting the last term on the right hand side with the new expression for the boundary term, we have

$\begin{array}{cc}\sum _j{b}_j& [\sum _i\frac{\partial T_i}{\partial t}\left(\underset{\mathrm{\Omega }}{\int }\rho C_v{N}_i{N}_{j}dV\right)+\sum _i\kappa T_i(\underset{\mathrm{\Omega }}{\int }\mathrm{\nabla }{N}_j\bullet \mathrm{\nabla }{N}_{i}dV)]=\\ & \sum _j{b}_j[\underset{\mathrm{\Omega }}{\int }{N}_{j}sdV-\kappa \underset{\mathrm{\Gamma }}{\int }\hat{h}{N}_{j}dA+\sum _i\kappa T_i\underset{\mathrm{\Gamma }}{\int }\hat{\alpha }{N}_i{N}_{j}dA]\end{array}$

Invoking the arbitrariness of $b_j$ and rearranging, we get

$\begin{array}{cc}\sum _i\frac{\partial T_i}{\partial t}\left(\underset{\mathrm{\Omega }}{\int }\rho C_v{N}_i{N}_{j}dV\right)& +\sum _i\kappa T_i(\underset{\mathrm{\Omega }}{\int }\mathrm{\nabla }{N}_j\bullet \mathrm{\nabla }{N}_{i}dV-\underset{\mathrm{\Gamma }}{\int }\hat{\alpha }{N}_i{N}_{j}dA)=\\ & \underset{\mathrm{\Omega }}{\int }s{N}_{j}dV-\underset{\mathrm{\Gamma }}{\int }\kappa \hat{h}{N}_{j}dA.\end{array}$

Define

$\begin{array}{cc}{M}_{ji}& :=\underset{\mathrm{\Omega }}{\int }\rho C_v{N}_j{N}_{i}dV\\ K_{ji}& :=\underset{\mathrm{\Omega }}{\int }\mathrm{\nabla }{N}_j\bullet \mathrm{\nabla }{N}_{i}dV-\underset{\mathrm{\Gamma }}{\int }\hat{\alpha }{N}_j{N}_{i}dA\\ f_j& :=\underset{\mathrm{\Omega }}{\int }s{N}_{j}dV-\underset{\mathrm{\Gamma }}{\int }\kappa \hat{h}{N}_{j}dA.\end{array}$

Then the finite element system of equations is

$\sum _i{M}_{ji}{\dot{T}}_i+\sum _i{K}_{ji}{T}_i=f_j.$

or, in matrix form

$𝐌\dot{𝐓}+𝐊𝐓=𝐟.$

Note that the $𝐊$ matrix and the $𝐟$ vector are different from the case where there are no Robin boundary conditions. The boundary terms in the discretized system of equations can be computed in the usual manner and will add terms to the diagonal of the $𝐊$ matrix.