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{aligned}&{\mathsf {Variational~IBVP~for~the~Heat~Equation}}\\&\\{\text{Find a function}}&~T(t)\in {\mathcal {S}}_{t},t\in [0,\tau ]~{\text{such that for all}}~w\in {\mathcal {V}}\\&\\&\int _{\Omega }\left(\rho ~C_{v}~{\frac {\partial T}{\partial t}}\right)~w~dV+\int _{\Omega }({\boldsymbol {\nabla }}w)\bullet ({\boldsymbol {\kappa }}\bullet {\boldsymbol {\nabla }}T)~dV=\int _{\Omega }s~w~dV-\int _{\Gamma _{q}}w~{\overline {q}}~dA\\&\int _{\Omega }w~\rho ~C_{v}~T(0)~dV=\int _{\Omega }w~\rho ~C_{v}~T_{0}~dV~.\end{aligned}}$

Let ${\mathcal {V}}_{h}$ be a finite dimensional approximation to ${\mathcal {V}}$ and let ${\mathcal {S}}_{h}$ be a finite dimensional approximation to ${\mathcal {S}}$. Let the weighting functions $w_{h}\in {\mathcal {V}}_{h}$ be such that $w_{h}=0$ on $\Gamma _{T}$. Similarly, any other function $v_{h}\in {\mathcal {V}}_{h}$ also goes to zero on the boundary $\Gamma _{T}$.

Assume that the trial solutions $T_{h}\in {\mathcal {S}}_{h}$ can be represented as

$T_{h}=v_{h}+{\overline {T}}_{h}~~~~~{\text{where}}~v_{h}\in {\mathcal {V}}_{h}~.$

The additional quantity ${\overline {T}}_{h}$ results in the satisfaction of the boundary condition $T={\overline {T}}$ on $\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

$\int _{\Omega }\left(\rho ~C_{v}~{\frac {\partial T_{h}}{\partial t}}\right)~w_{h}~dV+\int _{\Omega }({\boldsymbol {\nabla }}w_{h})\bullet ({\boldsymbol {\kappa }}\bullet {\boldsymbol {\nabla }}T_{h})~dV=\int _{\Omega }s~w_{h}~dV-\int _{\Gamma _{q}}w_{h}~{\overline {q}}~dA~.$

Replacing $T_{h}$ with $(v_{h}+{\overline {T}}_{h})$ we get

$\int _{\Omega }\left(\rho ~C_{v}~{\frac {\partial (v_{h}+{\overline {T}}_{h})}{\partial t}}\right)~w_{h}~dV+\int _{\Omega }({\boldsymbol {\nabla }}w_{h})\bullet ({\boldsymbol {\kappa }}\bullet {\boldsymbol {\nabla }}(v_{h}+{\overline {T}}_{h}))~dV=\int _{\Omega }s~w_{h}~dV-\int _{\Gamma _{q}}w_{h}~{\overline {q}}~dA~.$

After expanding and rearranging terms, we get

${\begin{aligned}\int _{\Omega }w_{h}~\rho ~C_{v}~{\frac {\partial v_{h}}{\partial t}}~dV+\int _{\Omega }({\boldsymbol {\nabla }}w_{h})\bullet ({\boldsymbol {\kappa }}\bullet {\boldsymbol {\nabla }}v_{h})~dV=&\int _{\Omega }w_{h}~s~dV-\int _{\Gamma _{q}}w_{h}~{\overline {q}}~dA\\&-\int _{\Omega }w_{h}~\rho ~C_{v}~{\frac {\partial {\overline {T}}_{h}}{\partial t}}~dV-\int _{\Omega }({\boldsymbol {\nabla }}w_{h})\bullet ({\boldsymbol {\kappa }}\bullet {\boldsymbol {\nabla }}{\overline {T}}_{h})~dV~.\end{aligned}}$

In compact notation

$\langle w_{h},~\rho ~C_{v}~{\dot {v}}\rangle _{h}+a(w_{h},~v_{h})=\langle w_{h},~s\rangle -\langle w_{h},~{\overline {q}}\rangle _{\Gamma }-\langle w_{h},~\rho ~C_{v}~{\dot {\overline {T}}}\rangle _{h}-a(w_{h},~{\overline {T}}_{h})~.$

