PDF (legal size)

## Review of FEM solution for the torsion problem of a rectangular cross section

Nov 27,2006   Compiled on September 4, 2021 at 3:50pm

### 1 Introduction

This is a review of the FEM solution to the torsion problem of a rectangular cross section beam. First a description of the problem is given, then a description of the FEM method is shown, followed by a simple numerical worked example.

### 2 The problem

The problem is to solve the Poisson 2D problem for rectangular cross section. This equation is the mathematical model for a beam under torsion as described in the following diagram.

Since we are using a triangle elements for the FEM mesh, the cross section mesh is the preferred mesh to use as shown in this diagram.

### 3 The big picture

Before going into the details of the FEM solution it might be useful to look at the big picture.

The following diagram shows more description of the methods.

### 4 Mathematical derivation

#### 4.1 Derivation of the symmetric weak form of the 2D Poisson equation

The following diagram shows the steps to obtain the symmetric weak for of the 2D Poisson PDE

#### 4.2 Converting the symmetric weak form equation from the global Cartesian coordinates system to natural coordinates system

Converting the integral equations from the global Cartesian coordinates system to what is called the natural coordinates system (the local coordinates system) is a standard step used in FEM.

”A local coordinates system that relies on the element geometry for its deﬁnition and whose coordinates range between zero and unity within the element is known a natural coordinates system. Such system have the property that one particular coordinate has unit value at one node of the element and zero value at the other nodes: its variation between nodes is linear”1

Integration of shape functions when they are written in the natural coordinates are simpliﬁed since the origin is now located on the element. These are the main reasons for changing from global coordinates to the natural coordinates. For simple geometries, one can avoid having to do this coordinates transformation, but in general and in practice it is the standard procedure to do.

I found that most of the technical and mathematical diﬃculties involved are in this step. So more details will be spend on this.

The global coordinates of the element is shown in this diagram

Given an equation or expression where the independent variables in the equation are $$x,y$$ (the global Cartesian coordinates system) and we wish to express this same equation using the independent variables $$\zeta ,\eta$$, then we perform coordinates transformations.

Given that $$x=x\left (\zeta ,\eta \right )$$ and $$y=y\left (\zeta ,\eta \right )$$ , we ﬁrst ﬁnd the diﬀerentials of the old coordinates system (i.e. $$dx,dy$$) in terms of the diﬀerentials of the new coordinates system ($$d\zeta ,d\eta$$)

The matrix that represents this mapping between the diﬀerentials in the old coordinates system and the new coordinates system us called the Jacobian (some books call the determinant of this matrix as the Jacobian). It is important to note that this mapping is between the diﬀerentials of the independent variables in the two coordinates system, and not between the variables themselves.

Hence we write $\ \left \{ \begin {array} [c]{c}dx\\ dy \end {array} \right \} =\overset {J}{\overbrace {\left [ \begin {array} [c]{cc}\frac {\partial x}{\partial \zeta } & \frac {\partial x}{\partial \eta }\\ \frac {\partial y}{\partial \zeta } & \frac {\partial y}{\partial \eta }\end {array} \right ] }}\left \{ \begin {array} [c]{c}d\zeta \\ d\eta \end {array} \right \}$ $$J$$ is also written as $J=\frac {\partial \left (x,y\right ) }{\partial \left (\zeta ,\eta \right ) }$ The main use for the Jacobian is in change of variables from one coordinates system to another, and also in performing area and volume integrals.

Hence, converting an integral from the global coordinates to the natural coordinates can be done as follows

$\int _{A} f(x,y) dx dy = \int _{-1}^{1}\int _{-1}^{1} g\left ( \zeta ,\eta \right ) \left \vert J\right \vert d\zeta d\eta$

When the natural coordinates are area coordinates (which is the case here), we should modify the above to become

$\int _{A}f\left (x,y\right ) dxdy=\int _{0}^{1}\int _{0}^{1-L_{1}}g\left ( L_{1},L_{2}\right ) \ \left \vert J\right \vert \ dL_{2}dL_{1}$

The area coordinates $$\left (L_{1},L_{2},L_{3}\right )$$ are illustrated in this diagram

It is important to realize that the shape functions $$N_{1},N_{2},N_{3}$$ used will be the same as the area coordinates.

Let us now start from the symmetric weak form equation, with the goal to convert it to the natural coordinates (see previous diagram for the derivation of this equation)

\begin {equation} {\displaystyle \int \limits _{\Gamma _{external\ boundary}}} \left (\frac {\partial u}{\partial x}v\right ) n_{x}\ +\left (\frac {\partial u}{\partial y}v\right ) n_{y}\ d\Gamma - {\displaystyle \int \limits _{\Omega _{j}}} \frac {\partial u}{\partial x}\frac {\partial v}{\partial x}\ +\frac {\partial u}{\partial y} \frac {\partial v}{\partial y} dx dy + {\displaystyle \int \limits _{\Omega _{j}}} fv dx dy =0\ \tag {1} \end {equation}

Since the ﬁrst integral above is carried along the boundaries of the whole domain itself (not along the boundaries of the individual elements themselves) and since we set the value of the test function $$v$$ to be zero at the boundaries of the domain, the ﬁrst part of the above integral is zero. Hence the above integral become

\begin {equation} -{\displaystyle \int \limits _{\Omega _{j}}} \frac {\partial u}{\partial x}\frac {\partial v}{\partial x}\ +\frac {\partial u}{\partial y}\frac {\partial v}{\partial y} dx dy + {\displaystyle \int \limits _{\Omega _{j}}} f v dx dy =0\ \tag {1} \end {equation}

