Aula Teórica 13 & 14
Evolution Equation. Examples and 1D
case
Reynolds Theorem
States that:


d
dV 
dV   v .n dA


t vc
dt sistema
surface
Based on this theorem we have deduced that the
difference between the total (or material) derivative and
the local time derivative is the convective derivative:

d


 uj
t
dt
x j
d



 uj
dt
t
x j
Why does properties change inside
a Material System?
• Because of Source and Sink terms (e.g.
chemical reactions, biological activity)
• Because of diffusion in presence of
concentration gradients.
Diffusion
Figures below represent 2 material systems, one fully white and the other
fully Black separated by a diaphragm. The top figures represent the
molecules (microscopic view) and the figures below the macroscopic view.
When the diaphragm is removed the molecules from both systems start to
mix and we start to see a grey zone between the two systems (b) at the end
everything will be grey (c). During situation (b) we there is a diffusive flux of
black molecules crossing the diaphragm section. This flux cannot be
advective because velocity is null.
(a)
(a)
(b)
(b)
(c)
(c)
Diffusivity
When the diaphragm is removed
molecules move randomly. The net flux
is the diffusive flux.
Cx
Cx+∆x
 d  cl  cl l ub
c
 d  l.ub
l
The flux of molecules in each sense is
proportional to the concentration and to
the individual random velocity:
But,
c
cl  cl l   l
l
Diffusivity is the product of the displacement
length and the molecule velocity. This velocityis in
fact the difference between the molecule velocity
and the average velocity of the molecules
accounted for in the advective term.
Ver texto sobre propriedades dos fluidos e do campo de velocidades
Diffusivity
• Diffusivity is definide as:
  l.ub
• Where ub is the molecule velocity part not resolved
(or included) in our velocity definition. In a laminar
flow is the brownian velocity while in a turbulent flow
is the turbulent velocity, a macroscopic velocity that
we can see in the tubulent eddies.
• l is the lenght of the displacement of a molecule
before being disturbed by another molecule (or of a
portion of fluid in a turbulent flow). When the
molecule hits another molecule it gets a new velocity.
• Diffusivity dimensions are: L2T 1
Diffusive Flux
• Is the flux due to difusivity and property
gradient:
 Dif 
 



c
n j dA
 c .n dA     

x j 
A

A

• The sense of the diffusive flux is opposit to the
sense of the gradient.
• Diffusive flux is nul if there is no gradient.
Evolution Equation
The Reynolds theorem states that:


d
dV 
dV   v .n dA


t vc
dt sistema
surface
Based on this theorem we have deduced that the
difference between the total (or material) derivative and
the local time derivative is the convective derivative:


Diffusion      .ndA
 










dV

S

S


v
.
n




.ndA
0
i


t vc
surface
Equação de evolução



t
x j



 u j  

x j


  So  Si 


Ou:
d




 uj
 
dt
t
x j
x j



 

x j


  So  S


• Diffusion plus advection
• Where is diffusive flux maximum?
Analysis of the evolution equation
dc c
c

 uj

dt t
x j x j
 c 

  S o  Si 
 x 
j 

• For the control volume:
c
u
x2
• What is the sign of 1 x ?
1
• What does that mean physically (how does the
advective flux vary with x1)?
• What is the relative value of diffusive flux in the lower
and upper faces of the control volume?
• How does diffusive flux along x2 contribute to the
concentration inside the control volume?
• If the flow is stationary and the material is
conservative what is the relation between advection
and diffusion? what is the divergence of the total flux
(advection + diffusion)?
c

 
t
x j
x1



c
u jc  
  So  Si 



x
j


Answers
•
•
•
•
•
•
•
c
0
The advective term is negative. The velocity is positive and
x1
In other words, the advective flux u1c decreases with x1, i.e., the quantity entering is bigger
than the quantity leaving (per unit of area).
The diffusive flux in the lower face is about zero. If located exactly over the symmetry line,
the flux will be exactly zero, because the gradient is null  dif   c
x2
On the upper face the diffusive flux is positive, i.e., the material is transported along the axis
x2 because the concentration is higher inside the control volume than above the control
volume.
The divergence of the vertical diffusive flux contributes to decrease the concentration inside
the volume.
The divergence of the horizontal diffusive flux is negative because the gradient on the left
side of the volume is higher than on the right side. Horizontal advection is however the main
mechanism to increase the concentration inside the control volume.
If the flow is stationary the horizontal advection plus the horizontal diffusion
balance
the
 
c 