Similarly, the initial condition can be written as

$\int _{\Omega }w_{h}~\rho ~C_{v}~T_{h}(0)~dV=\int _{\Omega }w_{h}~\rho ~C_{v}~T_{0}~dV~.$

Substituting for $T_{h}$ and expanding out terms, we have

$\int _{\Omega }w_{h}~\rho ~C_{v}~v_{h}(0)~dV=\int _{\Omega }w_{h}~\rho ~C_{v}~T_{0}~dV-\int _{\Omega }w_{h}~\rho ~C_{v}~{\overline {T}}_{h}(0)~dV~.$

In compact form,

$\langle w_{h},~\rho ~C_{v}~v_{h}(0)\rangle =\langle w_{h},~\rho ~C_{v}~T_{0}\rangle -\langle w_{h},~\rho ~C_{v}~{\overline {T}}_{h}(0)\rangle ~.$

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

${\text{(40)}}\qquad {\begin{aligned}&{\mathsf {Approximate~Variational~IBVP~for~the~Heat~Equation}}\\&\\{\text{Find a function}}&~T_{h}(t)=v_{h}+{\overline {T}}_{h}\in {\mathcal {S}}_{h},t\in [0,\tau ]~{\text{such that for all}}~w_{h}\in {\mathcal {V}}_{h}\\&\\\int _{\Omega }w_{h}~\rho ~C_{v}~{\frac {\partial v_{h}}{\partial t}}~dV&+\int _{\Omega }({\boldsymbol {\nabla }}w_{h})\bullet ({\boldsymbol {\kappa }}\bullet {\boldsymbol {\nabla }}v_{h})~dV=\int _{\Omega }w_{h}~s~dV-\int _{\Gamma _{q}}w_{h}~{\overline {q}}~dA\\&-\int _{\Omega }w_{h}~\rho ~C_{v}~{\frac {\partial {\overline {T}}_{h}}{\partial t}}~dV-\int _{\Omega }({\boldsymbol {\nabla }}w_{h})\bullet ({\boldsymbol {\kappa }}\bullet {\boldsymbol {\nabla }}{\overline {T}}_{h})~dV~.\\\int _{\Omega }w_{h}~\rho ~C_{v}~v_{h}(0)~dV&=\int _{\Omega }w_{h}~\rho ~C_{v}~T_{0}~dV-\int _{\Omega }w_{h}~\rho ~C_{v}~{\overline {T}}_{h}(0)~dV~.\end{aligned}}$

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 $\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.

Figure 1. FEM discretization for the heat conduction problem.

As we have seen before, a typical weighting function $w_{h}\in {\mathcal {V}}_{h}$ is assumed to have the form

${\text{(42)}}\qquad {w_{h}(\mathbf {x} )=\sum _{i=1}^{n}b_{i}N_{i}(\mathbf {x} )\equiv \sum _{i\in \eta -\eta ^{T}}b_{i}N_{i}(\mathbf {x} )~.}$

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 {\mathcal {V}}_{h}$ can be written as

${\text{(43)}}\qquad {v_{h}(\mathbf {x} ,t)=\sum _{i\in \eta -\eta ^{T}}T_{i}(t)N_{i}(\mathbf {x} )}$

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)}}\qquad {{\overline {T}}_{h}(\mathbf {x} ,t)=\sum _{i\in \eta ^{T}}{\overline {T}}_{i}(t)N_{i}(\mathbf {x} ),\qquad {\overline {T}}_{i}(t)={\overline {T}}(\mathbf {x} _{i},t)~.}$

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

