<<

. 51
( 67 .)



>>


N T ( ξi ) = n(i,:).



Using these shape functions the coordinate x( ξi ) within element e can be computed at
the i-th integration point:

x( ξi ) = N T ( ξi ) ∼e = n(i,:)*nodcoord.
x



Similarly, the value of the solution at the i-th integration point of element e is
obtained via

u( ξi ) = N T ( ξi ) ue = n(i,:)*nodu.




In a similar fashion the shape function derivatives with respect to the local coordinate
ξ are stored in an array, called dndxi.


Structure of the Finite Element code Typically the structure of a ¬nite element
programme is as follows.
(i) Pre-processing: mesh generation, boundary condition speci¬cation and parame-
ter declaration. This should provide the topology array top, the coordinate array
coord and a number of auxiliary arrays containing boundary conditions and
material parameters.
(ii) Based on the mesh and element types used, the index array pos can be computed.
(iii) Assembly of the coef¬cient matrix q = K and the element array rhs = f . Let

qe = K e and rhse = f e , then the assembly process in a MATLAB environment

would look like

% nelem: the number of elements

for ielem = 1:nelem

% compute qe and rhse

[qe,rhse]=<elementfunction>(ielem,coord,top,.....)

% get the location of the degrees of freedom
% in the solution array

ii = nonzeros(pos(ielem,:));

% add the element coefficient matrix
% and the element right-hand side array
253 14.9 Example

% to the total coefficient matrix q
% and the load array rhs
q(ii,ii) = q(ii,ii) + qe;
rhs(ii) = rhs(ii) + rhse;

end
(iv) Solution of the set of equations taking into account the boundary conditions.
(v) Post-processing based on the solution, for instance by computing associated quanti-
ties such as heat ¬‚uxes or stresses.




14.9 Example

As an example consider the diffusion problem with the following parameter set-
ting. We consider the domain : 0 ¤ x ¤ 1, with prescribed essential boundary
conditions at x = 0 and x = 1. These conditions are: u(0) = 0 and u(1) = 0. There
are no natural boundary conditions. The material constant satis¬es: c = 1 and the
source term: f = 1.
The domain is divided into ¬ve elements of equal length. Fig. 14.9 shows
the solution. The left part displays the computed solution u (solid line) as well as
the exact solution (dashed line). Remarkably, in this one-dimensional case with


0.14 0.5

0.4
0.12
0.3
0.1 0.2

0.1
0.08
u 0
p
0.06
“0.1

“0.2
0.04
“0.3
0.02
“0.4

0 “0.5
0.2 0.4 0.6 0.8 1
0 0 0.2 0.4 0.6 0.8 1
x x

Figure 14.9
Five element solution. Left: (solid line) approximate solution uh ( x), (dashed line) exact solution u( x).
Right: (solid line) approximate ¬‚ux ph ( x), (dashed line) exact ¬‚ux p( x).
Numerical solution of one-dimensional diffusion equation
254

the current choice of parameters, the nodal solutions are exact. The right part of
the ¬gure shows the computed ¬‚ux, say ¬‚ux p = c du/dx. Again, the solid line
denotes the computed ¬‚ux p, which is clearly discontinuous from one element
to the next, and the dashed line denotes the exact solution. The discontinuity of
the computed ¬‚ux ¬eld is obvious: the ¬eld u is piecewise linear, therefore the
derivative du/dx is piecewise constant. The ¬‚ux p does not necessarily have to be
piecewise constant: if the parameter c is a function of x, the ¬‚ux p will be varying
within an element.
Mesh re¬nement leads to an improved approximate solution. For instance,
using ten rather than ¬ve elements, yields the results depicted in Fig. 14.10. The
impact of a varying c, say c = 1 + x is depicted in Fig. 14.11.
Changing the interpolation order of the shape functions Ni from linear to
quadratic, also has a signi¬cant impact on the results, in particular on the qual-
ity of the ¬‚ux prediction. For the constant c case, the solution becomes exact.
Also for c = 1 + x a signi¬cant improvement can be observed, as depicted in
Fig. 14.13.




0.5
0.14

0.4
0.12
0.3
0.1 0.2

0.1
0.08
p 0
u
0.06
“0.1

“0.2
0.04
“0.3
0.02
“0.4

“0.5
0
0 0.2 0.4 0.6 0.8 1
0 0.2 0.4 0.6 0.8 1
x x

Figure 14.10
Ten element solution. Left: (solid line) approximate solution uh ( x), (dashed line) exact solution u( x).
Right: (solid line) approximate ¬‚ux ph ( x), (dashed line) exact ¬‚ux p( x).
255 14.9 Example

0.6
0.09

0.08
0.4
0.07
0.2
0.06
0
0.05
u p
0.04 “0.2
0.03
“0.4
0.02
“0.6
0.01

0 “0.8
0 0.2 0.4 0.6 0.8 1
0 0.2 0.4 0.6 0.8 1
x x

Figure 14.11
Five element solution for c = 1 + x. Left: (solid line) approximate solution uh ( x), (dashed line) exact
solution u( x). Right: (solid line) approximate ¬‚ux ph ( x), (dashed line) exact ¬‚ux p( x).



0.5
0.14

0.4
0.12
0.3
0.1 0.2

0.1
0.08
0
p
u
0.06
“0.1

“0.2
0.04
“0.3
0.02
“0.4

“0.5
0
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
x x

Figure 14.12
Solution for c = 1 using ¬ve quadratic elements. Left: (solid line) approximate solution uh ( x),
(dashed line) exact solution u( x). Right: (solid line) approximate ¬‚ux ph ( x), (dashed line) exact ¬‚ux
p( x).
Numerical solution of one-dimensional diffusion equation
256

0.6
0.09

0.08
0.4
0.07
0.2
0.06

0
0.05
p
u
0.04 “0.2

0.03
“0.4
0.02
“0.6
0.01

“0.8
0
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
x x

Figure 14.13
Solution for c = 1 + x using ¬ve quadratic elements. Left: (solid line) approximate solution uh ( x),
(dashed line) exact solution u( x). Right: (solid line) approximate ¬‚ux ph ( x), (dashed line) exact ¬‚ux
p( x).



Exercises

14.1 The method of weighted residuals can be used to ¬nd approximations of a
given function. Let f ( x) be a function that one would like to approximate
with a polynomial of the order n in a certain domain, say 0 ¤ x ¤ 1. Let
the polynomial be given by
h( x) = a0 + a1 x + · · · + an xn .
Ideally
g( x) = h( x) ’f ( x) = 0 for all x with 0 ¤ x ¤ 1.
In terms of the weighted residuals equation this may also be expressed as
1

<<

. 51
( 67 .)



>>