flw2teΒΆ

Purpose:

Compute element stiffness matrix for a triangular heat flow element.

_images/flw2te.svg
Syntax:
Ke = flw2te(ex, ey, ep, D)
[Ke, fe] = flw2te(ex, ey, ep, D, eq)
Description:

flw2te provides the element stiffness (conductivity) matrix Ke and the element load vector fe for a triangular heat flow element.

The element nodal coordinates \(x_1\), \(y_1\), \(x_2\) etc, are supplied to the function by ex and ey, the element thickness \(t\) is supplied by ep and the thermal conductivities (or corresponding quantities) \(k_{xx}\), \(k_{xy}\) etc are supplied by D.

\[\begin{split}\begin{array}{l} \mathbf{ex} = [\, x_1 \;\; x_2 \;\; x_3\,] \\ \mathbf{ey} = [\, y_1 \;\; y_2 \;\; y_3\,] \end{array} \qquad \mathbf{ep} = [\, t \,] \qquad \mathbf{D} = \begin{bmatrix} k_{xx} & k_{xy} \\ k_{yx} & k_{yy} \end{bmatrix}\end{split}\]

If the scalar variable eq is given in the function, the element load vector \(\mathbf{fe}\) is computed, using

\[\mathbf{eq} = [\, Q \,]\]

where \(Q\) is the heat supply per unit volume.

Theory:

The element stiffness matrix \(\mathbf{K}^e\) and the element load vector \(\mathbf{f}_l^e\), stored in Ke and fe, respectively, are computed according to

\[\mathbf{K}^e = (\mathbf{C}^{-1})^T \int_A \bar{\mathbf{B}}^T \mathbf{D} \bar{\mathbf{B}}\, t\, dA\, \mathbf{C}^{-1}\]
\[\mathbf{f}_l^e = (\mathbf{C}^{-1})^T \int_A \bar{\mathbf{N}}^T Q\, t\, dA\]

with the constitutive matrix \(\mathbf{D}\) defined by D.

The evaluation of the integrals for the triangular element is based on the linear temperature approximation \(T(x, y)\) and is expressed in terms of the nodal variables \(T_1\), \(T_2\) and \(T_3\) as

\[T(x, y) = \mathbf{N}^e \mathbf{a}^e = \bar{\mathbf{N}}\, \mathbf{C}^{-1} \mathbf{a}^e\]

where

\[\begin{split}\bar{\mathbf{N}} = [\, 1 \;\; x \;\; y\,] \qquad \mathbf{C} = \begin{bmatrix} 1 & x_1 & y_1 \\ 1 & x_2 & y_2 \\ 1 & x_3 & y_3 \end{bmatrix} \qquad \mathbf{a}^e = \begin{bmatrix} T_1 \\ T_2 \\ T_3 \end{bmatrix}\end{split}\]

and hence it follows that

\[\begin{split}\bar{\mathbf{B}} = \nabla \bar{\mathbf{N}} = \begin{bmatrix} 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} \qquad \nabla = \begin{bmatrix} \dfrac{\partial}{\partial x} \\ \dfrac{\partial}{\partial y} \end{bmatrix}\end{split}\]

Evaluation of the integrals for the triangular element yields

\[\mathbf{K}^e = (\mathbf{C}^{-1})^T \bar{\mathbf{B}}^T \mathbf{D} \bar{\mathbf{B}} \mathbf{C}^{-1} t A\]
\[\begin{split}\mathbf{f}_l^e = \frac{Q A t}{3} \begin{bmatrix} 1 \\ 1 \\ 1 \end{bmatrix}\end{split}\]

where the element area \(A\) is determined as

\[A = \frac{1}{2} \det \mathbf{C}\]