${\begin{aligned}\int _{\Omega }\left(\sum _{j}b_{j}N_{j}\right)~&\rho ~C_{v}~\left(\sum _{i}{\frac {\partial T_{i}}{\partial t}}N_{i}\right)~dV+\int _{\Omega }\left(\sum _{j}b_{j}{\boldsymbol {\nabla }}N_{j}\right)\bullet \left[\sum _{i}T_{i}\left({\boldsymbol {\kappa }}\bullet {\boldsymbol {\nabla }}N_{i}\right)\right]~dV=\\&\int _{\Omega }\left(\sum _{j}b_{j}N_{j}\right)~s~dV-\int _{\Gamma _{q}}\left(\sum _{j}b_{j}N_{j}\right)~{\overline {q}}~dA\\&-\int _{\Omega }\left(\sum _{j}b_{j}N_{j}\right)~\rho ~C_{v}~\left(\sum _{k}{\frac {\partial {\overline {T}}_{k}}{\partial t}}N_{k}\right)~dV\\&-\int _{\Omega }\left(\sum _{j}b_{j}{\boldsymbol {\nabla }}N_{j}\right)\bullet \left[\sum _{k}{\overline {T}}_{k}\left({\boldsymbol {\kappa }}\bullet {\boldsymbol {\nabla }}N_{k}\right)\right]~dV,\qquad i,j\in \eta -\eta ^{T},\quad k\in \eta ^{T}~.\end{aligned}}$

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

${\begin{aligned}\sum _{j}b_{j}&\left[\sum _{i}{\frac {\partial T_{i}}{\partial t}}\left(\int _{\Omega }\rho ~C_{v}~N_{i}~N_{j}~dV\right)+\sum _{i}T_{i}\left(\int _{\Omega }{\boldsymbol {\nabla }}N_{j}\bullet ({\boldsymbol {\kappa }}\bullet {\boldsymbol {\nabla }}N_{i})~dV\right)\right]=\\&\sum _{j}b_{j}\left[\int _{\Omega }N_{j}~s~dV-\int _{\Gamma _{q}}N_{j}~{\overline {q}}~dA-\sum _{k}{\frac {\partial {\overline {T}}_{k}}{\partial t}}\left(\int _{\Omega }\rho ~C_{v}~N_{j}~N_{k}~dV\right)\right.\\&-\left.\sum _{k}{\overline {T}}_{k}\left(\int _{\Omega }{\boldsymbol {\nabla }}N_{j}\bullet ({\boldsymbol {\kappa }}\bullet {\boldsymbol {\nabla }}N_{k})~dV\right)\right],\qquad i,j\in \eta -\eta ^{T},\quad k\in \eta ^{T}~.\end{aligned}}$

Using the usual argument about the arbitrary nature of $b_{j}$, we get

${\begin{aligned}\sum _{i}&{\frac {\partial T_{i}}{\partial t}}\left(\int _{\Omega }\rho ~C_{v}~N_{i}~N_{j}~dV\right)+\sum _{i}T_{i}\left(\int _{\Omega }{\boldsymbol {\nabla }}N_{j}\bullet ({\boldsymbol {\kappa }}\bullet {\boldsymbol {\nabla }}N_{i})~dV\right)=\\&\int _{\Omega }N_{j}~s~dV-\int _{\Gamma _{q}}N_{j}~{\overline {q}}~dA-\sum _{k}{\frac {\partial {\overline {T}}_{k}}{\partial t}}\left(\int _{\Omega }\rho ~C_{v}~N_{j}~N_{k}~dV\right)\\&-\sum _{k}{\overline {T}}_{k}\left(\int _{\Omega }{\boldsymbol {\nabla }}N_{j}\bullet ({\boldsymbol {\kappa }}\bullet {\boldsymbol {\nabla }}N_{k})~dV\right),\qquad i,j\in \eta -\eta ^{T},\quad k\in \eta ^{T}~.\end{aligned}}$

In compact notation,

