Computational Engineering

REPORT 1
-Juan Pablo Wilches Arboleda-
Consider that the flow of water in a porous medium is modelled by means of the partial differential equation:
𝛁 ∙ (𝑘𝛁𝑢) = 0 𝑖𝑛 Ω, (1)

where 𝑢(𝑥, 𝑦) is the piezometric head and 𝑘 is the permeability of the soil. Dirichlet boundary conditions
are assumed at the bottom of the trench (𝑢 = 𝑢𝑖𝑛𝑓 ) and at the terrain level (𝑢 = 𝑢𝑠𝑢𝑝 ). Homogeneous
Neumann boundary conditions are assumed for the rest of the boundary.

𝑢 = 𝑢𝑠𝑢𝑝
(4)
(3)
𝑢 = 𝑢𝑖𝑛𝑓
(5)
(2)
𝜕𝑢 𝜕𝑢
=0 =0 𝑯𝟐
𝜕𝑛 𝑯𝟏 (1)
𝜕𝑛

(6)

𝑩𝟏
𝑩𝟐
Figure 1. problem statement: trench excavated in a soil.

1) Derive the finite element discretization of the problem, starting from the strong form of the
problem to the weak and the discretized form of equation (1). Describe the matrix form of the
final algebraic system of equations to solve and detail the global and elemental forms of the
right-hand-side and the left-hand-side terms.
First, the strong form of the problem and its boundary conditions will have the following form:
𝛁 ∙ (𝑘𝛁𝑢) = 0 𝑖𝑛 Ω
𝜕𝑢
𝑢 = 𝑢𝑖𝑛𝑓 𝑖𝑛 Γ𝐷 ; 𝑢 = 𝑢𝑠𝑢𝑝 𝑖𝑛 Γ𝐷 ; = 0 𝑖𝑛 Γ𝑁
𝜕𝑛

 Weak form:
Aimed to obtain the weak form of the problem, let’s multiply by the test function 𝑤 at both sides of the
equality and integrating over the whole domain:
.
∫ 𝑤∇ ∙ (𝑘∇𝑢) 𝑑Ω = 0
Ω

Then, the integration by parts and the application of the Divergence theorem is applied:
. .
∫Ω ∇w ∙ (𝑘∇𝑢) 𝑑Ω = ∫Γ w ∙ 𝑘∇𝑢 𝑑Γ (2)

Splitting the boundary conditions (either Neuman or Dirichlet) over the domain Γ:
. . .
∫Γ w ∙ 𝑘∇𝑢 𝑑Γ = ∫Γ w ∙ 𝑘∇𝑢𝐷 𝑑Γ𝑁 + ∫Γ w ∙ 𝑘∇𝑢𝑁 𝑑Γ𝑁 (3)
𝐷 𝑁

As the test function is null at the Dirichlet conditions, the first term of the boundary conditions is equal to
0. Therefore, the terms from the equations (2) and (3) can be re-arrange as follows:
0
. . .
∫ w ∙ 𝑘∇𝑢 𝑑Γ = ∫ w ∙ 𝑘∇𝑢𝑁 𝑑Γ𝑁 = ∫ ∇w ∙ (𝑘∇𝑢) 𝑑Ω = 0
Γ Γ𝑁 Ω
 Computational Engineering
So, the final form weak form of the problem reads:
. .
∫Ω ∇w ∙ (𝑘∇𝑢) 𝑑Ω = ∫Γ w ∙ 𝑘∇𝑢𝑁 𝑑Γ𝑁 (4)
𝑁

 Space discretization:
The space can be discretized to approximate the continuous solution of the PDE problem, to perform that
the following approximations are considered for these functions:
𝑢~𝑢ℎ = ∑𝑁
𝐽=1 𝑁𝐽 (𝑥, 𝑦)𝑢𝐽 (𝑥, 𝑦) ; 𝑤~𝑤ℎ = 𝑁𝐼 (𝑥, 𝑦) (5)

Applying the approximation (5) into the weak form (4):

. .
∫Ω ∇𝑁𝐼 ∙ 𝑘∇ (∑𝑁
𝐽=1 𝑁𝐽 𝑢𝐽 ) 𝑑Ω = ∫Γ 𝑁𝐼 ∙ 𝑘∇𝑢𝑁 𝑑Γ𝑁 𝑁
. .
∑𝑁
𝐽=1 ∫Ω ∇𝑁𝐼 ∙ 𝑘∇𝑁𝐽 𝑢𝐽 𝑑Ω = ∫Γ 𝑁𝐼 ∙ 𝑘∇𝑢𝑁 𝑑Γ𝑁 (6)
𝑁

The linear system of equations to solve PDE problems has the following form:

