Computational Engineering

REPORT 2
-Juan Pablo Wilches Arboleda-
0. Problem statement

Along this report it is pretended to develop the following transient convection-diffusion equation according
to the parameters given. For instance, the problem to be model is the evolutions in time of the
concentration of a pollutant in a river. The following transient 1D equation is considered first:

𝜕𝑡 𝑢 + 𝑎𝑢𝑥 − 𝑣𝑢𝑥𝑥 = 𝑠 𝑥 ∈ (0,1) (1)

With source term 𝑠, convection velocity 𝑎 and diffusivity 𝑣. The problem is complemented with natural
boundary conditions, 𝑣𝑢𝑥 = 0, on both ends. Consider an initial condition 𝑢(𝑥, 0) = 𝑔(𝑥). The problem is
discretized with a uniform mesh of 100 elements and the final time is 5.

The following parameters are considered:

𝑎 = 1, 𝑠 = 1, 𝑣 = 0,01

1. Weak form of the problem

Aimed to obtain the weak form of the problem, the transient 1D equation (1) should be multiplied by the
test function 𝑤, which becomes 0 at both ends. Then, after multiplying by test function it proceeds to
integrate over all the domain (between 0 and 1).
1 1 1 1
∫ 𝑤𝑢𝑡 𝑑𝑥 + ∫ 𝑤𝑎𝑢𝑥 𝑑𝑥 − ∫ 𝑤𝑣𝑢𝑥𝑥 𝑑𝑥 = ∫ 𝑤 𝑠 𝑑𝑥 (2)
0 0 0 0

1
As it is known the boundary conditions on both ends of the equation (1), then the term ∫0 𝑤𝑣𝑢𝑥𝑥 𝑑𝑥 must
be integrated by parts in order to apply these conditions. The integration by parts follows a procedure such
as:

∫ 𝑢 ∙ 𝑑𝑣 = 𝑢 ∙ 𝑣 − ∫ 𝑣 ∙ 𝑑𝑢

𝑢=𝑤 ; 𝑑𝑣 = 𝑣𝑢𝑥𝑥 𝑑𝑥 ; 𝑑𝑢 = 𝑤𝑥 𝑑𝑥 ; 𝑣 = 𝑣𝑢𝑥

1 0 1
∫ 𝑤𝑣𝑢𝑥𝑥 𝑑𝑥 = 𝑤 ∙ 𝑣𝑢𝑥 |10 − ∫ 𝑣𝑢𝑥 𝑤𝑥 𝑑𝑥 (3)
0 0

As it was mentioned before, the natural boundary conditions set that 𝑣𝑢𝑥 = 0 on both ends. Therefore, this
term will be as:
1 1
∫ 𝑤𝑣𝑢𝑥𝑥 𝑑𝑥 = − ∫ 𝑣𝑢𝑥 𝑤𝑥 𝑑𝑥 (4)
0 0

So, the final form of the weak form will be:

1 1 1 1
∫ 𝑤𝑢𝑡 𝑑𝑥 + ∫ 𝑤𝑎𝑢𝑥 𝑑𝑥 + ∫ 𝑣𝑢𝑥 𝑤𝑥 𝑑𝑥 = ∫ 𝑤 𝑠 𝑑𝑥 (5)
0 0 0 0
 Computational Engineering
2. Spatial and time discretization in matrix form

In this kind of problems (transient convection-diffusion), 𝑢(𝑥, 𝑡) is approximated using separation of

variables through a sum of nodal values 𝑢𝑖 (𝑡), which depend on time, multiplied by shape functions that
only depend on space.
𝑛
ℎ (𝑥,
𝑢(𝑥, 𝑡)~𝑢 𝑡) = ∑ 𝑢𝑖 (𝑡)𝑁𝑖 (𝑥) (6)
𝑖=0

Since equation (6), 𝑢𝑡 and 𝑢𝑥 will be approximated to:

𝑛 𝑛 𝑛

𝑢𝑡 = ∑ 𝑢̇ 𝑗 (𝑡)𝑁𝑗 (𝑥) 𝑢𝑥 = ∑ 𝑢𝑗 (𝑡)𝑁𝑗 ′(𝑥) 𝑢𝑥𝑥 = ∑ 𝑢𝑗 (𝑡)𝑁𝑗 ′′(𝑥)

𝑖=0 𝑖=0 𝑖=0

Regarding the test function, it can be approximated by a test function as:

𝑤 = 𝑁𝑖 (𝑥)

Therefore, it is possible to substitute the approximations presented before into equation (5) such as:

1 𝑛 1 𝑛 1 𝑛 1
∫ 𝑁𝑖 (𝑥) ∑ 𝑢̇ 𝑗 (𝑡)𝑁𝑗 (𝑥) 𝑑𝑥 + ∫ 𝑁𝑖 (𝑥)𝑎 ∑ 𝑢𝑗 (𝑡)𝑁𝑗 ′(𝑥) 𝑑𝑥 + ∫ 𝑣 ∑ 𝑢𝑗 (𝑡)𝑁𝑗 ′(𝑥) 𝑁𝑖 ′(𝑥)𝑑𝑥 = ∫ 𝑁𝑖 (𝑥) 𝑠 𝑑𝑥
0 𝑖=0 0 𝑖=0 0 𝑖=0 0

Rearranging the terms to simplify it:

𝑛 1 1
∑[∫ (𝑁𝑖 (𝑥)𝑢̇ 𝑗 (𝑡)𝑁𝑗 (𝑥) + 𝑁𝑖 (𝑥)𝑎𝑢𝑗 (𝑡)𝑁𝑗′ (𝑥) + 𝑣 𝑢𝑗 (𝑡)𝑁𝑗′ (𝑥)𝑁𝑖′ (𝑥))𝑑𝑥] = ∫ 𝑁𝑖 (𝑥) 𝑠 𝑑𝑥 (7)
𝑖=0 0 0