In the following derivations, everything is done on an element $$j$$, hence all the $$u,v,$$ and element nodes coordinates $$x_{1},y_{1}$$, etc.. should have a superscript $$j$$ on, as in $$u^{j},v^{j}\,\ etc..$$. To make things easier to read, I will not put the superscript $$j$$ but will add it back at the end.

Consider ﬁrst the second integral from (1) Which can be rewritten as \begin {equation} I_{1}=-{\displaystyle \int \limits _{\Omega _{j}}} \left ( \begin {array} [c]{cc}\frac {\partial v}{\partial x} & \frac {\partial v}{\partial y}\end {array} \right ) \left \{ \begin {array} [c]{c}\frac {\partial u}{\partial x}\\ \frac {\partial u}{\partial y}\end {array} \right \} \textbf { dxdy}+{\displaystyle \int \limits _{\Omega _{j}}} fv dx dy \tag {2} \end {equation}

Consider the ﬁrst integral from above

\begin {equation} I={\displaystyle \int \limits _{\Omega _{j}}} \left ( \begin {array} [c]{cc}\frac {\partial v}{\partial x} & \frac {\partial v}{\partial y}\end {array} \right ) \left \{ \begin {array} [c]{c}\frac {\partial u}{\partial x}\\ \frac {\partial u}{\partial y}\end {array} \right \} dx dy \tag {2A} \end {equation}

The above is written with reference to the global coordinates system. However, We want our trial and test functions to be deﬁned in the natural coordinates system (where things are simpler). So we need a way to transform the above integral (2A) to the natural coordinates system.

Assume we have the mapping $$x=x\left (\zeta ,\eta \right )$$ and $$y=y\left ( \zeta ,\eta \right )$$ (we will see how to obtain this mapping below). This mapping tells us how the global coordinates themselves change as a function of the natural coordinates. Now we can use diﬀerentiation chain rule to see how the trial and test functions themselves change relative the global coordinates.

\begin {align*} \frac {\partial u}{\partial \zeta } & =\frac {\partial u}{\partial x}\frac {\partial x}{\partial \zeta }+\frac {\partial u}{\partial y}\frac {\partial y}{\partial \zeta }\\ \frac {\partial u}{\partial \eta } & =\frac {\partial u}{\partial x}\frac {\partial x}{\partial \eta }+\frac {\partial u}{\partial y}\frac {\partial y}{\partial \eta } \end {align*}

Similarly for the test function

\begin {align*} \frac {\partial v}{\partial \zeta } & =\frac {\partial v}{\partial x}\frac {\partial x}{\partial \zeta }+\frac {\partial v}{\partial y}\frac {\partial y}{\partial \zeta }\\ \frac {\partial v}{\partial \eta } & =\frac {\partial v}{\partial x}\frac {\partial x}{\partial \eta }+\frac {\partial v}{\partial y}\frac {\partial y}{\partial \eta } \end {align*}

To make things more clear, we rewrite the above using matrix notation. For the trial function

\begin {align} \left \{ \begin {array} [c]{c}\frac {\partial u}{\partial \zeta }\\ \frac {\partial u}{\partial \eta }\end {array} \right \} & =\left [ \begin {array} [c]{cc}\frac {\partial x}{\partial \zeta } & \frac {\partial y}{\partial \zeta }\\ \frac {\partial x}{\partial \eta } & \frac {\partial y}{\partial \eta }\end {array} \right ] \left \{ \begin {array} [c]{c}\frac {\partial u}{\partial x}\\ \frac {\partial u}{\partial y}\end {array} \right \} \nonumber \\ & =\left [ J\right ] \left \{ \begin {array} [c]{c}\frac {\partial u}{\partial x}\\ \frac {\partial u}{\partial y}\end {array} \right \} \tag {3} \end {align}

and similarly for the test function

\begin {align} \left \{ \begin {array} [c]{c}\frac {\partial v}{\partial \zeta }\\ \frac {\partial v}{\partial \eta }\end {array} \right \} & =\left [ \begin {array} [c]{cc}\frac {\partial x}{\partial \zeta } & \frac {\partial y}{\partial \zeta }\\ \frac {\partial x}{\partial \eta } & \frac {\partial y}{\partial \eta }\end {array} \right ] \left \{ \begin {array} [c]{c}\frac {\partial v}{\partial x}\\ \frac {\partial v}{\partial y}\end {array} \right \} \nonumber \\ & =\left [ J\right ] \left \{ \begin {array} [c]{c}\frac {\partial v}{\partial x}\\ \frac {\partial v}{\partial y}\end {array} \right \} \tag {4} \end {align}

From (3) and (4), we see the following inverse transformations

\begin {equation} \left \{ \begin {array} [c]{c}\frac {\partial u}{\partial x}\\ \frac {\partial u}{\partial y}\end {array} \right \} =\left [ J\right ] ^{-1}\left \{ \begin {array} [c]{c}\frac {\partial u}{\partial \zeta }\\ \frac {\partial u}{\partial \eta }\end {array} \right \} \tag {5} \end {equation}

and

\begin {equation} \left \{ \begin {array} [c]{c}\frac {\partial v}{\partial x}\\ \frac {\partial v}{\partial y}\end {array} \right \} =\left [ J\right ] ^{-1}\left \{ \begin {array} [c]{c}\frac {\partial v}{\partial \zeta }\\ \frac {\partial v}{\partial \eta }\end {array} \right \} \tag {6} \end {equation}

Now transpose the column vector in (6) to be a row vector because that is how it is laid out in the integral (2A), ( and remember to change the order when transposing a product)

$\left ( \begin {array} [c]{cc}\frac {\partial v}{\partial x} & \frac {\partial v}{\partial y}\end {array} \right ) =\left ( \begin {array} [c]{cc}\frac {\partial v}{\partial \zeta } & \frac {\partial v}{\partial \eta }\end {array} \right ) \left [ J\right ] ^{-T}$

