4.2.2 Discretization of the Thermal Equation

For convenience, the thermal equation (3.72) is rewritten here as

$$\symMatDens\symSpecHeat\dot\T - \ensuremath{\nabla\cdot{(\symTher... ...abla{\T}})}} - \symElecCond\left(\ensuremath{\nabla{\symElecPot}}\right)^2 = 0,$$ (4.33)

where the notation $\dot\T = \partial{\T}/\partial{t}$ is used. Following the same procedure as for Laplace's equation, one multiplies (4.33) by ${\symShapeFun}_{p}(\vec r)$ and integrates over the domain $\symDomain$, which together with Green's formula (4.26) and a Neumann boundary condition yields the weak formulation

$$\ensuremath{\int_{\symDomain}{\symMatDens\symSpecHeat\dot\T{\symS... ...bla{\symElecPot}}\right)^2 \symShapeFun_p}\ d{\symDomain}} = 0,\quad p=1,...,N.$$ (4.34)

The thermal equation (4.33) corresponds to a parabolic problem. Besides the spatial discretization, a discretization in time has also to be performed. A simple choice is the backward Euler method

$$\dot\T = \frac{\T^n - \T^{n-1}}{\Delta t_n},$$ (4.35)

where $\Delta t_n = t^n - t^{n-1}$ is the time step. Applying now the spatial discretization for the temperature variable

$$\T^n = \ensuremath{\sum_{j=1}^{4}{\T_j^n}}{\symShapeFun}_{j}(\vec r),$$ (4.36)

and the electric potential discretization (4.29), equation (4.34) is written for a single element as

$$\symMatDens\symSpecHeat \sum_{j=1}^{4}\T_j^n \int_{T} \symShapeFu... ...mShapeFun_j}}\ensuremath{\cdot}\ensuremath{\nabla{\symShapeFun_p}}\ d\symDomain$$

$$- \Delta t_n \symElecCond \sum_{i=1}^4\sum_{q=1}^4 \symElecPot_i^... ...nabla{\symShapeFun_q}}\right) \symShapeFun_p\ d\symDomain = 0,\qquad p=1,...,4.$$ (4.37)

Using the integral notations (4.31),

$$M_{jp} = \int_{T}\symShapeFun_j \symShapeFun_p\ d\symDomain = \de... ...\xi}\int_{0}^{1-\xi-\eta} \symShapeFun_j^t \symShapeFun_p^t\ d\zeta d\eta d\xi,$$ (4.38)

and

$$\Theta_{iqp} = \int_{T} \left(\ensuremath{\nabla{\symShapeFun_i}}\ensuremath{\cdot}\ensuremath{\nabla{\symShapeFun_q}}\right) \symShapeFun_p\ d\symDomain$$

$$=\det(\mathbf{J})\int_{0}^{1}\int_{0}^{1-\xi}\int_{0}^{1-\xi-\eta... ...f{\Lambda}\nabla^t \symShapeFun_q^t\right)\symShapeFun_p^t \ d\zeta d\eta d\xi,$$ (4.39)

(4.37) can be expressed as

$$\symMatDens\symSpecHeat \sum_{j=1}^{4}\T_j^n M_{jp} - \symMatDens... ...=1}^{4}\T_j^{n-1} M_{jp} + \Delta t_n \symThermCond \sum_{j=1}^{4}\T_j^n K_{jp}$$

$$- \Delta t_n \symElecCond \sum_{i=1}^4\sum_{q=1}^4 \symElecPot_i^n\symElecPot_q^n \Theta_{iqp} = 0,\qquad p=1,...,4.$$ (4.40)

This is the discrete system of equations, which has to be solved each time step in order to obtain the temperature at each node of the element.