𝐾 ∙ 𝑢̅ = 𝑓 ̿
Where the terms for the global matrix and the force vector in global form can found from equation (6):
. .
𝐾𝐼𝐽 = ∫Ω ∇𝑁𝐼 ∙ 𝑘∇𝑁𝐽 𝑑Ω ; 𝑓𝐼 = ∫Γ 𝑁𝐼 ∙ 𝑘∇𝑢𝑁 𝑑Γ𝑁 (7)
𝑁

As the final form of the system of equation is known for the problem presented, the isoparametric
transformation can be implemented according to the following equations:
𝑛 𝑛

𝑥 = ∑ 𝑥𝑖 𝑁𝑖 (𝜉, 𝜂) ; 𝑦 = ∑ 𝑦𝑖 𝑁𝑖 (𝜉, 𝜂)
𝑖=1 𝑖=1

1 1
𝑁1 (𝜉, 𝜂) = (𝜉 − 1)(𝜂 − 1) ; 𝑁2 (𝜉, 𝜂) = − (𝜉 + 1)(𝜂 − 1)
4 4
1 1
𝑁3 (𝜉, 𝜂) = (𝜉 + 1)(𝜂 + 1) ; 𝑁4 (𝜉, 𝜂) = − (𝜉 − 1)(𝜂 + 1) (8)
4 4

The Jacobian is described as:

𝑛 𝑛
⎡∑ 𝑥 𝜕𝑁𝑖
𝜕𝑥 𝜕𝑦 ∑ 𝑦𝑖
𝜕𝑁𝑖 ⎤
⎡ ⎤ ⎢ 𝑖
𝜕𝜉 𝜕𝜉 ⎥
⎢ 𝜕𝜉 𝜕𝜉 ⎥ 𝑖=1 𝑖=1
𝑱 (𝜉, 𝜂) = ⎢
= 𝑛 𝑛
⎥
⎢𝜕𝑥 𝜕𝑦⎥ ⎢ 𝜕𝑁𝑖 𝜕𝑁𝑖 ⎥
⎣𝜕𝜂 𝜕𝜂⎦ ⎢∑ 𝑥𝑖 ∑ 𝑥𝑖 ⎥
⎣ 𝑖=1 𝜕𝜂 𝑖=1
𝜕𝜂 ⎦

∇𝑥𝑦 = 𝑱−1 ∇𝜂𝜉 ; 𝑑𝑥𝑑𝑦 = |𝑱|𝑑𝜉𝑑𝜂

Once the isoparametric is introduced, it is possible to find the elemental form of matrix K:
1
𝐾𝐼𝐽 = ∫ 𝑘 ∇𝑁𝐼 (𝑥, 𝑦) ∇𝑁𝐽 (𝑥, 𝑦) 𝑑𝑥𝑑𝑦
0
1
𝐾𝑖𝑗 = ∫ 𝑘 𝑱−𝟏 ∇𝑁𝑎 (𝜉, 𝜂) 𝑱−𝟏 ∇𝑁𝑏 (𝜉, 𝜂) |𝑱| 𝑑𝜉𝑑𝜂
−1

= ∑41 𝑘 𝑱−𝟏 (𝑧𝑔 ) ∇𝑁𝑎 (𝑧𝑔 ) 𝑱−𝟏 (𝑧𝑔 ) ∇𝑁𝑏 (𝑧𝑔 ) |𝑱| 𝑤𝑔 (9)

It is important to highlight that 𝑧𝑔 are the Gauss points and 𝑤𝑔 are the weights defined as:
1 1 1 1
𝑧𝑔 = [− − ] ; 𝑤𝑔 = [1 1 1 1]
√3 √3 √3 √3
 Computational Engineering
2) Complete the questions asked in the computer assignment “Equilibrium problems”. When
identifying specific terms, make references to the equations derived in the section (a) above.
a) Discuss the physical meaning of the boundary condition for each of the six segments in the
boundary.
According to figure 1, the boundary conditions are numbered and described as follows:
 Condition 1: this boundary is part of the Neuman conditions as it belongs to the continuous medium of
the saturated soil. It is exactly the middle of the excavated trench.
 Condition 2: It is the bottom of the trench, where is applied one of the Dirichlet conditions, which means
that the piezometric head is known: 𝑢𝑖𝑛𝑓 .
 Condition 3: Represents the lateral wall of the trench where the is not water flow (Neuman condition).
 Condition 4: It is the top of the trench, where is applied one of the Dirichlet conditions, which means
that the piezometric head is known: 𝑢𝑠𝑢𝑝 .
 Condition 5: This part of the soil is assumed very far from the area of influence, then it is considered
𝜕𝑢
that any water flux appears ( = 0).
𝜕𝑛
 Condition 6: The bed of the saturated soil is considered as impervious, where there isn’t any water flux.
b) Identify in the file main.m the code lines where the following operations are done:
 Definition of parameter values of the problem: In order to define the main parameters of the problem,