vertical diffusion. The divergence of both fluxes would be zero x  u c  x   0


If the velocity was increased, the concentration inside the control volume would tend to
increase because the quantity entering would increase and thus the quantity leaving would
have to increase too.
j
•
j
j
Case of concentration
Consider two parallel plates and a property with a parabolic type distribution (blue)
with maximum value at the center. Assume a stationary flow and a conservative
property.
a) Draw a control volume and indicate the fluxes (sense and relative magnitude).
b) Where would the diffusive flux be maximum? Could this property be a
concentration? What kind of property could it be?
c) What is the sign of the material time derivative?
d) If the material is conservative (no sink or source), can the profile be completely
developed?
y
c

ux

x x x
a)
b)
c)
d)
 c  

 
 x x  x y
 c 


 x 
y 

Diffusive flux increase with “r”. Advective flux depends on the longitudinal gradient, not yet known.
Property gradient is maximum at the boundary. This means that there the diffusive flux is maximum and
consequently the property can pass through the boundary. The property can be a concentration only if there is
adsortion at the boundary. It could be a temperature or momentum as well.
Property is being lost across the boundary and consequently there is a longitudinal negative gradient , unless if
there is a source. The total derivative is negative (and so is the advective derivative)
No. Without a source there would be a negative longitudinal gradient. Otherwise it would depend on the source.
Nitrogen Transport in the Baltic Sea
• The figure shows the
distribution of nitrate
discharged by the river
Oder after 4, 8, 12 and
16 years of emission.
• In the region next to
the river mouth the
concentration becomes
constant after 8 years.
What does that mean?
• Is advection important
in this system?
1D case
• Vamos supor o caso de um reservatório com a forma de um
canal rectangular de 10 m de largura e 200 de comprimento,
com velocidade nula. A concentração é elevada na zona
central e nula na generalidade do canal. O material é
conservativo.
• Escreva a equação que rege a evolução da concentração.
dc c
c

 uj

dt t
x j x j
• E se existisse decaimento?.
 c 

  S o  Si 
 x 
j 

With first order decay
dc c   c 

    kc
dt t x  x 
• How will the concentration evolve in each case?
• In case of decay, concentration
are lower and tend to zero.
t0
t1
t2
t∞
If there is a lateral discharge?
• Concentration will grow because of the
discharge, will get homogenized because of
diffusion and will decay because of decay.
• When discharge balances total decay, the
system will reach equilibrium.
Solution of the problem
• The first case (pure diffusion) has analytical
solution, but not the others.
• How to solve the problem numerically?
Let’s split the channel into
elementary volumes
Cell i-1
Cell i
Cell i+1
And apply the equation to each of them:


 

cdV     c.n dA   kc dVol

t vc
If the property can be considered uniform inside each cell and along each surface:
 cit  t  cit
Vol 
 t
Vol  Ax

 ci  ci 1 
 ci 1  ci
   A 
  A 
 x 
 x


  Vol kc 

Where A is the area of the cross section between the elementary volumes,
assumed constant in the academic example.
Numerical solution
 cit  t  cit
Vol 
 t
Vol  Ax

 ci  ci 1 
 ci 1  ci 
   A 
  A 
  Vol kc 
 x 
 x 

• Dividing all the equation by the volume:
 cit  t  cit

 t

 ci  ci 1 
 ci 1  ci
   
  
2
2

x

x





  kc

• In this equation we have two variables, cit  t , cit and an extra
variable, c, which time we have not defined. cit Are the actual
t  t
c
concentrations and i are the concentrations that we want
to calculate. The concentration c is the concentration used to
calculate the fluxes and decay. What is it?
Temporal descretisation
• The flux equation gives quantity per unit of
time. When we equation as:
 cit  t  cit

 t

c c 
c c
    i 2i 1     i 1 2 i
 x 
 x


  kc

• We are computing the amount that crossed
the surface during a time period (t) and
dividing it by the length of the time. The most
convenient time to allocate to c is the middle
of the time interval:
 cit  cit  t
c  
2




The equation becomes
 cit  t  cit

 t

 k t  t
  ci  cit 
 2


  cit  t  cit1t    cit1t  cit  t 

2
 cit  t  cit

 t


  cit  cit1    cit1  cit
   
  
2
2
2

x
2

x




x
2
  
 2
   cit  t  cit1t
  
x 2
 2
x
2


   cit1t  cit  t
  
x 2
 2
 kcit  t
  cit  cit1    cit1  cit
 
  
  
