4.15 Solve \(\nabla ^{2}u=f\left ( x,y\right ) \) on square using FEM. Assume \(u=0\) on boundaries and solve using \(f\left ( x,y\right ) =xy\) and also \(f\left ( x,y\right ) =1\) and
\(f\left ( x,y\right ) =3\left ( \cos (4 \pi x)+\sin (3 \pi y)\right )\)
Use Galerkin method and weak form, Using a bilinear trial function. Let width of square
be 1.
Solution
Using as an example with 9 elements to illustrate the method. The program below can be
called with different number of elements.
The trial function is bilinear
\begin{equation} \tilde {u}=c_{1}+c_{2}x+c_{3}y+c_{4}xy \tag {1}\end{equation}
Looking at one element, and using local coordinates systems with element having width \(2a\)
and height \(2b\) gives
Evaluating the trial function (1) at each corner node of the above element gives
\begin{align*} \tilde {u}_{1} & =c_{1}-ac_{2}-bc_{3}+abc_{4}\\ \tilde {u}_{2} & =c_{1}+ac_{2}-bc_{3}-abc_{4}\\ \tilde {u}_{3} & =c_{1}+ac_{2}+by_{3}+abc_{4}\\ \tilde {u}_{4} & =c_{1}-ac_{2}+bc_{3}+abc_{4}\end{align*}
Or
\[\begin {Bmatrix} \tilde {u}_{1}\\ \tilde {u}_{2}\\ \tilde {u}_{3}\\ \tilde {u}_{4}\end {Bmatrix} =\begin {pmatrix} 1 & -a & -b & ab\\ 1 & a & -b & -ab\\ 1 & a & b & ab\\ 1 & -a & b & ab \end {pmatrix}\begin {Bmatrix} c_{1}\\ c_{2}\\ c_{3}\\ c_{4}\end {Bmatrix} \]
Hence
\begin{align}\begin {Bmatrix} c_{1}\\ c_{2}\\ c_{3}\\ c_{4}\end {Bmatrix} & =\begin {pmatrix} 1 & -a & -b & ab\\ 1 & a & -b & -ab\\ 1 & a & b & ab\\ 1 & -a & b & ab \end {pmatrix} ^{-1}\begin {Bmatrix} \tilde {u}_{1}\\ \tilde {u}_{2}\\ \tilde {u}_{3}\\ \tilde {u}_{4}\end {Bmatrix} \nonumber \\ & =\frac {1}{4ab}\begin {pmatrix} ab & ab & ab & ab\\ -b & b & b & -b\\ -a & -a & a & a\\ 1 & -1 & 1 & -1 \end {pmatrix}\begin {Bmatrix} \tilde {u}_{1}\\ \tilde {u}_{2}\\ \tilde {u}_{3}\\ \tilde {u}_{4}\end {Bmatrix} \tag {2}\end{align}
coor = {{-a, -b}, {a, -b}, {a, b}, {-a, b}};
uVal = {u1, u2, u3, u4}; eq = c1 + c2 x + c3 y + c4 x y;
r = MapThread [eq /. {u -> #1, x -> #2[[1]], y -> #2[[2]]}&, {uVal, coor}];
mat = Last@Normal[ CoefficientArrays [r, {c1, c2, c3, c4}]];
MatrixForm [ Inverse [ mat ] // MatrixForm
Substituting (2) back into (1) and rearranging terms results in
\begin{gather*} \begin {aligned} \tilde {u} & =c_{1}+c_{2}x+c_{3}y+c_{4}xy\\ & =\frac {1}{4ab}\left ( \left ( ab-bx-ay+xy\right ) \tilde {u}_{1}+\left ( ab+bx-ay-xy\right ) \tilde {u}_{2}+\left ( ab+bx+ay+xy\right ) +\left ( ab-bx+ay-xy\right ) \right ) \end {aligned} \end{gather*}
Since \(ab\) is the \(\frac {1}{4}\) of the area of element, then the above becomes
\begin{gather*} \tilde {u}=\frac {1}{A}\left ( \left ( ab-bx-ay+xy\right ) \tilde {u}_{1}+\left ( ab+bx-ay-xy\right ) \tilde {u}_{2}+\left ( ab+bx+ay+xy\right ) +\left ( ab-bx+ay-xy\right ) \right ) \end{gather*}
The above can now be written in term of what is called shape functions\begin{gather*} \tilde {u}(x,y)=N_{1}(x,y)\tilde {u_{1}}+N_{2}(x,y)\tilde {u_{2}}N_{3}(x,y)\tilde {u_{3}}+N_{4}(x,y)\tilde {u_{4}}\end{gather*} Where\begin{align*} N_{1} & =\frac {1}{A}\left ( ab-bx-ay+xy\right ) =\frac {1}{A}\left ( a-x\right ) \left ( b-y\right ) \\ N_{2} & =\frac {1}{A}\left ( ab+bx-ay-xy\right ) =\frac {1}{A}\left ( a+x\right ) \left ( b-y\right ) \\ N_{3} & =\frac {1}{A}\left ( ab+bx+ay+xy\right ) =\frac {1}{A}\left ( a+x\right ) \left ( b+y\right ) \\ N_{4} & =\frac {1}{A}\left ( ab-bx+ay-xy\right ) =\frac {1}{A}\left ( a-x\right ) \left ( b+y\right ) \end{align*}
Now that the shape functions are found, the next step is to determine the local element
stiffness matrix. This can be found from the weak form integral over the area of the
element.
\begin{equation} I_{i}=\int \limits _{x=-a}^{x=a}\int \limits _{y=-b}^{y=b}w\left ( \nabla ^{2}u-f\left ( x,y\right ) \right ) \,dxdy \tag {3}\end{equation}
Where \(w\left ( x,y\right ) \) is the test function. For Galerkin method, the test function \(w_{i}=\frac {d\tilde {u}}{du_{i}}=N_{i}\left ( x,y\right ) \) . Hence
\begin{gather*} w\left ( x,y\right ) =\begin {Bmatrix} N_{1}\left ( x,y\right ) \\ N_{2}\left ( x,y\right ) \\ N_{3}\left ( x,y\right ) \\ N_{4}\left ( x,y\right ) \end {Bmatrix} =\begin {Bmatrix} \frac {1}{A}\left ( a-x\right ) \left ( b-y\right ) \\ \frac {1}{A}\left ( a+x\right ) \left ( b-y\right ) \\ \frac {1}{A}\left ( a+x\right ) \left ( b+y\right ) \\ \frac {1}{A}\left ( a-x\right ) \left ( b+y\right ) \end {Bmatrix} \end{gather*}
Using weak form, integration by parts is applied to (2)
\begin{gather*} I_{i}=\int \limits _{x=-a}^{x=a}\int \limits _{y=-b}^{y=b}\left ( -\frac {\partial w}{\partial x_{i}}\frac {\partial \tilde {u}}{\partial x_{i}}-wf\left ( x,y\right ) \right ) \,dxdy+\overbrace {\int \limits _{\Gamma }w\frac {\partial \tilde {u}}{\partial n}}^{\text {goes to zero}} \end{gather*}
The term \(\int \limits _{\Gamma }w\frac {\partial \tilde {u}}{\partial n}\) is the integration over the boundary of the element. Since there is only an
essential boundary condition over all the boundaries (this is the given Dirchilet boundary
condition), \(w=0\) on the boundary and this integral vanishes. There is no natural boundary
conditions for this example.
For those elements not on the external edge of the overall grid (i.e. internal elements), each
contribution to this integral will cancel from the adjacent internal element. What this
means is that the above reduces to just the first integral
\begin{align*} I_{i} & =\int \limits _{x=-a}^{x=a}\int \limits _{y=-b}^{y=b}\left ( -\frac {\partial w}{\partial x_{i}}\frac {\partial \tilde {u}}{\partial x_{i}}-wf\left ( x,y\right ) \right ) \,dxdy\\ & =\int \limits _{x=-a}^{x=a}\int \limits _{y=-b}^{y=b}-\begin {Bmatrix} \frac {\partial N_{1}}{\partial x}\\ \frac {\partial N_{2}}{\partial x}\\ \frac {\partial N_{3}}{\partial x}\\ \frac {\partial N_{4}}{\partial x}\end {Bmatrix}\begin {Bmatrix} \frac {\partial N_{1}}{\partial x} & \frac {\partial N_{2}}{\partial x} & \frac {\partial N_{3}}{\partial x} & \frac {\partial N_{4}}{\partial x}\end {Bmatrix}\begin {Bmatrix} \tilde {u}_{1}\\ \tilde {u}_{2}\\ \tilde {u}_{3}\\ \tilde {u}_{4}\end {Bmatrix} -\begin {Bmatrix} \frac {\partial N_{1}}{\partial y}\\ \frac {\partial N_{2}}{\partial y}\\ \frac {\partial N_{3}}{\partial y}\\ \frac {\partial N_{4}}{\partial y}\end {Bmatrix}\begin {Bmatrix} \frac {\partial N_{1}}{\partial y} & \frac {\partial N_{2}}{\partial y} & \frac {\partial N_{3}}{\partial y} & \frac {\partial N_{4}}{\partial y}\end {Bmatrix}\begin {Bmatrix} \tilde {u}_{1}\\ \tilde {u}_{2}\\ \tilde {u}_{3}\\ \tilde {u}_{4}\end {Bmatrix} \\ & -\int \limits _{x=-a}^{x=a}\int \limits _{y=-b}^{y=b}\begin {Bmatrix} N_{1}\left ( x,y\right ) \\ N_{2}\left ( x,y\right ) \\ N_{3}\left ( x,y\right ) \\ N_{4}\left ( x,y\right ) \end {Bmatrix} f\left ( x,y\right ) \,dxdy \end{align*}
Hence
\begin{gather*} \begin {aligned} I_{i} & =\int \limits _{x=-a}^{x=a}\int \limits _{y=-b}^{y=b}-\begin {pmatrix} \frac {\partial N_{1}}{\partial x}\frac {\partial N_{1}}{\partial x}+\frac {\partial N_{1}}{\partial y}\frac {\partial N_{1}}{\partial y} & \frac {\partial N_{1}}{\partial x}\frac {\partial N_{2}}{\partial x}+\frac {\partial N_{1}}{\partial y}\frac {\partial N_{2}}{\partial y} & \frac {\partial N_{1}}{\partial x}\frac {\partial N_{3}}{\partial x}+\frac {\partial N_{1}}{\partial y}\frac {\partial N_{3}}{\partial y} & \frac {\partial N_{1}}{\partial x}\frac {\partial N_{4}}{\partial x}+\frac {\partial N_{1}}{\partial y}\frac {\partial N_{4}}{\partial y}\\ \frac {\partial N_{2}}{\partial x}\frac {\partial N_{1}}{\partial x}+\frac {\partial N_{2}}{\partial y}\frac {\partial N_{1}}{\partial y} & \frac {\partial N_{2}}{\partial x}\frac {\partial N_{2}}{\partial x}+\frac {\partial N_{2}}{\partial y}\frac {\partial N_{2}}{\partial y} & \frac {\partial N_{2}}{\partial x}\frac {\partial N_{3}}{\partial x}+\frac {\partial N_{2}}{\partial y}\frac {\partial N_{3}}{\partial y} & \frac {\partial N_{2}}{\partial x}\frac {\partial N_{4}}{\partial x}+\frac {\partial N_{2}}{\partial y}\frac {\partial N_{4}}{\partial y}\\ \frac {\partial N_{3}}{\partial x}\frac {\partial N_{1}}{\partial x}+\frac {\partial N_{3}}{\partial y}\frac {\partial N_{1}}{\partial y} & \frac {\partial N_{3}}{\partial x}\frac {\partial N_{2}}{\partial x}+\frac {\partial N_{3}}{\partial y}\frac {\partial N_{2}}{\partial y} & \frac {\partial N_{3}}{\partial x}\frac {\partial N_{3}}{\partial x}+\frac {\partial N_{3}}{\partial y}\frac {\partial N_{3}}{\partial y} & \frac {\partial N_{3}}{\partial x}\frac {\partial N_{4}}{\partial x}+\frac {\partial N_{3}}{\partial y}\frac {\partial N_{4}}{\partial y}\\ \frac {\partial N_{4}}{\partial x}\frac {\partial N_{1}}{\partial x}+\frac {\partial N_{4}}{\partial y}\frac {\partial N_{1}}{\partial y} & \frac {\partial N_{4}}{\partial x}\frac {\partial N_{2}}{\partial x}+\frac {\partial N_{4}}{\partial y}\frac {\partial N_{2}}{\partial y} & \frac {\partial N_{4}}{\partial x}\frac {\partial N_{3}}{\partial x}+\frac {\partial N_{4}}{\partial y}\frac {\partial N_{3}}{\partial y} & \frac {\partial N_{4}}{\partial x}\frac {\partial N_{4}}{\partial x}+\frac {\partial N_{4}}{\partial y}\frac {\partial N_{4}}{\partial y}\end {pmatrix}\begin {Bmatrix} \tilde {u}_{1}\\ \tilde {u}_{2}\\ \tilde {u}_{3}\\ \tilde {u}_{4}\end {Bmatrix} \\ & -\int \limits _{x=-a}^{x=a}\int \limits _{y=-b}^{y=b}\begin {Bmatrix} N_{1}\left ( x,y\right ) \\ N_{2}\left ( x,y\right ) \\ N_{3}\left ( x,y\right ) \\ N_{4}\left ( x,y\right ) \end {Bmatrix} f\left ( x,y\right ) \,dxdy \end {aligned} \end{gather*}
In the above, we have the from \(I_{i}=\int \limits _{\Omega }-k_{i}\left \{ \tilde {u}\right \} d\Omega -\int \limits _{\Omega }\left \{ \tilde {u}\right \} f_{i}d\Omega \) , hence the element stiffness matrix is the first integral,
and the force vector comes from the second integral.
The integration is now carried out to obtain the element stiffness matrix. This was done
using Mathematica. The local stiffness matrix for element \(i\) is
\begin{align*} k_{i} & =\int \limits _{x=-a}^{x=a}\int \limits _{y=-b}^{y=b}\frac {\partial w}{\partial x_{i}}\frac {\partial \tilde {u}}{\partial x_{i}}dxdy\\ & =\begin {pmatrix} \frac {a^{2}+b^{2}}{3ab} & \frac {a}{6b}-\frac {b}{3a} & -\frac {a^{2}+b^{2}}{6ab} & -\frac {a}{3b}+\frac {b}{6a}\\ \frac {a}{6b}-\frac {b}{3a} & \frac {a^{2}+b^{2}}{3ab} & -\frac {a}{3b}+\frac {b}{6a} & -\frac {a^{2}+b^{2}}{6ab}\\ -\frac {a^{2}+b^{2}}{6ab} & -\frac {a}{3b}+\frac {b}{6a} & \frac {a^{2}+b^{2}}{3ab} & \frac {a}{6b}-\frac {b}{3a}\\ -\frac {a}{3b}+\frac {b}{6a} & -\frac {a^{2}+b^{2}}{6ab} & \frac {a}{6b}-\frac {b}{3a} & \frac {a^{2}+b^{2}}{3ab}\end {pmatrix} \end{align*}
For example, for element of width 1 and height 1, then \(a=\frac {1}{2},b=\frac {1}{2}\) and the above becomes
\[ k_{i}= \begin {pmatrix} \frac {2}{3} & -\frac {1}{6} & -\frac {1}{3} & -\frac {1}{6}\\ -\frac {1}{6} & \frac {2}{3} & -\frac {1}{6} & -\frac {1}{3}\\ -\frac {1}{3} & -\frac {1}{6} & \frac {2}{3} & -\frac {1}{6}\\ -\frac {1}{6} & -\frac {1}{3} & -\frac {1}{6} & \frac {2}{3}\end {pmatrix} \]
ClearAll [ a , b];
area = 4 a b;
phi = (1/area) {(a - x) (b - y),
(a + x) (b - y),
(a + x) (b + y),
(a - x) (b + y)};
vx = {{ D [phi[[1]], x]},
{ D [phi[[2]], x]},
{ D [phi[[3]], x]},
{ D [phi[[4]], x]}};
vy = {{ D [phi[[1]], y]},
{ D [phi[[2]], y]},
{ D [phi[[3]], y]},
{ D [phi[[4]], y]}};
k = Integrate [vx. Transpose [vx]+vy. Transpose [vy],{x, -a, a}, {y, -b, b}]
Hence
\[ I_{i}=-k_{i}\begin {Bmatrix} \tilde {u}_{1}\\ \tilde {u}_{2}\\ \tilde {u}_{3}\\ \tilde {u}_{4}\end {Bmatrix} -\int \limits _{x=-a}^{x=a}\int \limits _{y=-b}^{y=b}\begin {Bmatrix} N_{1}\left ( x,y\right ) \\ N_{2}\left ( x,y\right ) \\ N_{3}\left ( x,y\right ) \\ N_{4}\left ( x,y\right ) \end {Bmatrix} f\left ( x,y\right ) \,dxdy \]
Now the integration of the force vector is carried out.
\begin{align*} I_{i}^{f} & =\int \limits _{x=-a}^{x=a}\int \limits _{y=-b}^{y=b}\begin {Bmatrix} N_{1}\left ( x,y\right ) \\ N_{2}\left ( x,y\right ) \\ N_{3}\left ( x,y\right ) \\ N_{4}\left ( x,y\right ) \end {Bmatrix} f\left ( x,y\right ) \,dxdy\\ & =\int \limits _{x=-a}^{x=a}\int \limits _{y=-b}^{y=b}\begin {Bmatrix} \frac {1}{A}\left ( a-x\right ) \left ( b-y\right ) \\ \frac {1}{A}\left ( a+x\right ) \left ( b-y\right ) \\ \frac {1}{A}\left ( a+x\right ) \left ( b+y\right ) \\ \frac {1}{A}\left ( a-x\right ) \left ( b+y\right ) \end {Bmatrix} f\left ( x,y\right ) \,dxdy \end{align*}
In the above, the integration depends on the physical location of the element. We used a
local coordinate system initially to obtain the shape functions. This has now to be
transformed to the global coordinates system. Taking the center of each element as \(\left ( x_{0},y_{0}\right ) \) then
we need to replace \(a\rightarrow x_{0}+a\) and \(b\rightarrow y_{0}+b\) everywhere in the above integral. We did not have to do this for
finding the local stiffness matrix, since that did not depend on the physical location of the
element (it was constant). But for the integration of the force function, we need to do this
mapping. In the code, the center of each element is found, and the replacement is
done.