Now we are ready to convert the integral $$I_{2}$$ in eq (2A) to the natural coordinates system (these are area coordinates, notice the integral limits and the order of integration)

\begin {align} I_{2} & ={\displaystyle \int \limits _{\Omega _{j}}} \left ( \begin {array}[c]{cc}\frac {\partial v}{\partial x} & \frac {\partial v}{\partial y}\end {array} \right ) \left \{ \begin {array}[c]{c}\frac {\partial u}{\partial x}\\ \frac {\partial u}{\partial y}\end {array} \right \} dx dy \nonumber \\ & ={\displaystyle \int \limits _{0}^{1}} {\displaystyle \int \limits _{0}^{1-\zeta }} \left ( \begin {array}[c]{cc}\frac {\partial v}{\partial \zeta } & \frac {\partial v}{\partial \eta }\end {array} \right ) \left [ J\right ] ^{-T}\left [ J\right ] ^{-1}\left \{ \begin {array}[c]{c} \frac {\partial u}{\partial \zeta }\\ \frac {\partial u}{\partial \eta }\end {array} \right \} \det \left [ J\right ] d\eta d\zeta \tag {7} \end {align}

Where we used the standard relationship that $dxdy=\det \left [ J\right ] d\eta d\zeta$

Remember to put $$d\eta$$ ﬁrst before $$d\zeta$$ since the inner limit is on $$\zeta .$$

Now that we have everything in the natural area coordinates system, we can do the integration. One small point left, which is to determine the diﬀerentials involved in (7).

For this we now need to decide on the actual form of the trial and test functions and on the mapping between the global and the natural coordinates system. The following explains this part, we will come back to the above integral once we have obtained the diﬀerentials $$\frac {\partial v}{\partial \zeta },\frac {\partial v}{\partial \eta },\frac {\partial u}{\partial \zeta },\frac {\partial u}{\partial \eta }$$ and determined the Jacobian.

The following diagram shows the linear transformation we will use. This is a standard transformation where the natural coordinates are called the area coordinates described more below.

We see from the above diagram that

\begin {align} x & =x_{1}^{j}+\left (x_{2}^{j}-x_{1}^{j}\right ) \zeta +\left (x_{3}^{j}-x_{1}^{j}\right ) \eta \ \tag {8}\\ y & =y_{1}^{j}+\left (y_{2}^{j}-y_{1}^{j}\right ) \zeta +\left (y_{3}^{j}-y_{1}^{j}\right ) \eta \ \nonumber \end {align}

From the above we obtain the following diﬀerentials

\begin {align*} \frac {\partial x}{\partial \zeta } & =\left (x_{2}-x_{1}\right ) \\ \frac {\partial x}{\partial \eta \ } & =\left (x_{3}-x_{1}\right ) \\ \frac {\partial y}{\partial \zeta } & =\left (y_{2}-y_{1}\right ) \\ \frac {\partial y}{\partial \eta \ } & =\left (y_{3}-y_{1}\right ) \end {align*}

Now, we consider the trial and test functions. based on the above transformation shown in eq (8), We see that the linear trial and test functions can also be written in similar transformation

\begin {align} u^{j} & =u_{1}^{j}+\left (u_{2}^{j}-u_{1}^{j}\right ) \zeta +\left ( u_{3}^{j}-u_{1}^{j}\right ) \eta \nonumber \\ v^{j} & =v_{1}^{j}+\left (v_{2}^{j}-v_{1}^{j}\right ) \zeta +\left ( v_{3}^{j}-v_{1}^{j}\right ) \eta \tag {9} \end {align}

Again, immediately, we obtain the following diﬀerentials from the above expressions

\begin {align*} \frac {\partial u}{\partial \zeta } & =\left (u_{2}-u_{1}\right ) \\ \frac {\partial u}{\partial \eta \ } & =\left (u_{3}-u_{1}\right ) \\ \frac {\partial v}{\partial \zeta } & =\left (v_{2}-v_{1}\right ) \\ \frac {\partial v}{\partial \eta \ } & =\left (v_{3}-v_{1}\right ) \end {align*}

Hence the Jacobian can now be evaluated (see eq(3) for reference)

\begin {align} \left [ J\right ] & =\left [ \begin {array} [c]{cc}\frac {\partial x}{\partial \zeta } & \frac {\partial y}{\partial \zeta }\\ \frac {\partial x}{\partial \eta } & \frac {\partial y}{\partial \eta }\end {array} \right ] \nonumber \\ & =\left [ \begin {array} [c]{cc}\left (x_{2}-x_{1}\right ) & \left (y_{2}-y_{1}\right ) \\ \left (x_{3}-x_{1}\right ) & \left (y_{3}-y_{1}\right ) \end {array} \right ] \tag {10} \end {align}

And its inverse is

$\left [ J\right ] ^{-1}=\left [ \begin {array} [c]{cc}\left (y_{3}-y_{1}\right ) & \left (y_{1}-y_{2}\right ) \\ \left (x_{1}-x_{3}\right ) & \left (x_{2}-x_{1}\right ) \end {array} \right ] \frac {1}{\det \left [ J\right ] }$

And

$\left [ J\right ] ^{-T}=\left [ \begin {array} [c]{cc}\left (y_{3}-y_{1}\right ) & \left (x_{1}-x_{3}\right ) \\ \left (y_{1}-y_{2}\right ) & \left (x_{2}-x_{1}\right ) \end {array} \right ] \frac {1}{\det \left [ J\right ] }$

Now that we have all the diﬀerentials needed, we can now go back to the integral in eq (7) and compute it:

\begin {align} I_{2} & ={\displaystyle \int \limits _{0}^{1}} {\displaystyle \int \limits _{0}^{1-\zeta }} \left ( \begin {array} [c]{cc}\frac {\partial v}{\partial \zeta } & \frac {\partial v}{\partial \eta }\end {array} \right ) \left [ J\right ] ^{-T}\left [ J\right ] ^{-1}\left \{ \begin {array} [c]{c}\frac {\partial u}{\partial \zeta }\\ \frac {\partial u}{\partial \eta }\end {array} \right \} \det \left [ J\right ] d\eta d\zeta \nonumber \\ & \nonumber \\ & ={\displaystyle \int \limits _{0}^{1}} {\displaystyle \int \limits _{0}^{1-\zeta }} \left [ \begin {array} [c]{cc}\left (v_{2}-v_{1}\right ) & \left (v_{3}-v_{1}\right ) \end {array} \right ] \left [ J\right ] ^{-T}\left [ J\right ] ^{-1}\left \{ \begin {array} [c]{c}\left (u_{2}-u_{1}\right ) \\ \left (u_{3}-u_{1}\right ) \end {array} \right \} \det \left [ J\right ] d\eta d\zeta \nonumber \\ & \nonumber \\ & ={\displaystyle \int \limits _{0}^{1}} {\displaystyle \int \limits _{0}^{1-\zeta }} \left [ \begin {array} [c]{ccc}v_{1} & v_{2} & v_{3}\end {array} \right ] \left \{ \begin {array} [c]{cc}-1 & -1\\ 1 & 0\\ 0 & 1 \end {array} \right \} \left [ J\right ] ^{-T}\left [ J\right ] ^{-1}\left [ \begin {array} [c]{ccc}-1 & 1 & 0\\ -1 & 0 & 1 \end {array} \right ] \left \{ \begin {array} [c]{c}u_{1}\\ u_{2}\\ u_{3}\end {array} \right \} \det \left [ J\right ] d\eta d\zeta \nonumber \\ & \nonumber \\ & =\left [ \begin {array} [c]{ccc}v_{1} & v_{2} & v_{3}\end {array} \right ] \overset {\text {This is the local stiffness matrix for element j}}{\overbrace {\left ( {\displaystyle \int \limits _{0}^{1}} {\displaystyle \int \limits _{0}^{1-\zeta }} \left \{ \begin {array} [c]{cc}-1 & -1\\ 1 & 0\\ 0 & 1 \end {array} \right \} \left [ J\right ] ^{-T}\left [ J\right ] ^{-1}\left [ \begin {array}[c]{ccc}-1 & 1 & 0\\ -1 & 0 & 1 \end {array} \right ] \det \left [ J\right ] d\eta d\zeta \right ) }}\left \{ \begin {array}[c]{c}u_{1}\\ u_{2}\\ u_{3}\end {array} \right \} \tag {11} \end {align}

Now we can evaluate $$K^{j}$$.

$K^{j}={\displaystyle \int \limits _{0}^{1}} {\displaystyle \int \limits _{0}^{1-\zeta }} \left \{ \begin {array} [c]{cc}-1 & -1\\ 1 & 0\\ 0 & 1 \end {array} \right \} \left [ J\right ] ^{-T}\left [ J\right ] ^{-1}\left [ \begin {array}[c]{ccc}-1 & 1 & 0\\ -1 & 0 & 1 \end {array} \right ] \det \left [ J\right ] d\eta d\zeta$

The integrand is

\begin {align*} I_{3} & =\left \{ \begin {array} [c]{cc}-1 & -1\\ 1 & 0\\ 0 & 1 \end {array} \right \} \left [ J\right ] ^{-T}\left [ J\right ] ^{-1}\left [ \begin {array}[c]{ccc}-1 & 1 & 0\\ -1 & 0 & 1 \end {array} \right ] \det \left [ J\right ] \\ & =\left \{ \begin {array}[c]{cc}-1 & -1\\ 1 & 0\\ 0 & 1 \end {array} \right \} \overset {J^{-T}}{\overbrace {\left [ \begin {array}[c]{cc}\left (y_{3}-y_{1}\right ) & \left (x_{1}-x_{3}\right ) \\ \left (y_{1}-y_{2}\right ) & \left (x_{2}-x_{1}\right ) \end {array} \right ] \frac {1}{\det \left [ J\right ] }}}\overset {J^{-1}}{\overbrace {\left [ \begin {array}[c]{cc}\left (y_{3}-y_{1}\right ) & \left (y_{1}-y_{2}\right ) \\ \left (x_{1}-x_{3}\right ) & \left (x_{2}-x_{1}\right ) \end {array} \right ] \frac {1}{\det \left [ J\right ] }}}\left [ \begin {array}[c]{ccc}-1 & 1 & 0\\ -1 & 0 & 1 \end {array} \right ] \det \left [ J\right ] \\ & =\frac {1}{\det \left [ J\right ] }\left [ \begin {array} [c]{ccc}b_{1}^{2}+c_{1}^{2} & b_{2}b_{1}+c_{2}c_{1} & b_{3}b_{1}+c_{3}c_{1}\\ b_{2}b_{1}+c_{2}c_{1} & b_{2}^{2}+c_{2}^{2} & b_{3}b_{2}+c_{3}c_{2}\\ b_{3}b_{1}+c_{3}c_{1} & b_{3}b_{2}+c_{3}c_{2} & b_{3}^{2}+c_{3}^{2}\end {array} \right ] \end {align*}

Where

$b_{1}=y_{2}-y_{3},\ b_{2}=y_{3}-y_{1},\ b_{3}=y_{1}-y_{2},\ c_{1}=x_{3}-x_{2},\ c_{2}=x_{1}-x_{3},\ c_{3}=x_{2}-x_{1}$