2
2
2  x  2  x 2

t
 kcit
 
 2
 2t k  t  t t t  t t t
 2t k  t t t
t  t
c

1


c

c

c

 ci 
c


1 
i 1
i
2
2 i 1
2 i 1
2
2 i 1
2x 2
2

x
2
2

x
2

x
2

x
2
2

x





D t  t 
k
D
D
k
D

ci 1  1  D  cit  t  cit1t  cit1  1  D  cit  cit1
2
2
2
2
2
2


This equation has 3 unknowns. The all channel requires the
resolution of a system of equations in each time step.
Explicit calculation
• If we had assumed:
c  cit
• We would have got:
cit  t  Dcit1  1  2 D  k cit  Dcit1
• That can be solved explicitly, but has stability
problems.
Implicit calculation
• If we had assumed:
c  cit  t
• We would have got:
 Dcit1t  1  2 D  k cit  t  Dci1t  cit
• That has less computation than the “average”
approach, but still requires the resolution of a
system of equations.
Results (D<0.5)
D=
DT=
20
difusividate
1
k=
0
i-3
0
0
100
dx
Time
0.25
i-2
i-1
i
i+1
i+2
i+3
0
0
0
0
1
0
0
0
0
100
0
0
0
0.25
0.5
0.25
0
0
0
200
0
0
0.06
0.25
0.38
0.25
0.06
0
0
300
0
0.02
0.09
0.23
0.31
0.23
0.09
0.02
0
400
0
0.03
0.11
0.22
0.27
0.22
0.11
0.03
0
500
0.01
0.03
0.12
0.21
0.25
0.21
0.12
0.04
0.01
The results evolve as expected.
0
Results (D=0.75)
D=
DT=
300
dx
20
difusividate
Time
1
k=
i-3
0
0
0.75
0
i-2
i-1
i
i+1
i+2
i+3
0
0
0
0
1
0
0
0
0
300
0
0
0
0.75
-0.5
0.75
0
0
0
600
0
0
0.56
-0.8
1.38
-0.8
0.56
0
0
900
0
0.42
-0.8
1.83
-1.8
1.83
-0.8
0.42
0
1200
0.32
-0.8
2.11
-2.9
3.65
-2.9
2.11
-0.8
0.32
1500
-0.79
0.71
-3.9
5.77
-6.2
5.77
-3.9
2.24
-0.8
0
The results are strange. Concentration bigger than initial are
obtained. Negative concentrations are also obtained….
The method is unstable
Stability
• The method became unstable because the
parenthesis becomes negative when D>0.5;
cit t  Dcit1  1  2D  k cit  Dcit1
• When the parenthesis changes sign, the effect of ci at time t
changes. If it is positive the larger is the initial concentration
the larger is the subsequent concentration. If it is negative the
larger is the initial concentration the smaller is the
subsequent, which is physically impossible.
• This is the limitation of the explicit methods. They are
conditionally stable. In pure diffusion problems the D (the
diffusion number) must be smaller than 0.5. If other
processes exist (e.g. k≠0) D must be smaller.
Implicit methods stability
• Implicit (as well as semi-implicit) methods are
more difficult to program, but do not have
stability limitations allowing larger values for
D, meaning that we can combine small grid
sizes with large time steps, while in explicit
methods the time step is associated to the
square of the spatial step.
How large should diffusivity be?
• We have to reassess the concept of diffusivity:
“Diffusivity is the product of a nonresolved velocity part by the length of the
displacement of the fluid due to that
velocity”.
Value of diffusivity in our problem
• When we assume the velocity to be uniform in
the cross section area, we are neglecting
spatial variability. If the profile was parabolic
the average velocity would be half of the
maximum velocity. We would be neglecting a
eddies with velocity similar to the average
velocity and radius of half channel width.
• Diffusivity should be of the order of: UL/2.
Stability
• The system gets unstable if D>0.5
D
t
x 2
x
Ut 
2

Uxt
 0 .5
2
x
• This equation shows that the length of the
displacement due to diffusion must be smaller
than half of the cell length.
Initial Conditions
• Initial values to be provided in each cell
Boundary conditions
• Boundary conditions specify how the fluid
interacts with the surrounding environment.
– Solid boundaries
– Open boundaries,
• In geophysics boundaries can be vertical of
horizontal (Top and bottom). Gravity make
them important.
• One can impose property’s values or fluxes
(advection and/or diffusion).
Input data
•
•
•
•
Geometry,
Flow properties,
Process parameters,
Execution parameters.