${\sum _{i}\langle N_{j},~\rho C_{v}N_{i}\rangle ~{\dot {T_{i}}}+\sum _{i}a(N_{j},N_{i})~T_{i}=\langle N_{j},~s\rangle -\langle N_{j},~{\overline {q}}\rangle _{\Gamma }-\sum _{k}\langle N_{j},~\rho C_{v}N_{k}\rangle ~{\dot {\overline {T}}}_{k}-\sum _{k}a(N_{j},N_{k})~{\overline {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{aligned}M_{ij}&:=\langle N_{i},~\rho C_{v}N_{j}\rangle \\K_{ij}&:=a(N_{i},N_{j})\\f_{j}&:=\langle N_{j},~s\rangle -\langle N_{j},~{\overline {q}}\rangle _{\Gamma }-\sum _{k}\langle N_{j},~\rho C_{v}N_{k}\rangle ~{\dot {\overline {T}}}_{k}-\sum _{k}a(N_{j},N_{k})~{\overline {T}}_{k}~.\end{aligned}}$

Then we get

${\sum _{i}M_{ji}~{\dot {T_{i}}}+\sum _{i}K_{ji}~T_{i}=f_{j}~.}$

In matrix form,

${\text{(45)}}\qquad {\mathbf {M} ~{\dot {\mathbf {T} }}+\mathbf {K} ~\mathbf {T} =\mathbf {f} ~.}$

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

${\begin{aligned}\mathbf {M} &=\sum _{e=1}^{E}\mathbf {M} ^{e}\\\mathbf {K} &=\sum _{e=1}^{E}\mathbf {K} ^{e}\\\mathbf {f} &=\sum _{e=1}^{E}\mathbf {f} ^{e}\end{aligned}}$

where

${\begin{aligned}M_{ij}^{e}&=\int _{\Omega ^{e}}\rho ~C_{v}~N_{i}~N_{j}~dV~,\\K_{ij}^{e}&=\int _{\Omega ^{e}}{\boldsymbol {\nabla }}N_{i}\bullet ({\boldsymbol {\kappa }}\bullet {\boldsymbol {\nabla }}N_{j})~dV\equiv \int _{\Omega ^{e}}\mathbf {B} _{i}^{T}~\mathbf {D} ~\mathbf {B} _{j}~dV~,\\f_{j}^{e}&=\int _{\Omega ^{e}}N_{j}~s~dV-\int _{\Gamma _{q}^{e}}N_{j}~{\overline {q}}~dA-\sum _{k=1}^{n^{e}}\left[M_{jk}^{e}{\dot {\overline {T}}}_{k}-K_{jk}^{e}{\overline {T}}_{k}\right]\end{aligned}}$

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 $\mathbf {T} _{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{aligned}M_{ij}^{e}&=\int _{\Omega ^{e}}\rho ~C_{v}~N_{i}~N_{j}~dV~,\\K_{ij}^{e}&=\int _{\Omega ^{e}}{\boldsymbol {\nabla }}N_{i}\bullet ({{\boldsymbol {\kappa }}\bullet {{\boldsymbol {\nabla }}N_{j}})}~dV\equiv \int _{\Omega ^{e}}\mathbf {B} _{i}^{T}~\mathbf {D} ~\mathbf {B} _{j}~dV~,\\f_{j}^{e}&=\int _{\Omega ^{e}}N_{j}~s~dV-\int _{\Gamma _{q}^{e}}N_{j}~{\overline {q}}~dA-\sum _{k=1}^{n^{e}}\left[M_{jk}^{e}{\dot {\overline {T}}}_{k}-K_{jk}^{e}{\overline {T}}_{k}\right]\end{aligned}}$

Isoparametric Map

Isoparametric map.

Coordinate Transformation

${\begin{aligned}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{aligned}}$

Integrating Stiffness Matrix

Stiffness matrix terms:

$K_{ij}^{e}=\int _{\Omega ^{e}}{\boldsymbol {\nabla }}N_{i}\bullet ({\boldsymbol {\kappa }}\bullet {\boldsymbol {\nabla }}N_{j})~dV$

Index notation:

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

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

$K_{ij}^{e}=\int _{\Omega ^{e}}\left({\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}}\right)~dx~dy$

Transformation

Chain Rule:

${\begin{aligned}{\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{aligned}}$

Matrix notation:

${\begin{bmatrix}{\frac {\partial N_{i}}{\partial \xi }}\\\\{\frac {\partial N_{i}}{\partial \eta }}\end{bmatrix}}={\begin{bmatrix}{\frac {\partial x}{\partial \xi }}&&{\frac {\partial y}{\partial \xi }}\\\\{\frac {\partial x}{\partial \eta }}&&{\frac {\partial y}{\partial \eta }}\end{bmatrix}}{\begin{bmatrix}{\frac {\partial N_{i}}{\partial x}}\\\\{\frac {\partial N_{i}}{\partial y}}\end{bmatrix}}=\mathbf {J} ^{e}{\begin{bmatrix}{\frac {\partial N_{i}}{\partial x}}\\\\{\frac {\partial N_{i}}{\partial y}}\end{bmatrix}}$

Inverse Transformation

${\begin{bmatrix}{\frac {\partial N_{i}}{\partial x}}\\\\{\frac {\partial N_{i}}{\partial y}}\end{bmatrix}}=(\mathbf {J} ^{e})^{-1}{\begin{bmatrix}{\frac {\partial N_{i}}{\partial \xi }}\\\\{\frac {\partial N_{i}}{\partial \eta }}\end{bmatrix}}~;~~\mathbf {J} ^{*}=(\mathbf {J} ^{e})^{-1}~{\text{exists if}}~\det(\mathbf {J} ^{e})\neq 0~.$

with

${\begin{aligned}{\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{aligned}}$

Returning to the integration of the stiffness matrix terms:

$K_{ij}^{e}=\int _{\Omega ^{e}}\left({\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}}\right)~dx~dy$

In parent element coordinates:

${\begin{aligned}K_{ij}^{e}&=\int _{-1}^{1}\left(\kappa _{11}\left(J_{11}^{*}{\frac {\partial N_{i}}{\partial \xi }}+J_{12}^{*}{\frac {\partial N_{i}}{\partial \eta }}\right)+\left(J_{11}^{*}{\frac {\partial N_{j}}{\partial \xi }}+J_{12}^{*}{\frac {\partial N_{j}}{\partial \eta }}\right)+\right.\\&\left.\kappa _{22}\left(J_{21}^{*}{\frac {\partial N_{i}}{\partial \xi }}+J_{22}^{*}{\frac {\partial N_{i}}{\partial \eta }}\right)+\left(J_{21}^{*}{\frac {\partial N_{j}}{\partial \xi }}+J_{22}^{*}{\frac {\partial N_{j}}{\partial \eta }}\right)\right)~d\xi ~d\eta \end{aligned}}$

Gaussian Quadrature

$\int _{-1}^{1}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

$\int _{0}^{1}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,\ \int _{0}^{1}f(x)\,dx=1=C_{0}+C_{1}\ \qquad \ Equation1,\qquad f(x)=x,\ \int _{0}^{1}f(x)\,dx={\frac {1}{2}}=C_{1}x_{1}\ \qquad \ Equation2.\qquad f(x)=x^{2},\ \int _{0}^{1}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}}.\qquad c_{1}={\frac {3}{4}}.\qquad c_{0}={\frac {1}{4}}.$

