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{align} & \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{align} }$

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{align} \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{align}

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{align} & \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{align} }$

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{align} \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{align}

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

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

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

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

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{align} 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{align}

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{align} \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{align}

where

\begin{align} M^e_{ij} & = \int_{\Omega^e} \rho~C_v~N_i~N_j~dV~,\\ K^e_{ij} & = \int_{\Omega^e}\boldsymbol{\nabla} N_i\bullet(\boldsymbol{\kappa}\bullet\boldsymbol{\nabla} N_j)~dV \equiv \int_{\Omega^e} \mathbf{B}^T_i~\mathbf{D}~\mathbf{B}_j~dV ~, \\ f^e_j & = \int_{\Omega^e} N_j~s~dV - \int_{\Gamma^e_q} N_j~\overline{q}~dA - \sum_{k=1}^{n^e} \left[ M^e_{jk} \dot{\overline{T}}_k - K^e_{jk} \overline{T}_k \right] \end{align}

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 $\mathbf{M}$ , $\mathbf{K}$ , $\mathbf{f}$

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

\begin{align} M^e_{ij} & = \int_{\Omega^e} \rho~C_v~N_i~N_j~dV~,\\ K^e_{ij} & = \int_{\Omega^e}\boldsymbol{\nabla} N_i\bullet({\boldsymbol{\kappa}\bullet{\boldsymbol{\nabla} N_j})}~dV \equiv \int_{\Omega^e} \mathbf{B}^T_i~\mathbf{D}~\mathbf{B}_j~dV ~, \\ f^e_j & = \int_{\Omega^e} N_j~s~dV - \int_{\Gamma^e_q} N_j~\overline{q}~dA - \sum_{k=1}^{n^e} \left[ M^e_{jk} \dot{\overline{T}}_k - K^e_{jk} \overline{T}_k \right] \end{align}

Isoparametric Map

Isoparametric map.

Coordinate Transformation

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

Integrating Stiffness Matrix

Stiffness matrix terms:

K^e_{ij} = \int_{\Omega^e}\boldsymbol{\nabla} N_i\bullet(\boldsymbol{\kappa}\bullet\boldsymbol{\nabla} N_j)~dV

Index notation:

K^e_{ij} = \int_{\Omega^e} N_{i,m}~\kappa_{mn}~N_{j,n}~dV

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

K^e_{ij} = \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{align} \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{align}

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) \ne 0 ~.

with

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

Returning to the integration of the stiffness matrix terms:

K^e_{ij} = \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{align} K^e_{ij} & = \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{align}

Gaussian Quadrature

\int_{-1}^1 F_{ij}(\xi,\eta)~d\xi~d\eta = \sum_{I=1}^M\sum_{J=1}^NF_{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 = Co f(0)+ C1 f(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 \ Equation 1, \qquad f (x)= x,\ \int_{0}^1 f (x)\,dx= \frac{1}{2} = C_1x_1\ \qquad \ Equation 2. \qquad f (x)= x^2,\ \int_{0}^1 f (x)\,dx= \frac{1}{3} = C_1x_1^2.

Equation 2 and 3 will yield.

\frac{c_1x_1}{c_1x_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{align} \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{align}

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{align} \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{align}

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

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

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

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

Define

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

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.

This article is issued from Wikiversity - version of the Friday, October 26, 2012. The text is available under the Creative Commons Attribution/Share Alike but additional terms may apply for the media files.