Hence $K^{j}={\displaystyle \int \limits _{0}^{1}} {\displaystyle \int \limits _{0}^{1-\zeta }} \frac {1}{\det \left [ J\right ] }\left [ \begin {array}[c]{ccc}b_{1}^{2}+c_{1}^{2} & b_{2}b_{1}+c_{2}c_{1} & b_{3}b_{1}+c_{3}c_{1}\\ b_{2}b_{1}+c_{2}c_{1} & b_{2}^{2}+c_{2}^{2} & b_{3}b_{2}+c_{3}c_{2}\\ b_{3}b_{1}+c_{3}c_{1} & b_{3}b_{2}+c_{3}c_{2} & b_{3}^{2}+c_{3}^{2}\end {array} \right ] d\eta d\zeta$

But the integrand is a constant, hence we take it out of the integral

$K^{j}=\frac {1}{\det \left [ J\right ] }\left [ \begin {array} [c]{ccc}b_{1}^{2}+c_{1}^{2} & b_{2}b_{1}+c_{2}c_{1} & b_{3}b_{1}+c_{3}c_{1}\\ b_{2}b_{1}+c_{2}c_{1} & b_{2}^{2}+c_{2}^{2} & b_{3}b_{2}+c_{3}c_{2}\\ b_{3}b_{1}+c_{3}c_{1} & b_{3}b_{2}+c_{3}c_{2} & b_{3}^{2}+c_{3}^{2}\end {array} \right ] {\displaystyle \int \limits _{0}^{1}} {\displaystyle \int \limits _{0}^{1-\zeta }} d\eta d\zeta$

Now we evaluate $${\displaystyle \int \limits _{0}^{1}} {\displaystyle \int \limits _{0}^{1-\zeta }} d\eta d\zeta$$

\begin {align*} {\displaystyle \int \limits _{0}^{1}} {\displaystyle \int \limits _{0}^{1-\zeta }} d\eta d\zeta & ={\displaystyle \int \limits _{0}^{1}} \left ( {\displaystyle \int \limits _{0}^{1-\zeta }} d\eta \right ) d\zeta \\ & ={\displaystyle \int \limits _{0}^{1}} \left (\eta \right ) _{0}^{1-\zeta } d\zeta \\ & ={\displaystyle \int \limits _{0}^{1}} 1-\zeta d\zeta \\ & =\left (\zeta -\frac {\zeta ^{2}}{2}\right )_{0}^{1}\\ & =1-\frac {1}{2}\\ & =\frac {1}{2} \end {align*}

Hence

$K^{j}=\frac {1}{2\det \left [ J\right ] }\left [ \begin {array} [c]{ccc}b_{1}^{2}+c_{1}^{2} & b_{2}b_{1}+c_{2}c_{1} & b_{3}b_{1}+c_{3}c_{1}\\ b_{2}b_{1}+c_{2}c_{1} & b_{2}^{2}+c_{2}^{2} & b_{3}b_{2}+c_{3}c_{2}\\ b_{3}b_{1}+c_{3}c_{1} & b_{3}b_{2}+c_{3}c_{2} & b_{3}^{2}+c_{3}^{2}\end {array} \right ]$

But from (10) we see that $\det \left [ J\right ] =x_{3}\left (y_{1}-y_{2}\right ) +x_{1}\left ( y_{2}-y_{3}\right ) +x_{2}\left (y_{3}-y_{1}\right )$

and the area of a triangle with corners at $$\left (x_{1},y_{1}\right ),\left (x_{2},y_{2}\right ) ,\left (x_{3},y_{3}\right )$$ is given by

$A=\frac {1}{2}\det \left ( \begin {array}[c]{ccc}x_{1} & y_{1} & 1\\ x_{2} & y_{2} & 1\\ x_{3} & y_{3} & 1 \end {array} \right ) =\frac {1}{2}\left (x_{1}\left (y_{2}-y_{3}\right ) +x_{2}\left ( y_{3}-y_{1}\right ) +x_{3}\left (y_{1}-y_{2}\right ) \right )$

Hence we get $\det \left (\left [ J\right ] ^{j}\right ) =2A^{j}$

Therefore we can replace $$\det \left (\left [ J\right ] ^{j}\right )$$ by $$2A^{j}$$ everywhere. Rewrite the local stiﬀness matrix in terms of the local element area:

$K^{j}=\frac {1}{4A^{j}}\left [ \begin {array} [c]{ccc}b_{1}^{2}+c_{1}^{2} & b_{2}b_{1}+c_{2}c_{1} & b_{3}b_{1}+c_{3}c_{1}\\ b_{2}b_{1}+c_{2}c_{1} & b_{2}^{2}+c_{2}^{2} & b_{3}b_{2}+c_{3}c_{2}\\ b_{3}b_{1}+c_{3}c_{1} & b_{3}b_{2}+c_{3}c_{2} & b_{3}^{2}+c_{3}^{2}\end {array} \right ]$

Now that we have $$K^{j}$$ we plug it back into eq (11)\begin {align*} I_{2} & =\left [ \begin {array} [c]{ccc}v_{1}^{j} & v_{2}^{j} & v_{3}^{j}\end {array} \right ] K^{j}\left \{ \begin {array} [c]{c}u_{1}^{j}\\ u_{2}^{j}\\ u_{3}^{j}\end {array} \right \} \\ & =\left [ \begin {array} [c]{ccc}v_{1}^{j} & v_{2}^{j} & v_{3}^{j}\end {array} \right ] \frac {1}{4A^{j}}\left [ \begin {array} [c]{ccc}b_{1}^{2}+c_{1}^{2} & b_{2}b_{1}+c_{2}c_{1} & b_{3}b_{1}+c_{3}c_{1}\\ b_{2}b_{1}+c_{2}c_{1} & b_{2}^{2}+c_{2}^{2} & b_{3}b_{2}+c_{3}c_{2}\\ b_{3}b_{1}+c_{3}c_{1} & b_{3}b_{2}+c_{3}c_{2} & b_{3}^{2}+c_{3}^{2}\end {array} \right ] \left \{ \begin {array}[c]{c}u_{1}^{j}\\ u_{2}^{j}\\ u_{3}^{j}\end {array} \right \} \end {align*}