Hence

$\int _{0}^{1}f(x)\,dx={\frac {1}{4}}f(0)+{\frac {3}{4}}f({\frac {2}{3}})$

Now,

$f(x)=x^{3}.\qquad \int _{0}^{1}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 ~({\boldsymbol {\nabla }}T\cdot \mathbf {n} )=h\qquad \qquad {\text{on}}~\Gamma ~.$

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

$\int _{\Gamma }w~{\overline {q}}~dA$

where

${\overline {q}}=({\boldsymbol {\kappa }}\cdot {\boldsymbol {\nabla }}T)\cdot \mathbf {n} ~.$.

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

${\overline {q}}=\kappa ~({\boldsymbol {\nabla }}T\cdot \mathbf {n} )~.$

Now, we can write the Robin boundary condition as

${\boldsymbol {\nabla }}T\cdot \mathbf {n} ={\cfrac {h}{\beta }}-{\cfrac {\alpha }{\beta }}~T\equiv {\hat {h}}-{\hat {\alpha }}~T\qquad \qquad {\text{on}}~\Gamma ~.$

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

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

Then the boundary term takes the form

$\int _{\Gamma }w~{\overline {q}}~dA=\int _{\Gamma }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}(\mathbf {x} )=\sum _{j=1}^{n}b_{j}N_{j}(\mathbf {x} )~;~~T_{h}=v_{h}(\mathbf {x} ,t)=\sum _{i=1}^{n}T_{i}(t)N_{i}(\mathbf {x} )$