Performing a space and time discretization will lead to: 𝑴 ∙ 𝒖̇ + 𝑪 ∙ 𝒖 + 𝑲 ∙ 𝒖 = 𝒇. Where each of them can
be found in equation (7):
1 1 1
(𝑒) (𝑒) (𝑒)
𝑀𝑖𝑗 = ∫ 𝑁𝑖 (𝑥)𝑁𝑗 (𝑥) 𝑑𝑥 ; 𝐶𝑖𝑗 = ∫ 𝑁𝑖 (𝑥)𝑎𝑁𝑗′ (𝑥)𝑑𝑥 ; 𝐾𝑖𝑗 = ∫ 𝑁𝑖′ (𝑥)𝑣𝑁𝑗′ (𝑥)𝑑𝑥 ;
0 0 0

1
𝑓𝑖 (𝑒) = ∫ 𝑁𝑖 (𝑥) 𝑠 𝑑𝑥
0

It is important to highlight those matrices 𝑴, 𝑪, 𝑲 and vector 𝒇 are the matrices already assembled (global
notation), containing the contributions of all the nodes of the system. These matrices are already defined
in the notes of the subject as:

ℎ 2 1 𝑣 1 −1
𝑴 = 𝔸𝑴𝒊𝒋 (𝒆) = 𝔸 ( [ ]) ; 𝑲 = 𝔸𝑲𝒊𝒋 (𝒆) = 𝔸 ( [ ])
6 1 2 ℎ −1 1
𝑎 −1 1 (𝒆) 𝑠ℎ 1
𝑪 = 𝔸𝑪𝒊𝒋 (𝒆) = 𝔸 ( [ ]) 𝒇 = 𝔸𝒇𝒊 = 𝔸( [ ])
2 −1 1 2 1

As a method for the time discretization and the time step have to be chosen, it is important to take into
consideration the possible numerical instabilities that can appear even if the problem is convective or
 Computational Engineering
diffusive dominant. For this reason, it is introduced the Péclet number, which is defined in the following
equation:

𝑎ℎ
𝑃𝑒 = (8)
2𝑣

When the Péclet number is greater than 1, instabilities start to appear, meaning that the problem us
convective dominant. There are some methods to deal with instabilities and improve the behavior of the
model. One of the possible solutions can be the reduction of the element size; however, this solution can
drive to high computational cost. To solve this issue, other methos are preferred such as the artificial
diffusion method or the Streamline Upward Petrov-Galerkin (SUPG) method.

To deal with the transient problem of this report, it is considered an explicit method, more specifically the
Explicit Euler Method, which offers a lower computational cost than the implicit method if the stability of
the problem is achieved. As the Explicit method of Euler is unconditionally stable, which means that the
time step condition must be accomplished considering a lumped mass matrix 𝑴𝑳 :

ℎ2
Δ𝑡 ≤ ∆𝑡𝑐𝑟𝑖𝑡𝑖𝑐𝑎𝑙 = (9)
2

Aimed to obtain the element size to model the problem, it is imposed a Péclet number equal to 0,75 (such
that possible instabilities are avoided):

2 ∙ 𝑃𝑒 ∙ 𝑣 2 ∙ 0,75 ∙ 0,01
ℎ= = = 𝟎, 𝟎𝟏𝟓 (10)
𝑎 1

Then, with an element size of 0,015 it is guaranteed that numerical instabilities can appear in the modeling.
On the other hand, to know the time step it should be computed the critical time step of the Explicit Euler
Method.

ℎ2 0,0152
∆𝑡𝑐𝑟𝑖𝑡𝑖𝑐𝑎𝑙 = = = 1,125 ∙ 10−4 − −−→ ∆𝒕 = 𝟎, 𝟓 ∙ 𝟏𝟎−𝟒 (11)
2 2

In order to do not work with critical time step it is prudent to use a lower value. For instance, a value of
𝟎, 𝟓 ∙ 𝟏𝟎−𝟒 is chosen for this parameter.

3. Implementation of MATLAB code

To implement the finite element solution of the unsteady problem with linear finite elements through
MATLAB code. Let’s consider the initial condition is a Gaussian hill centered at x = 0.4 with a maximum of
𝑈𝑚𝑎𝑥 (𝑥, 0) = 1. The governing equation that considered an initial condition Gaussian hill is:

(𝑥 − 𝑏)2
𝑓(𝑥) = 𝑎 ∙ exp (− ) (12)
2𝑐 2

Where the parameter a is the height of the curve peak, so 𝑎 = 𝑈max (𝑥, 0) = 1, 𝑏is the position of the center,
so 𝑏 = 0.4, 𝑎𝑛𝑑 𝑐 is the standard deviation, which controls the width of the bell. In this case, the value has
been chosen arbitrary to 𝑐 = 0,1 aimed to capture properly the evolution of it.
 Computational Engineering

Figura 1. Explicit Euler solution for the transient problem having Gaussian hill initial condition.

From figure 1 it can be seen; the unsteady problem is solved through the Explicit Euler Method. The
advantage of this method, as it was already mentioned, is that provides robust results with no excessive
computational cost if and only if the model is stable. In this problem, the time step and element size are
the ones obtained in the chapter 2 (ℎ = 0,015; ∆𝑡 = 0,5 ∙ 10−4 ).
To conclude, it is possible to say the parabola is tending to zero which, physically, could mean that the
system is tending to the steady state regime. Also, the peak is decreasing as long as the time advances,
the peak is moving to the center of the plot (x=0,5).