And now that we completed this integral we go back to eq (2) :

$I_{1}=- \overset {\text {This is the integral we just completed above}}{\overbrace {{\displaystyle \int \limits _{\Omega _{j}}} \left ( \begin {array} [c]{cc}\frac {\partial v}{\partial x} & \frac {\partial v}{\partial y}\end {array} \right ) \left \{ \begin {array} [c]{c}\frac {\partial u}{\partial x}\\ \frac {\partial u}{\partial y}\end {array} \right \} dA }}+{\displaystyle \int \limits _{\Omega _{j}}}fv dA$

We need to work on the second integral above $${\displaystyle \int \limits _{\Omega _{j}}} fv dA$$ and transform it to the natural coordinates.

\begin {align*} I & ={\displaystyle \int \limits _{\Omega _{j}}} f^{j} v^{j} dx dy\\ & =\left ( \begin {array}[c]{ccc}v_{1}^{j} & v_{2}^{j} & v_{3}^{j}\end {array} \right ) {\displaystyle \int \limits _{\Omega _{j}}} f^{j} dx dy\\ & =\left ( \begin {array}[c]{ccc}v_{1}^{j} & v_{2}^{j} & v_{3}^{j}\end {array} \right ) \left \{ \begin {array} [c]{c}Q_{1}^{j}\\ Q_{2}^{j}\\ Q_{3}^{j}\end {array} \right \} \ \end {align*}

Where \begin {align*} Q_{1}^{j} & =\int _{0}^{1}\int _{0}^{1-\zeta }f_{1}^{j}N_{1}\det \left [ J\right ] ^{j} d\eta d\zeta \\ & =\int _{0}^{1}\int _{0}^{1-\zeta }f_{1}^{j}\left (1-\zeta -\eta \right ) \det \left [ J\right ] ^{j} d\eta d\zeta \\ & =f_{1}^{j}\det \left [ J\right ] ^{j}\int _{0}^{1}\left (\int _{0}^{1-\zeta }\left (1-\zeta -\eta \right ) d\eta \right ) d\zeta \\ & =f_{1}^{j}\det \left [ J\right ] ^{j}\int _{0}^{1}\left (\eta -\zeta \eta -\frac {\eta ^{2}}{2}\right ) _{0}^{1-\zeta }d\zeta \\ & =f_{1}^{j}\det \left [ J\right ] ^{j}\int _{0}^{1}\left (\left ( 1-\zeta \right ) -\zeta \left (1-\zeta \right ) -\frac {\left (1-\zeta \right ) ^{2}}{2}\right ) d\zeta \\ & =f_{1}^{j}\det \left [ J\right ] ^{j}\int _{0}^{1}\left (1-\zeta -\zeta +\zeta ^{2}-\frac {\left (1+\zeta ^{2}-2\zeta \right ) }{2}\right ) d\zeta \\ & =f_{1}^{j}\frac {\det \left [ J\right ] ^{j}}{2}\int _{0}^{1}\left ( 2-4\zeta +2\zeta ^{2}-1-\zeta ^{2}+2\zeta \right ) d\zeta \\ & =f_{1}^{j}\frac {\det \left [ J\right ] ^{j}}{2}\int _{0}^{1}\left ( 1-2\zeta +\zeta ^{2}\right ) d\zeta \\ & =f_{1}^{j}\frac {\det \left [ J\right ] ^{j}}{2}\left (\zeta -2\frac {\zeta ^{2}}{2}+\frac {\zeta ^{3}}{3}\right ) _{0}^{1}\\ & =f_{1}^{j}\frac {\det \left [ J\right ] ^{j}}{2}\left (1-1+\frac {1}{3}\right ) \\ & =f_{1}^{j}\det \left [ J\right ] ^{j}\left (\frac {1}{6}\right ) \end {align*}

And

\begin {align*} Q_{2}^{j} & =\int _{0}^{1}\int _{0}^{1-\zeta }f_{2}^{j}N_{2}\det \left [ J\right ] ^{j}d\eta d\zeta \\ & =\det \left [ J\right ] ^{j}\ f_{2}^{j}\int _{0}^{1}\int _{0}^{1-\zeta }\zeta d\eta d\zeta \\ & =\det \left [ J\right ] ^{j}\ f_{2}^{j}\int _{0}^{1}\left (\int _{0}^{1-\zeta }\zeta d\eta \right ) d\zeta \\ & =\det \left [ J\right ] ^{j}\ f_{2}^{j}\int _{0}^{1}\left (\zeta \eta \right ) _{0}^{1-\zeta }d\zeta \\ & =\det \left [ J\right ] ^{j}\ f_{2}^{j}\int _{0}^{1}\left (\zeta \left ( 1-\zeta \right ) \right ) d\zeta \\ & =\det \left [ J\right ] ^{j}\ f_{2}^{j}\int _{0}^{1}\left (\zeta -\zeta ^{2}\right ) d\zeta \\ & =\det \left [ J\right ] ^{j}\ f_{2}^{j}\left (\frac {\zeta ^{2}}{2}-\frac {\zeta ^{3}}{3}\right ) _{0}^{1}\\ & =\det \left [ J\right ] ^{j}\ f_{2}^{j}\left (\frac {1}{2}-\frac {1}{3}\right ) \\ & =\det \left [ J\right ] ^{j}\ f_{2}^{j}\left (\frac {1}{6}\right ) \end {align*}