the script “main.m” establish the domain dimensions, permeability, reference element and tolerance
applied in finding the mesh boundaries in lines 9, 10, 13, 16 and 19, respectively.
 Computation of shape functions in the reference space: The shape functions in the reference space are
computed in the line 24, where the function “ReferenceElement2D” is called to perform calculation. In
the first chapter of the report these functions can be found in equation (8).
 Mesh creation: The mesh is created in line 27 and 28. In line 27, the code asks for a value of h, which
is the desired element size of the mesh. On the other hand, line 18 calls the function “createMeshtrench”
that creates the mesh according to the initial parameters and the size of the element h.
 Identification of nodes belonging to the trench boundary: In line 38 the nodes belonging to the trench
boundary are found through MATLAB function “find”.
 Discretization, integration and assembling of equation (1): The process of discretization, integration and
assembling of equation (1) take place in line 51 and 52, where the function “compute_system” is called.
Within this function is coded another function named “calcul_elemental” which complements the task
of discretization and integration. In chapter 1 this process can be found in equation (6) and (7).
 Application of Dirichlet boundary conditions depicted in figure 1: The Dirichlet boundary condition are
applied in lines from 55 to 57 to compute the force vector f. In chapter 1 this process can be found in
equation (3).
 Application of Neumann boundary conditions depicted in figure 1: No in the code. However, these
conditions are applied in chapter 1 in equation (3).
 Solution of the linear system of equations: The solution of the linear system of equations can be found
in line 60, with the reduced matrix K and the force vector f.
c) More specifically, regarding the file compute_system.m, identify the code lines and structures
(e.g. loops) that you consider essential to compute the following elements:
𝜕𝑁𝑖 𝜕𝑁𝑖
 Gradient of shape function Ni in the reference coordinates ∇𝑁𝑖 (𝜉. 𝜂) = ( , ): Lines 33 and 34. To
𝜕𝜉 𝜕𝜂
perform the gradients is used a for loop structure.
𝜕𝑁𝑖 𝜕𝑁𝑖
 Gradient of shape function Ni in the physical coordinates ∇𝑁𝑖 (𝑥. 𝑦) = ( , ): Line 40. To perform the
𝜕𝑥 𝜕𝑥
gradients is used a for loop structure.
 Computational Engineering
.
 Component (i, j) of the matrix: 𝑘𝑖𝑗 = ∫Ω ∇𝑁𝑖 ∇𝑁𝑗 𝑑Ω for the element Ωe: Line 43. For loop
𝑒
 Global matrix K: Line 14. It is needed to use a for loop to assess the matrix along each element node.
d) Solve the problem for 𝑯𝟏 = 𝒖𝒊𝒏𝒇 = 𝟓 𝒎, 𝑯𝟐 = 𝒖𝒔𝒖𝒑 = 𝟏𝟎 𝒎, 𝑩𝟏 = 𝟓 𝒎, 𝑩𝟐 = 𝟐𝟎 𝒎 and 𝒌 = 𝟏𝟎−𝟒 𝒎/𝒔.
Compute the flow in the bottom of the trench Q as accurately as possible. What is the achieved
accuracy?
Aimed to perform the flow in the bottom of the trench Q as accurately as possible, the element mesh size
(h) is decreased until the error converges to an acceptable value, under 1% is desirable. In the table 1 is
shown the different values of h that were used, and the results obtained for the water flow Q.
h (m) Q (m3/s) Error 0,14
h vs error 0,12
5 0,000537
0,1
2,5 0,000474 0,117
0,08
1 0,00043 0,093 0,06
0,5 0,000411 0,042 0,04

0,25 0,0004 0,027 0,02

0
0,125 0,000393 0,016 6 5 4 3 2 1 0
0,1 0,000392 0,004
Table 1. Results with different values of h Figure 2. Values of the element mesh size vs the relative error.

From the figure 2, it can be seen how the relative error of the flow calculation decreases if the value of the
element mesh size decreases too, it seems to be a quadratic convergence after the value of 2,5 m for the
h. An error lower than 1% (0,4%) is reached with an element mesh size of 0,1 m, which is quite reliable to
perform the computation.
e) What changes would you do to compute Q more accurately than in d) with the same computing
resources?
Considering an equilibrium between increasing the accuracy in the computation of Q and the
computational cost, it is proposed to perform a discretization of the mesh by parts, which consists of create
a more discretized mesh, lower element mesh size, close to the corner where the vertical wall of trench
and the soil are in contact with each other and, on the other hand, a higher element mesh size at the rest
of the soil were the gradient of water flow is not as high as in the previous area. This solution will probably
increase the accuracy in the computation of Q keeping the same computing resources of d statement.

Figure 3. Piezometric head in the soil