Aula Teórica 9&10
Equação de Evolução. Exemplos. Caso
1D
Evolution Equation


d
dV 
dV   v .n dA


t vc
dt sistema
surface
The rate of change inside the system is the
(Production-Destruction) + (diffusion exchange).
Designating Production – Destruction by (Sources –
Sinks) and knowing that:


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




 
  S o  S i 



x
j 

• 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 

x2
• For the control volume:
c
• What is the sign of u1 ?
x1
• 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)?
x1
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 concentration 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?
e) If it was a concentration without adsortion, what would be the final profile?
y
c

ux

x x x
a)
b)
c)
d)
e)
f)
 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. The total
derivative is negative (and so is the advective derivative)
yes the system could be stationary, but the fluid would loose property along the flow.
If it is stationary and conservative what would be the sign of the longitudinal gradient component?
If it was a concentration this profile could not maintain. The gradient should become null at the solid boundary (no
diffusive flux). The final profile would be uniform. Transversal diffusion would homogenize it.
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?
• The Reynolds Theorem:


d
dV 
dV   v .n dA


t vc
dt sistema
surface
• The lagrangian formulation:


 
d
dV  (So  Si)    c.n d A

dt sistema
surface
The Eulerian formulation


 


dV    u .n dA    c.n dA  ( So  Si )

t vc
surface
surface



 


dV    u .n    c.n dA  ( So  Si )

t vc
surface
• This evolution equation states that “the rate
of change inside a control volume balances
the fluxes, plus sources minus sinks”.
• In differential form the equation becomes:
c
c

 u j

t
x j x j
c


t
x j
 c

 x
j


 u j c   c

x j


  S o  Si 



  S o  Si 


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?
Flux calculation
• 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”.
Definição de velocidade
Cx
Cx+∆x
• A figura representa moléculas de dois fluidos em
repouso. A velocidade mede o volume de
moléculas que passa por unidade de área.
• Se a velocidade for nula, o volume que passa num
sentido é igual ao que passa no sentido contrário.
Difusão
• Mas as moléculas têm movimento browniano
e por isso - num fluido - estão sempre a mudar
de posição relativa.
• Se as moléculas que estão de um lado da
superfície forem iguais às que estão do outro
lado, o saldo é estatisticamente nulo.
• Se a concentração for diferente, então existirá
um saldo com um fluxo resultante orientado
da concentração maior para a menor.
Fluxo difusivo por unidade de área
 d  cl  cl l ub
c
cl  cl l   l
l
c
 d  l.ub
l
Na direcção “x”:
 dx
c
 
x
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.