And

\begin {align*} Q_{3}^{j} & =\int _{0}^{1}\int _{0}^{1-\zeta }f_{3}^{j}N_{3}\det \left [ J\right ] ^{j}d\eta d\zeta \\ & =\det \left [ J\right ] ^{j}\ f_{3}^{j}\int _{0}^{1}\int _{0}^{1-\zeta }\eta d\eta d\zeta \\ & =\det \left [ J\right ] ^{j}\ f_{3}^{j}\int _{0}^{1}\left (\frac {\eta ^{2}}{2}\right ) _{0}^{1-\zeta }d\zeta \\ & =\det \left [ J\right ] ^{j}\ f_{3}^{j}\int _{0}^{1}\frac {\left ( 1-\zeta \right ) ^{2}}{2}d\zeta \\ & =\det \left [ J\right ] ^{j}\ f_{3}^{j}\int _{0}^{1}\frac {\left ( 1-2\zeta +\zeta ^{2}\right ) }{2}d\zeta \\ & =\frac {\det \left [ J\right ] ^{j}}{2}\ f_{3}^{j}\left (\zeta -\frac {2\zeta ^{2}}{2}+\frac {\zeta ^{3}}{3}\right ) _{0}^{1}\\ & =\frac {\det \left [ J\right ] ^{j}}{2}\ f_{3}^{j}\left (1-1+\frac {1}{3}\right ) \\ & =\det \left [ J\right ] ^{j}\ f_{3}^{j}\left (\frac {1}{6}\right ) \end {align*}

Hence

\begin {align*} I_{3} & =\left ( \begin {array} [c]{ccc}v_{1}^{j} & v_{2}^{j} & v_{3}^{j}\end {array} \right ) \left \{ \begin {array} [c]{c}Q_{1}^{j}\\ Q_{2}^{j}\\ Q_{3}^{j}\end {array} \right \} \\ & =\left ( \begin {array} [c]{ccc}v_{1}^{j} & v_{2}^{j} & v_{3}^{j}\end {array} \right ) \left \{ \begin {array} [c]{c}\det \left [ J\right ] ^{j}f_{1}^{j}\left (\frac {1}{6}\right ) \\ \det \left [ J\right ] ^{j}\ f_{2}^{j}\left (\frac {1}{6}\right ) \\ \det \left [ J\right ] ^{j}\ f_{3}^{j}\left (\frac {1}{6}\right ) \end {array} \right \} \\ & =\frac {\det \left [ J\right ] ^{j}}{6}\left ( \begin {array} [c]{ccc}v_{1}^{j} & v_{2}^{j} & v_{3}^{j}\end {array} \right ) \left \{ \begin {array} [c]{c}f_{1}^{j}\\ \ f_{2}^{j}\\ \ f_{3}^{j}\end {array} \right \} \end {align*}

But $$\det \left [ J\right ] ^{j}=2A^{j}$$ Hence the above becomes

$I_{3}=\frac {A^{j}}{3}\left ( \begin {array} [c]{ccc}v_{1}^{j} & v_{2}^{j} & v_{3}^{j}\end {array} \right ) \left \{ \begin {array} [c]{c}f_{1}^{j}\\ \ f_{2}^{j}\\ \ f_{3}^{j}\end {array} \right \}$ Now we have the integral in eq (2) completed. We now have our local equations completed. Here it is. We next need to assemble them.\begin {align*} - \frac {1}{4A^{j}}\left [ \begin {array} [c]{ccc}v_{1}^{j} & v_{2}^{j} & v_{3}^{j}\end {array} \right ] \left [ \begin {array} [c]{ccc}K_{1,1}^{j} & K_{1,2}^{j} & K_{1,3}^{j}\\ K_{2,1}^{j} & K_{2,2}^{j} & K_{2,3}^{j}\\ K_{3,1}^{j} & K_{3,2}^{j} & K_{3,3}^{j}\end {array} \right ] \left \{ \begin {array}[c]{c}u_{1}^{j}\\ u_{2}^{j}\\ u_{3}^{j}\end {array} \right \} +\overset {\text {This is the above integral we just completed}}{\overbrace {\frac {A^{j}}{3}\left ( \begin {array}[c]{ccc}v_{1}^{j} & v_{2}^{j} & v_{3}^{j}\end {array} \right ) \left \{ \begin {array}[c]{c}f_{1}^{j}\\ \ f_{2}^{j}\\ \ f_{3}^{j}\end {array} \right \} }} & =0 \\ -\frac {1}{4A^{j}}\left [ \begin {array} [c]{ccc}K_{1,1}^{j} & K_{1,2}^{j} & K_{1,3}^{j}\\ K_{2,1}^{j} & K_{2,2}^{j} & K_{2,3}^{j}\\ K_{3,1}^{j} & K_{3,2}^{j} & K_{3,3}^{j}\end {array} \right ] \left \{ \begin {array} [c]{c}u_{1}^{j}\\ u_{2}^{j}\\ u_{3}^{j}\end {array} \right \} +\frac {A^{j}}{3}\left \{ \begin {array}[c]{c} f_{1}^{j}\\ \ f_{2}^{j}\\ \ f_{3}^{j}\end {array} \right \} & =0 \\ \frac {1}{4A^{j}}\left [ \begin {array} [c]{ccc}K_{1,1}^{j} & K_{1,2}^{j} & K_{1,3}^{j}\\ K_{2,1}^{j} & K_{2,2}^{j} & K_{2,3}^{j}\\ K_{3,1}^{j} & K_{3,2}^{j} & K_{3,3}^{j}\end {array} \right ] \left \{ \begin {array}[c]{c}u_{1}^{j}\\ u_{2}^{j}\\ u_{3}^{j}\end {array} \right \} & =\frac {1}{3}A^{j}\left \{ \begin {array}[c]{c}f_{1}^{j}\\ \ f_{2}^{j}\\ \ f_{3}^{j}\end {array} \right \} \end {align*}

