<<

. 57
( 67 .)



>>

‚x ‚x ‚y ‚y


Step 3 Introduce a discretization for both the weighting function w and the
unknown u, so within each element
n
wh | = Ni ( x, y) we,i = N T ( x, y) we (16.17)
∼ ∼
e
i=1

n
uh | = Ni ( x, y) ue,i = N T ( x, y) ue . (16.18)


e
i=1
The three-dimensional convection-diffusion equation
282

Step 4 Substitution of this discretization into Eq. (16.16) yields
‚wh ‚uh ‚wh ‚uh
∇wh · c∇uh = c +
‚x ‚x ‚y ‚y
‚N T we ‚N T ue ‚N T w ‚N T u
+ ∼ ∼ e ∼ ∼e
=c ∼∼
∼∼
‚x ‚x ‚y ‚y
T ‚N ‚N ‚N ‚N T
T
= c we +∼∼
∼ ∼
u. (16.19)
‚y ‚y ∼e
‚x ‚x





Step 5 Using the discretization in Eq. (16.15) yields
Nel
‚N ‚N T ‚N ‚N T
+∼∼
wT ∼ ∼
c d ue
∼e
‚x ‚x ‚y ‚y ∼
e
e=1
Nel
= + wT N n · c∇u d
wT Nf d . (16.20)
∼e ∼e
∼ ∼
e e
e=1

The element matrix is given by
‚N ‚N T ‚N ‚N T
Ke = +∼∼
∼ ∼
c d, (16.21)
‚x ‚x ‚y ‚y
e

and the element column by

f= + N n · c∇u d .
Nf d (16.22)
∼ ∼
∼e
e e

Using this notation, Eq. (16.20) may be written as

Nel Nel
ue =
wT K e wT f e . (16.23)
∼e ∼e
∼ ∼
e=1 e=1

Following a similar procedure as outlined in Chapter 14, this may be rearranged
into
w T K u = wT f , (16.24)

∼ ∼ ∼

which should hold for all w. This ¬nally leads to:


Ku = f . (16.25)
∼ ∼

Similar to the situation described in Section 14.6 the column u contains an

unknown and a known part depending on the essential boundary conditions. In the
nodes, where essential boundary conditions are prescribed, the associated external
loads in the right-hand side integrals are unknown. However, the set equations can
be partitioned as has been done in Chapter 14 to arrive at a set that can be solved.
283 16.6 Convection-diffusion equation

16.6 Convection-diffusion equation

Assuming isotropic diffusion, the convection-diffusion equation is given by
‚u
+ v · ∇u = ∇ · ( c∇u) + f , (16.26)
‚t
with v the convective velocity. This equation should hold on the spatial domain
during a certain period of time, say S = [ 0, T]. Initial boundary conditions must
be speci¬ed:
u( x, t = 0) = uini ( x) in , (16.27)
as well as essential and natural boundary conditions:
u = U at (16.28)
u


n · c∇u = P at p. (16.29)
The weak form is obtained analogously to the procedure of Section 16.4, giving
‚u
+ wv · ∇u + ∇w · ( c∇u)
w d
‚t
= + wn · c∇u d .
wf d (16.30)

Spatial discretization is performed in a two-dimensional con¬guration by intro-
ducing
wh | = N T ( x, y) we , (16.31)
∼ ∼
e

and
uh | = N T ( x, y) ue . (16.32)


e

For a particular element e, the individual integrals of Eq. (16.30) can be
converted to:
‚u due
= wT M e ∼
w d , (16.33)
∼e
‚t dt
e

with:
Me = NNT d . (16.34)
∼∼
e

Further, by using
v = v x e x + vy e y , (16.35)
it can be written:

wv · ∇u d = wT C e u e , (16.36)
∼e ∼
e
The three-dimensional convection-diffusion equation
284

with
‚N T ‚N T
Ce = + vy ∼

N vx d. (16.37)
‚x ‚y

e

The remaining integrals from Eq. (16.30) follow from the previous section.
Therefore, Eq. (16.30) may be formulated as
Nel Nel
du
M e ∼e + ( Ce + K e ) ue =
wT wT f e , (16.38)
∼e ∼e
∼ ∼
dt
e=1 e=1
which must be satis¬ed for all admissible weighting values, hence after the
assembly process the following set of equations results:
du
+ ( C + K) u = f .

M (16.39)
∼ ∼
dt
Application of the θ-scheme for temporal discretization in the time increment
[ tn , tn+1 ] results in:
1
+ θC + θK un+1
M ∼
t
1
un + f θ ,
=M ’( 1 ’ θ) C’( 1 ’ θ) K (16.40)
∼ ∼
t
t = tn+1 ’ tn and
with
f θ = θf n+1 + ( 1 ’ θ) f n . (16.41)
∼ ∼ ∼




16.7 Isoparametric elements and numerical integration

To compute the element matrices M e , Ce and K e and the element array f e the shape

functions N and the shape function derivatives with respect to the coordinates

x and y, need to be available. To de¬ne the shape functions it is convenient to
introduce a local coordinate system.
In particular for an arbitrarily shaped quadrilateral element, for instance as
depicted in Fig. 16.4(a) it is dif¬cult or impossible to de¬ne explicitly the shape
functions with respect to the global coordinates x and y. However, de¬ning the
shape functions in a unit square domain, as depicted in Fig. 16.4(b) is straight-
forward. The local coordinate system de¬ned in Fig. 16.4(b) is chosen such that
’1 ¤ ξ ¤ 1 and ’1 ¤ · ¤ 1. Node one, for instance, has local coordinates
ξ = ’1 and · = ’1. Remember that the shape function of a certain node has
to be one at the spatial location of this node and zero at the spatial location of all
other nodes. Using the local coordinate system allows an elegant de¬nition of the
shape functions of a four-node quadrilateral element:
285 16.7 Isoparametric elements and numerical integration

·
3 (“1, 1) (1, 1)
4 4 3

ξ
ey 1
2

2
1
ex (1, “1)
(“1, “1)
(a) Quadrilateral element with (b) Quadrilateral element with
respect to a global coordinate system respect to a local coordinate system

Figure 16.4
Quadrilateral element with respect to global and local coordinate systems.


1
N1 = (1 ’ ξ ) ( 1 ’ ·)
4
1
N2 = (1 + ξ ) ( 1 ’ ·)
4 (16.42)
1
N3 = (1 + ξ ) ( 1 + ·)
4
1
N4 = (1 ’ ξ ) ( 1 + ·) .
4
An element having these shape functions is called a bi-linear element. Along the
edges of the element the shape functions are linear with respect to either ξ or ·.
Within the element, however, the shape functions are bi-linear with respect to ξ
and ·. For instance:
1
N1 =
(1 ’ ξ ’ · + ξ ·) . (16.43)
4
Fig. 16.5 shows N1 visualized as a contour plot.

<<

. 57
( 67 .)



>>