Plugging these expressions into the boundary term leads to

${\begin{aligned}\sum _{j}b_{j}\int _{\Gamma }N_{j}~{\overline {q}}~dA&=\sum _{j}b_{j}\int _{\Gamma }N_{j}~\kappa ~({\hat {h}}-{\hat {\alpha }}~\sum _{i}T_{i}N_{i})~dA\\&=\sum _{j}b_{j}\left(\int _{\Gamma }\kappa ~{\hat {h}}~N_{j}~dA-\sum _{i}T_{i}~\int _{\Gamma }\kappa ~{\hat {\alpha }}~N_{i}~N_{j}~dA\right)\end{aligned}}$

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{aligned}\sum _{j}b_{j}&\left[\sum _{i}{\frac {\partial T_{i}}{\partial t}}\left(\int _{\Omega }\rho ~C_{v}~N_{i}~N_{j}~dV\right)+\sum _{i}\kappa ~T_{i}\left(\int _{\Omega }{\boldsymbol {\nabla }}N_{j}\bullet {\boldsymbol {\nabla }}N_{i}~dV\right)\right]=\\&\sum _{j}b_{j}\left[\int _{\Omega }N_{j}~s~dV-\int _{\Gamma }N_{j}~{\overline {q}}~dA\right]\end{aligned}}$

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

${\begin{aligned}\sum _{j}b_{j}&\left[\sum _{i}{\frac {\partial T_{i}}{\partial t}}\left(\int _{\Omega }\rho ~C_{v}~N_{i}~N_{j}~dV\right)+\sum _{i}\kappa ~T_{i}\left(\int _{\Omega }{\boldsymbol {\nabla }}N_{j}\bullet {\boldsymbol {\nabla }}N_{i}~dV\right)\right]=\\&\sum _{j}b_{j}\left[\int _{\Omega }N_{j}~s~dV-\kappa ~\int _{\Gamma }{\hat {h}}~N_{j}~dA+\sum _{i}\kappa ~T_{i}~\int _{\Gamma }{\hat {\alpha }}~N_{i}~N_{j}~dA\right]\end{aligned}}$

Invoking the arbitrariness of $b_{j}$ and rearranging, we get

${\begin{aligned}\sum _{i}{\frac {\partial T_{i}}{\partial t}}\left(\int _{\Omega }\rho ~C_{v}~N_{i}~N_{j}~dV\right)&+\sum _{i}\kappa ~T_{i}\left(\int _{\Omega }{\boldsymbol {\nabla }}N_{j}\bullet {\boldsymbol {\nabla }}N_{i}~dV-\int _{\Gamma }{\hat {\alpha }}~N_{i}~N_{j}~dA\right)=\\&\int _{\Omega }s~N_{j}~dV-\int _{\Gamma }\kappa ~{\hat {h}}~N_{j}~dA~.\end{aligned}}$

Define

${\begin{aligned}M_{ji}&:=\int _{\Omega }\rho ~C_{v}~N_{j}~N_{i}~dV\\K_{ji}&:=\int _{\Omega }{\boldsymbol {\nabla }}N_{j}\bullet {\boldsymbol {\nabla }}N_{i}~dV-\int _{\Gamma }{\hat {\alpha }}~N_{j}~N_{i}~dA\\f_{j}&:=\int _{\Omega }s~N_{j}~dV-\int _{\Gamma }\kappa ~{\hat {h}}~N_{j}~dA~.\end{aligned}}$

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

$\mathbf {M} ~{\dot {\mathbf {T} }}+\mathbf {K} ~\mathbf {T} =\mathbf {f} ~.$

Note that the $\mathbf {K}$ matrix and the $\mathbf {f}$ 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 $\mathbf {K}$ matrix.

← Solution of Poisson equation

Finite elements

Heat equation time integration →