Where $\left [ \begin {array} [c]{ccc}K_{1,1}^{j} & K_{1,2}^{j} & K_{1,3}^{j}\\ K_{2,1}^{j} & K_{2,2}^{j} & K_{2,3}^{j}\\ K_{3,1}^{j} & K_{3,2}^{j} & K_{3,3}^{j}\end {array} \right ] =\left [ \begin {array} [c]{ccc}b_{1}^{2}+c_{1}^{2} & b_{2}b_{1}+c_{2}c_{1} & b_{3}b_{1}+c_{3}c_{1}\\ b_{2}b_{1}+c_{2}c_{1} & b_{2}^{2}+c_{2}^{2} & b_{3}b_{2}+c_{3}c_{2}\\ b_{3}b_{1}+c_{3}c_{1} & b_{3}b_{2}+c_{3}c_{2} & b_{3}^{2}+c_{3}^{2}\end {array} \right ]$

#### 4.3 Note on the shape functions

Looking at the trial function in eq (9), repeated here$u^{j}=u_{1}^{j}+\left (u_{2}^{j}-u_{1}^{j}\right ) \zeta +\left (u_{3}^{j}-u_{1}^{j}\right ) \eta \$ Hence we see that\begin {align} u^{j} & =u_{1}^{j}+u_{2}^{j}\zeta -u_{1}^{j}\zeta +u_{3}^{j}\eta -u_{1}^{j}\eta \nonumber \\ & =u_{1}^{j}\overset {N_{1}}{\overbrace {(1-\zeta -\eta )}}+u_{2}^{j}\overset {N_{2}}{\overbrace {\zeta }}+u_{3}^{j}\overset {N_{3}}{\overbrace {\eta }} \tag {2} \end {align}

And since we are looking for a trial function to be of the form $$u_{1}\left ( Basis_{1}\right ) +u_{2}\left (Basis_{2}\right ) +u_{3}\left (Basis_{3}\right )$$ we see from the above that the 3 shape or basis functions are the following

\begin {align*} N_{1} & =1-\zeta -\eta \\ N_{2} & =\zeta \\ N_{3} & =\eta \end {align*}

And since we are using the Galerkin method, where the test function uses the same basis functions as the trial function, we can write the test function as

$v^{j}=v_{1}^{j}N_{1}+v_{2}^{j}N_{2}+v_{3}^{j}N_{3}$

### 5 Assembly of the global stiﬀness matrix

The global stiﬀness matrix $$K$$ is always square and symmetric and positive deﬁnite. (At least for structural analysis). Recall a positive deﬁnite matrix $$K$$ is one such that for any nonzero vector $$x$$ we always have $$x^{\ast }Ax>0$$ where $$x^{\ast }$$ is the conjugate of $$x.$$ Properties of positive deﬁnite matrix is that all its eigenvalues are positive, and it has positive determinant, and hence a positive deﬁnite matrix is always invertible.

In addition, the global stiﬀness matrix is banded. This means that all non-zero elements are found along bands close to the main diagonal of the matrix. Within the band itself, some values can be zero.

The width of the band is a function of the numbering of the nodes used. Diﬀerent node numbering can result in smaller band width. We want to have as small a band width as possible to take advantage of some numerical methods that can utilize banded matrices.

Band width can be reduced if we keep the node numbering in each element as close as possible to each others.

Now that we have found the local stiﬀness matrix $$K^{j}$$ for element $$j$$ we can assemble the global stiﬀness matrix as shown in this diagram. The direct stiﬀness construction method is used. This is explained in the following diagram

### 6 Assembly of the global load vector

This follows in similar fashion as above. The 3 elements Load vector $$\left \{ f\right \}$$ for element $$j$$ is added to the entries of the global load vector $$\left \{ F\right \}$$ using the node numbering mapping.

### 7 Modiﬁcation of the ﬁnal global stiﬀness matrix and load vectors and ﬁnal solution

Now we have the following equation

$\left [ K\right ] \left \{ u\right \} =\left \{ F\right \}$

Where $$K$$ is the assembled global stiﬀness matrix and $$\left \{ F\right \}$$ is the assembled global load vector. Before we solve for $$\left \{ u\right \}$$, which is the stress function at all the nodes, we must modify $$K$$ and $$F$$ to take care of the given boundary conditions. I attach below 2 pages from a book which gives a good explanation and small example on this point.

Now that we have the modiﬁed $$K^{\ast }$$ and $$F^{\ast }$$ you can solve for $$\left \{ u\right \}$$ using your favorite linear equations solver.

### 8 Conclusion

The following are the main steps in solving the torsion problem described in this report.

### 9 References

1. Methods of computer modeling in Engineering & the sciences. Volume 1. Satya N. Atluri. Tech Science Press

2. Lecture notes, MAE 207. Spring and Fall 2006. UCI. Instructor: Professor Atluri SN.

3. Applied ﬁnite element analysis. Larry Segerlind.

4. The ﬁnite element method for engineers. Kenneth Huebner.

5. Mathematical methods in the physical sciences. 2nd edition. Mary Boas.

6. Fellow students reports and code from MAE 207 projects: Roy Culver , Paul Nylandres, Q Wang. see other related reports on my MAE 207 class web page

1The FEM for engineers. Kenneth Huebner. 1975