Department of Physics Umeå University
Swedish Defence Research Agency (FOI) 26 March 2015
Tobias Brännvall Spring 2015
Master’s thesis, 30 hp Physical Engineering, 300 hp
Source Term Estimation in the
Atmospheric Boundary Layer
Using the adjoint of the Reynolds Averaged Scalar
Transport equation
Tobias Brännvall
Abstract
This work evaluates whether the branch of Reynolds Averaging in Computational Fluid Dynamics can be used to, based on real field measurements, find the source of the measured gas in question. The method to do this is via the adjoint to the Reynolds Averaged Scalar Transport equation, explained and derived herein. Since the Inverse is only as good as the main equation, forward runs are made to evaluate the turbulence model.
Reynolds Averaged Navier Stokes is solved in a domain containing 4 cubes in a 2x2 grid, generating a velocity field for said domain. The turbulence model in question is a union of two modifications to the standard two equation k-ε model in order to capture blunt body turbulence but also to model the atmospheric boundary layer. This field is then inserted into the Reynolds Averaged Scalar Transport equation and the simulation is compared to data from the Environmental Flow wind tunnel in Surrey. Finally the adjoint scalar transport is solved, both for synthetic data that was generated in the forward run, but also for the data from EnFlo.
It was discovered that the turbulent Schmidt number plays a major role in capturing the dispersed gas, three different Schmidt numbers were tested, the standard 0.7, the unconventional 0.3 and a height dependent Schmidt number. The widely accepted value of 0.7 did not capture the dispersion at all and gave a huge model error. As such the adjoint scalar transport was solved for 0.3 and a height dependent Schmidt number.
Institutionen för fysik Umeå Universitet
Totalförsvarets Forskningsinstitut (FOI) 26 Mars 2015
Tobias Brännvall Vt 2015 Examensarbete, 30 hp Civilingenjör-Teknisk Fysik, 300 hp
Källtermsuppskattning i det
atmosfäriska gränsskiktet
Med hjälp av den adjungerade Reynolds tidsmedlade
Skalära Transportekvationen
Tobias Brännvall
Sammanfattning
Detta arbete utvärderar hurvida Reynolds medelvärdesmodellering inom flödessimuleringar kan användas till att finna källan till en viss gas baserat på verkliga mätningar ute i fält. Metoden går ut på att använda den adjungerade ekvationen till Reynolds tidsmedlade skalära transportekvationen, beskriven och härledd häri. Då bakåtmodellen bygger på framåtmodellen, måste såleds framåtmodellen utvärderas först.
Navier-Stokes ekvationer med en turbulensmodell löses i en domän, innehållandes 4 kuber i en 2x2 orientering, för vilken en hastighetsprofil erhålles. Turbulensmodellen som användes är en union av två olika k-ε modeller, där den ena fångar turbulens runt tröga objekt och den andra som modellerar atmosfäriska gränsskiktet. Detta fält används sedan i framåtmodellen av skalära transportekvationen, som sedan jämförs med körningar från EnFlo windtunneln i Surrey. Slutligen testkörs även den adjungerade ekvationen, både för syntetiskt data genererat i framåtkörningen men även för data från EnFlo tunneln.
Då det visade sig att det turbulenta Schmidttalet spelar stor roll inom spridning i det atmosfäriska gränsskiktet, gjordes testkörningar med tre olika Schmidttal, det normala 0.7, det väldigt låga talet 0.3 samt ett höjdberoende Schmidttal. Det visade sig att det vanligtvis använda talet 0.7 inte alls lyckas fånga spridningen tillfredställande och gav ett stort modellfel. Därför löstes den adjungerade ekvationen för 0.3 samt för ett höjdberoende Schmidttal.
Acknowledgements
Table of content
Nomenclature ... 1 Introduction ... 1 1.1 Background ... 1 1.2 Problem Description ... 1 1.3 Purpose ... 1 2 Theory ... 2 2.1 Reynolds averaging ... 22.2 The incompressible Reynolds Averaged Navier Stokes (RANS) ... 3
2.3 Reynolds averaged scalar transport equation ... 4
2.4 The Adjoint ... 5
2.5 Adjoint of Reynolds averaged transport equation ... 6
2.6 A Measurement and the adjoint duality ... 8
2.7 Source strength ... 9
2.8 The scalar product in a computational grid ... 10
2.9 Source term estimator ... 11
2.10 The Atmospheric Boundary Layer (ABL) ... 12
2.11 K-ε Model for the atmospheric boundary layer ... 12
2.12 Boundary conditions ... 14
2.13 The Turbulent Schmidt number ... 15
2.14 Formal definition of the Adjoint Equation ... 16
3 Outline of numerical solution in OpenFOAM ... 17
3.1 Implementation of the Adjoint equation ... 17
3.2 Numerical Schemes and Algorithms ... 18
3.3 Pre-processing - Parallel computing ... 19
3.4 Post processing ... 19
4 The domain ... 20
4.1 The simple array ... 20
4.2 Mesh ... 20
5 Results ... 22
5.1 Wind profile and turbulence – mimicking EnFlo ... 23
5.2 Flow over the simple array ... 25
5.4 Synthetic data test ... 31
5.5 ENFLO simple array data ... 33
5.5.1 Five measurements ... 33 5.5.2 Three measurements ... 34 6 Discussion ... 38 7 Concluding words ... 41 8 Bibliography ... 42 A. Appendix ... i
A.1 Dispersion comparisons of different Schmidt numbers ... i
Nomenclature
𝐴 Differential operator
𝑐 Concentration of neutral gas 𝐶 Mean Concentration of neutral gas 𝑐′ Fluctuating Concentration of neutral gas 〈𝑐𝑗〉 Measured concentration
𝐶∗ Adjoint concentration
𝐶1 Turbulence model constant
𝐶2 Turbulence model constant 𝐶𝜇 Turbulence model constant
𝐷 Diffusivity of neutral gas
𝐷𝑡 Turbulent diffusivity of Reynolds averaged neutral gas 𝛿 Dirac delta function
𝜖 Turbulent energy dissipation 𝜙, 𝜓 Variable
𝜙′, 𝜓′ Fluctuating part of variable
𝛷, 𝛹 Mean quantity of variable ℎ𝑟𝑒𝑓 Reference height
𝜅 Von Karman constant 𝑘 Turbulent kinetic energy
𝑘0 Constant describing standard deviation
𝐾 Sum of diffusivity and turbulent diffusivity of Reynolds averaged neutral gas 𝐿 Transport differential operator
𝐿∗ Adjoint Transport differential operator
𝜈 Fluid viscosity
𝜈𝑡 Turbulent viscosity of fluid
𝛺 Spatial Domain, Ω = (𝑥, 𝑦, 𝑧) = (𝑥1, 𝑥2, 𝑥3)
𝛺𝑘 Node k in spatial domain 𝛺𝑗 Measurement position
𝛺̃ Vorticity of fluid
𝑃(𝑆̂) Probability distribution of source 𝑝 Pressure
𝑃 Mean pressure 𝑝′ Fluctuating pressure
𝑃𝑘 Turbulent production rate
𝑄 Source strength
𝑄̂ Least square estimation of source strength 𝜌 Fluid density
𝑆̂ Source position 𝑠 Source of neutral gas
𝑆 Reynolds averaged source of neutral gas 𝑆∗ Adjoint source
𝑆̃ Strain rate of fluid
𝑆𝑐 Schmidt number of neutral gas
𝑆𝑐𝑡 Turbulent Schmidt number of Reynolds averaged neutral gas 𝜎𝑘 Prandtl number for turbulent kinetic energy
𝜎𝜖 Prandtl number for turbulent dissipation 𝑡 Forward time variable
𝑢𝑖 Velocity in i = (1,2,3) = (x,y,z) direction
𝑈𝑖 Reynolds Averaged Velocity in i = (1,2,3) = (x,y,z) direction 𝑢𝑖′ Fluctuating velocity
𝑢∗ Friction velocity
𝑢𝑟𝑒𝑓 Reference velocity
𝑣 , 𝑤 Real function
𝑥𝑖 Spatial coordinate (1,2,3) = (x,y,z) 𝜉 Normalization factor
𝑦+ Dimensionless distance to wall
1
1 Introduction
This master thesis is a computational fluid dynamics (CFD) focused work with regards to inverse modelling. The main focus is regarding the inverse to the Reynolds averaged scalar transport, but with some quantitative analysis regarding the forward model and how it behaves compared to wind tunnel measurements. This work was done under the Swedish Defence Research Agency (FOI) in a cooperative research project regarding dispersion of gases in the atmosphere. Due to the source term having both strength and a position, along with few measurements in the real world, makes the problem of Source Term estimation (STE) an ill-posed problem.
1.1 Background
Normally when performing dispersion modelling, the source of the pollutant is known. However if it is not known it needs to be estimated based from field measurements. One real life example of this occurring is the Chernobyl accident, which happened the 26th of April in 1986 during which radioactive observations were performed all over Europe. The nuclear power plant Forsmark in Sweden, had their alarms set off from too high radiation levels on the 28th and started an evacuation as they thought it was their power plant that had an accident as the soviet did not make the accident public prior to this. Had source term estimation been an active subject at the time, it would have been possible to perform back-trajectory simulations as demonstrated by Robertson. [1]
Currently there is much research being performed to create well performing estimator models which can narrow down the position (both in space and in time) of the source. Dispersion in urban environments requires advanced CFD model, which in turn leads to more complex estimation models. One project in where these models can be expanded upon Is in the MODITIC project, an European Defence Agency (EDA) project with collaborators from Sweden, Norway, England and France. This project exists in order to further enhance the knowledge of dispersion in urban environments, both forward modelling and backtracking.
1.2 Problem Description
The task is to test if Reynolds averaged turbulence models could be used in source term estimation, i.e. finding the source of a gas plume that has been detected by field measurements. This can be done using the adjoint of the scalar transport, then testing both for synthetic data and also for data supplied from the Environmental Flow (EnFlo) wind tunnel at the University of Surrey, in Gilford, England. The work is to be implemented in OpenFOAM 2.3.0, and the gas in question is a neutral gas, meaning it has the same density as the carrier gas (in this case air) and as such gravity does not affect it.
1.3 Purpose
2
2 Theory
2.1 Reynolds averaging
The main idea behind Reynolds time-averaging is to decompose any variable 𝜙(𝑥, 𝑡) into a sum of a mean component 𝛷(𝑥, 𝑡), and a fluctuating component 𝜙′(𝑥, 𝑡) [2] [3]
𝜙(𝑥, 𝑡) = 𝛷(𝑥, 𝑡) + 𝜙′(𝑥, 𝑡). (2.1.1)
The Reynolds time average over the time τ is defined as
𝛷(𝑥, 𝑡) = 〈𝜙(𝑥, 𝑡)〉 = 1
2𝜏∫ 𝜙(𝑥, 𝑡)𝑑𝑡
𝑡+𝜏
𝑡−𝜏 ,
(2.1.2)
Under the assumption that the time scale of the turbulent fluctuations is much less than τ, and that 𝜏 is much less than the time scale relative to the mean flow. Within the topic of averaging, the ergodic hypothesis is used, which states that the ensemble average, time average and volume average are all equal. As such the average defined in (2.1.2) is still considered time dependent, see [3] [4] [5] for more definitions and properties of the average and ergodicity. Under the ergodic assumption for the average, the following rules apply to the Reynolds time average.
1. The time average of a constant k is equal to the constant
𝑘 = 〈𝑘〉 = 𝑘. (2.1.3)
2. The time average of a time average quantity is the average itself
〈𝛷(𝑥, 𝑡)〉 = 〈〈𝛷(𝑥, 𝑡)〉〉 = 𝛷(𝑥, 𝑡). (2.1.4) 3. The time average is distributive as it is a definite integral, as such
〈𝜙(𝑥, 𝑡) + 𝜓(𝑥, 𝑡)〉 = 〈𝜙(𝑥, 𝑡)〉 + 〈𝜓(𝑥, 𝑡)〉 = 𝛷(𝑥) + 𝛹(𝑥). (2.1.5) 4. The time average of a mean quantity times a fluctuating quantity is zero as it acts as a
constant times the timed average of a fluctuating quantity
〈𝜙′(𝑥, 𝑡)𝛹(𝑥, 𝑡)〉 = 〈𝜙′(𝑥, 𝑡)〉𝛹(𝑥, 𝑡) = 0. (2.1.6) 5. The time average of the product of two variable quantities 𝜙(𝑥, 𝑡) and 𝜓(𝑥, 𝑡) is
〈𝜙𝜓〉 = 〈(𝛷 + 𝜙′)(𝛹 + 𝜓′)〉 = 〈𝛷𝛹〉 + 〈𝛷𝜓′〉 + 〈𝜙′𝛹〉 + 〈𝜙′𝜓′〉 = 𝛷𝛹 + 〈𝜙′𝜓′〉. (2.1.7)
Where in the last equality (2.1.6) was used. 6. For a derivative, the average is given by
3
2.2 The incompressible Reynolds Averaged Navier Stokes (RANS)
The incompressible momentum equation and continuity equation respectively reads, [3] 𝜕𝑢𝑖 𝜕𝑡 + 𝜕 𝜕𝑥𝑗(𝑢𝑖𝑢𝑗) = − 1 𝜌 𝜕𝑝 𝜕𝑥𝑖+ 𝜕 𝜕𝑥𝑗[𝜈 ( 𝜕𝑢𝑖 𝜕𝑥𝑗+ 𝜕𝑢𝑗 𝜕𝑥𝑖)], (2.2.1) 𝜕𝑢𝑖 𝜕𝑥𝑖 = 0. (2.2.2)
Here 𝑢𝑖 is the velocity in the i:th direction, 𝜌 the density and 𝜈 is the kinematic viscosity of the fluid.
Here the tensor notation is used such that repeated indices indicate summation. Now, decompose each variable into a mean value, and a fluctuating value such that 𝑢𝑖= 𝑈𝑖+ 𝑢𝑖′ and 𝑝 = 𝑃 + 𝑝′ and
then take the time average of the momentum and continuity equations, remembering that density and viscosity are constants,
𝜕〈𝑈𝑖+𝑢𝑖′〉 𝜕𝑡 + 𝜕 𝜕𝑥𝑗〈(𝑈𝑗+ 𝑢𝑗 ′)(𝑈 𝑖+ 𝑢𝑖′)〉 = −1𝜌𝜕〈𝑃+𝑝 ′〉 𝜕𝑥𝑖 + 𝜕 𝜕𝑥𝑗[𝜈 ( 𝜕〈𝑈𝑖+𝑢𝑖′〉 𝜕𝑥𝑗 + 𝜕〈𝑈𝑗+𝑢𝑗′〉 𝜕𝑥𝑖 )], (2.2.3) 𝜕〈𝑈𝑖+ 𝑢𝑖′〉 𝜕𝑥𝑖 = 0. (2.2.4)
Note that for each term, the fluctuating part will be zero using (2.1.8). Aside from the convective term seen below, using (2.1.7) yield
𝜕 𝜕𝑥𝑗〈(𝑈𝑗+ 𝑢𝑗′)(𝑈𝑖+ 𝑢𝑖′)〉 = 𝜕 𝜕𝑥𝑗(𝑈𝑖𝑈𝑗) + 𝜕 𝜕𝑥𝑗〈𝑢𝑖′𝑢𝑗′〉. (2.2.5) As such the final form for the momentum and continuity equation respectively is
𝜕𝑈𝑖 𝜕𝑡 + 𝜕 𝜕𝑥𝑗(𝑈𝑖𝑈𝑗) = − 1 𝜌 𝜕𝑃 𝜕𝑥𝑖+ 𝜕 𝜕𝑥𝑗[𝜈 ( 𝜕𝑈𝑖 𝜕𝑥𝑗+ 𝜕𝑈𝑗 𝜕𝑥𝑖) − 〈𝑢𝑖′𝑢𝑗′〉], (2.2.6) 𝜕𝑈𝑗 𝜕𝑥𝑖 = 0. (2.2.7)
Here the term 〈𝑢𝑖′𝑢𝑗′〉 is commonly referred to as the Reynolds stress tensor. The problem to resolve
turbulent flows have now shifted into trying to model this tensor. Boussinesq hypothesized in 1877 that the turbulent stresses are proportional to the mean velocity gradient, with a closure coefficient conceptually known as the turbulent eddy viscosity, 𝜈𝑡. i.e. the Reynolds stresses are modelled as
−〈𝑢𝑖′𝑢𝑗′〉 = 𝜈𝑡[
𝜕𝑈𝑖 𝜕𝑥𝑗+
𝜕𝑈𝑗
𝜕𝑥𝑖], (2.2.8)
4 𝜕𝑈𝑖 𝜕𝑡 + 𝜕 𝜕𝑥𝑗(𝑈𝑖𝑈𝑗) = − 1 𝜌 𝜕𝑃 𝜕𝑥𝑖+ 𝜕 𝜕𝑥𝑗[(𝜈 + 𝜈𝑡) ( 𝜕𝑈𝑖 𝜕𝑥𝑗 + 𝜕𝑈𝑗 𝜕𝑥𝑖)]. (2.2.9) Like that, in order to resolve the turbulent eddies, the problem has been shifted into trying to model the turbulent eddy viscosity, for which many such turbulence models exist. The main problem is that there is no real universal turbulence model, but different models fit for different scenarios. This is indeed a problem, as it enforces a need of user knowledge of flow structures in order to determine whether a model is suitable.
2.3 Reynolds averaged scalar transport equation
The transport equation for a passive scalar c of concentration [kg m⁄ 3] for a neutral gas is [3] 𝜕𝑐 𝜕𝑡+ 𝜕 𝜕𝑥𝑖(𝑐𝑢𝑖) − 𝜕 𝜕𝑥𝑖(𝐷 𝜕𝑐 𝜕𝑥𝑖) = 𝑠. (2.3.1)
here s is the source term and D is the diffusion coefficient. As for the momentum equation, we decompose the variables into a mean value and a fluctuating value; 𝑢𝑖 = 𝑈𝑖+ 𝑢𝑖′, 𝑐 = 𝐶 + 𝑐′ and
𝑠 = 𝑆 + 𝑠′ then substitute into (2.3.1) and take the Reynolds time average as previously 𝜕〈𝐶 + 𝑐′〉 𝜕𝑡 + 𝜕 𝜕𝑥𝑖〈(𝐶 + 𝑐′)(𝑈𝑖+ 𝑢𝑖′)〉 − 𝜕 𝜕𝑥𝑖(𝐷 𝜕〈𝐶 + 𝑐′〉 𝜕𝑥𝑖 ) = 〈𝑆 + 𝑠′〉. (2.3.2) As before using (2.1.7) together with (2.1.8) we attain
𝜕𝐶 𝜕𝑡 + 𝜕 𝜕𝑥𝑖[𝐶𝑈𝑖] − 𝜕 𝜕𝑥𝑖(𝐷 𝜕𝐶 𝜕𝑥𝑖− 〈𝑐′𝑢𝑖′〉) = 𝑆. (2.3.3) Now, using a similar hypothesis as the one formulated by Boussinesq in (2.2.8), one may hypothesize that the average of the fluctuating terms is proportional to the gradient of the mean concentration, where the closing constant is referred to as the turbulent diffusion constant,
−〈𝑐′𝑢 𝑖 ′〉 = 𝐷 𝑡 𝜕𝐶 𝜕𝑥𝑖. (2.3.4)
Here the relation between the diffusivity and turbulent diffusivity to viscosity and turbulent viscosity is 𝐷 = 𝜈/𝑆𝑐 𝐷𝑡= 𝜈𝑡/𝑆𝑐𝑡 respectively, where 𝑆𝑐 is the Schmidt number, a fluid dependent number
which describes the ratio between viscosity and molecular diffusion. That is the ratio of the transfer of momentum (viscosity) and the mass transfer (molecular diffusion). 𝑆𝑐𝑡 Is then analogously the
turbulent Schmidt number, a model constant which needs to be determined. The final transport equation for the mean concentration C is then, together with the continuity equation
5
2.4 The Adjoint
Consider the real function v(x) with the following properties as stated by Marchuk [6]
a) v(x) is defined on the domain Ω and is continuous and differentiable in all interior points of the domain
b) v(x) is square summable on Ω together with the derivatives 𝑑𝑣 𝑑𝑥⁄ and 𝑑2𝑣 𝑑⁄ 2𝑥. that is, they exist and are finite;
∫ [(𝑑 2𝑣 𝑑𝑥2) 2 + (𝑑𝑣 𝑑𝑥) 2 + 𝑣2(𝑥)] 𝑑𝑥 𝛺 < ∞ (2.4.1)
Consider now the differential operator A on to the function v(x) equal to some function f, together with a boundary operator 𝐴𝑏 and boundary function 𝑓𝑏.
𝐴𝑣 = 𝑓, 𝑥 ∈ 𝛺 (2.4.2)
𝐴𝑏𝑣 = 𝑓𝑏, 𝑥 ∈ 𝜕𝛺 (2.4.3)
The functions v(x) with the property in (2.4.1) forms a subspace of the Hilbert space of square integrable functions. Introduce an inner product for the functions v(x) and w(x) in the Hilbert space H = L2(Ω) defined on Ω, denoted by
(𝑣, 𝑤) = ∫ 𝑣𝑤
𝛺 𝑑𝑥. (2.4.4)
Let now the operator A affect v(x), the result is a new function Av which is also defined on Ω and is in H. for this inner product, there exist an adjoint operator 𝐴∗ to the function w(x) for all v(x), w(x) satisfying a) and b) such that
(𝐴𝑣, 𝑤) = (𝑣, 𝐴∗𝑤). (2.4.5)
The equality in (2.4.5) is according to Marchuk [6] referred to as the Lagrange identity. Now assume that 𝐴∗𝑤 solves the equation below, with corresponding boundary conditions
𝐴∗𝑤 = 𝑓∗, 𝑥 ∈ 𝛺 (2.4.6)
𝐴𝑏∗𝑤 = 𝑓
𝑏∗, 𝑥 ∈ 𝜕𝛺. (2.4.7)
When applying the scalar product the attain that
(𝐴𝑣, 𝑤) = (𝑓, 𝑤) (2.4.8)
(𝐴∗𝑤, 𝑣) = (𝑓∗, 𝑣). (2.4.9)
Take equation (2.4.9)and subtract from(2.4.8);
6 According to (2.4.5) the left hand side is equal to zero; hence we get what Marchuk refers to as the adjoint duality representation
(𝑓, 𝑤) = (𝑓∗, 𝑣) (2.4.11)
2.5 Adjoint of Reynolds averaged transport equation
Denote the sum of diffusion and turbulent diffusion as 𝐾 = 𝜈 𝑆𝑐⁄ +𝜈𝑡 𝑆𝑐
𝑡
⁄ = 𝐷 + 𝐷𝑡, then equation (2.3.5) can be written as
𝐿𝐶 = 𝑆. (2.5.1)
Where the differential operator L is
𝐿 = 𝑑 𝑑𝑡+ 𝑈𝑖 𝜕 𝜕𝑥𝑖− 𝜕 𝜕𝑥𝑖𝐾 𝜕 𝜕𝑥𝑖 (2.5.2)
Now introduce a test function 𝐶∗ and define a scalar product between (2.5.1) and 𝐶∗ as
(𝐿𝐶, 𝐶∗) = ∫ 𝑑𝑡𝑇 0
∫ (𝐿𝐶)𝐶∗
𝛺 𝑑𝛺.
(2.5.3)
Then the scalar product of (2.5.1) and 𝐶∗ is
(𝐿𝐶, 𝐶∗) = (𝑆, 𝐶∗) (2.5.4) or explicitly ∫ 𝑑𝑡 𝑇 0 ∫ [ 𝜕𝐶 𝜕𝑡 + 𝑈𝑖 𝜕𝐶 𝜕𝑥𝑖− 𝜕 𝜕𝑥𝑖(𝐾 𝜕𝐶 𝜕𝑥𝑖)] 𝐶 ∗𝑑𝛺 𝛺 = ∫ 𝑑𝑡 𝑇 0 ∫ 𝑆𝐶 ∗𝑑𝛺 𝛺 . (2.5.5)
Since the product inside the integral is distributive we may write (2.5.5) as
∫ 𝑑𝑡𝑇 0 ∫ [𝐶 ∗𝜕𝐶 𝜕𝑡 + 𝐶∗𝑈𝑖 𝜕𝐶 𝜕𝑥𝑖− 𝐶∗ 𝜕 𝜕𝑥𝑖(𝐾 𝜕𝐶 𝜕𝑥𝑖)] 𝑑𝛺 𝛺 = ∫ 𝑑𝑡𝑇 0 ∫ 𝑆𝐶∗𝑑𝛺 𝛺 , (2.5.6)
7 I. ∫ ∫ 𝜕𝐶 𝜕𝑡𝐶∗𝑑𝑡 𝑇 0 𝑑𝛺 𝛺 = ∫ ∫ (𝜕 𝜕𝑡(𝐶∗𝐶) − 𝜕𝐶∗ 𝜕𝑡 𝐶) 𝑑𝑡 𝑇 0 𝑑𝛺 = 𝛺 = ∫ [𝐶∗𝐶] 0 𝑇𝑑𝛺 𝛺 − ∫ ∫ 𝜕𝐶 ∗ 𝜕𝑡 𝐶𝑑𝑡 𝑇 0 𝑑𝛺. 𝛺 (2.5.7) II. ∫ ∫ 𝑈𝑖 𝜕𝐶 𝜕𝑥𝑖𝐶∗𝑑𝛺 𝛺 𝑑𝑡 𝑇 0 = ∫ ∫ ( 𝜕 𝜕𝑥𝑖(𝑈𝑖𝐶𝐶∗) − 𝑈𝑖 𝜕𝐶∗ 𝜕𝑥𝑖𝐶) 𝛺 𝑑𝛺𝑑𝑡 𝑇 0 = = ∫ ∫ 𝑈𝑖𝐶𝐶∗𝑑𝜕𝛺 𝜕𝛺 𝑑𝑡 𝑇 0 − ∫ ∫ 𝑈𝑖𝜕𝐶∗ 𝜕𝑥𝑖 𝐶𝑑𝛺 𝛺 𝑑𝑡. 𝑇 0 (2.5.8) III. − ∫ ∫ 𝜕 𝜕𝑥𝑖(𝐾 𝜕𝐶 𝜕𝑥𝑖) 𝐶∗𝑑𝛺 𝛺 𝑑𝑡 𝑇 0 = − ∫ ∫ ( 𝜕 𝜕𝑥𝑖(𝐾 𝜕𝐶 𝜕𝑥𝑖𝐶 ∗) − 𝐾𝜕𝐶 𝜕𝑥𝑖 𝜕𝐶∗ 𝜕𝑥𝑖 𝑑Ω) Ω 𝑑𝑡 𝑇 0 = − ∫ ∫ [𝜕𝑥𝜕 𝑖(𝐾 𝜕𝐶 𝜕𝑥𝑖𝐶∗) − ( 𝜕 𝜕𝑥𝑖(𝐾 𝜕𝐶∗ 𝜕𝑥𝑖 𝐶) − 𝜕 𝜕𝑥𝑖(𝐾 𝜕𝐶∗ 𝜕𝑥𝑖) 𝐶) 𝑑Ω] Ω 𝑑𝑡 𝑇 0 = − ∫ ∫ [𝜕𝑥𝜕 𝑖(𝐾 𝜕𝐶 𝜕𝑥𝑖𝐶∗− 𝐾 𝜕𝐶∗ 𝜕𝑥𝑖 𝐶) + 𝜕 𝜕𝑥𝑖(𝐾 𝜕𝐶∗ 𝜕𝑥𝑖) 𝐶] 𝑑Ω Ω 𝑑𝑡 𝑇 0 = − ∫ 𝑑𝑡𝑇 0 ∫ 𝜕 𝜕𝑥𝑖(𝐾 𝜕𝐶∗ 𝜕𝑥𝑖) 𝐶 𝑑𝛺 𝛺 + ∫ 𝑑𝑡𝑇 0 ∫ (𝐾𝜕𝐶∗ 𝜕𝑥𝑖𝐶 − 𝐾 𝜕𝐶 𝜕𝑥𝑖𝐶∗) 𝑑𝜕𝛺 𝜕𝛺 . (2.5.9)
8 Assume now that 𝐶∗(Ω, 𝑡) is a solution to
𝐿∗𝐶∗= 𝑆∗ (2.5.13)
with the adjoint operator 𝐿∗ defined as
𝐿∗= − 𝜕 𝜕𝑡− 𝑈𝑖 𝜕 𝜕𝑥𝑖− 𝜕 𝜕𝑥𝑖(𝐾 𝜕 𝜕𝑥𝑖) (2.5.14)
together with the boundary conditions
𝐿∗𝑏𝐶∗= 𝑆
𝑏∗, (2.5.15)
where 𝐿∗𝑏 is the boundary operator, which will be defined such that (2.5.12) is equal to zero. Then,
for the scalar product (𝐿𝐶, 𝐶∗) there exist the adjoint operator 𝐿∗ such that (𝐿𝐶, 𝐶∗) + (𝐿
𝑏𝐶, 𝐶∗)𝑏= (𝐶, 𝐿∗𝐶∗) + (𝐶, 𝐿∗𝑏𝐶∗)𝑏. (2.5.16)
And similar to what was attained in (2.4.11) there exist an adjoint duality representation of (2.5.16) (𝑆, 𝐶∗) + (𝐿
𝑏𝐶, 𝐶∗)𝑏 = (𝐶, 𝑆∗) + (𝐶, 𝐿∗𝑏𝐶∗)𝑏. (2.5.17)
The scalar product (. , . )𝑏denotes the surface integral over the boundary, involving the boundary
conditions. The advantage of having this view on the problem as opposed to how Marchuk, is that this allows for a source on the boundary, as opposed to only allowing it to exist within the domain.
2.6 A Measurement and the adjoint duality
Consider a measurement 〈𝑐𝑗〉 from a set of measurements of size n, 𝑗 ∈ [1, 𝑛]. In general the
measurement will be a time and spaced average. If the measurement can be considered a point measure in space, i.e. the measurement 〈𝑐𝑗〉 occurs in a point Ω𝑗∈ Ω, corresponding to a spatially
perfect measurement, Define the measurement as
〈𝑐𝑗〉 = 1 𝑇∫ 𝑑𝑡 𝑇 0 ∫ 𝑐(𝛺𝑗)𝑑𝛺 𝛺 =1 𝑇∫ 𝑑𝑡 𝑇 0 ∫ 𝑐𝛿(𝛺 − 𝛺𝑗)𝑑𝛺 𝛺 = 𝑐(𝛺𝑗). (2.6.1)
Now, let 𝑆𝑗∗: = 𝛿(Ω − Ω𝑗) in the adjoint equation in (2.5.13), then the right hand part of the scalar
product from (2.5.17) reads
(𝐶, 𝑆𝑗∗) = ∫ 𝑑𝑡 𝑇 0 ∫ 𝐶𝛿(𝛺 − 𝛺𝑗)𝑑𝛺 𝛺 = 𝑇𝐶(𝛺𝑗) = 𝑇〈𝑐𝑗〉. (2.6.2)
Using the adjoint duality (𝑆, 𝐶𝑗∗) = (𝐶, 𝑆𝑗∗) + (𝐶, 𝐿∗𝑏𝐶𝑗∗)𝑏− (𝐿𝑏𝐶, 𝐶𝑗∗)𝑏 a possibility arises to locate
the source S, by solving the adjoint transport equation for the scalar field 𝐶𝑗∗ for each measurement
9 assumption of a point measure means that any spatial error is assumed to be zero, i.e. that the position of the measurement is indeed a fixed point.
2.7 Source strength
Not only is the position of the source a problem itself, but the strength of the source must be taken care of. Consider the source 𝑆 = 𝑄𝑆̂, where 𝑄[𝑘𝑔] is the strength of the source and 𝑆̂[𝑚−3𝑠−1] is a normalized source density, describing both the shape and the position of the source. This shape can be general; it may be a point source in time and space (an explosion), a continuous in time point in space source, a surface source on a boundary or even a volume source inside the domain. The scalar product in (2.5.17) will then become
(𝐶𝑗∗, 𝑆) = (𝐶
𝑗∗, 𝑄𝑆̂) = 𝑄(𝐶𝑗∗, 𝑆̂), (2.7.1)
as such the adjoint duality will be
𝑄(𝐶𝑗∗, 𝑆̂) = (𝐶, 𝑆
𝑗∗) + (𝐶, 𝐿∗𝑏𝐶𝑗∗)𝑏− (𝐿𝑏𝐶, 𝐶𝑗∗)𝑏. (2.7.2)
10
2.8 The scalar product in a computational grid
There are many different ways of viewing the scalar product defined in (2.7.1), for instance [7] views it as a optimization over a set of measures, what we will do is view it as a Bayesian method of reconstructing the source, namely to view it as a probability distribution. Consider the domain Ω discretized into m grid cells such that Ω ≈⋃𝑚 Ω𝑖
𝑖=1 . Since the source may have any shape inside the
domain, be it a line source, area source, volume source or even a point source. Consider the source or as an instantaneous source in space, i.e. a Dirac function in space
𝑆̂ = 𝛿(Ω − Ω𝑘), (2.8.1)
or a boundary source
𝑆̂𝑏 = 𝛿(Ωb− Ωb𝑘), (2.8.2)
where Ω𝑏𝑘 is a point at the boundary, as opposed to Ωk which is an interior point.
Now the scalar product reads
(𝐶𝑗∗, 𝑆̂) = ∫ 𝑑𝑡 𝑇 0 ∫ 𝐶𝑗∗𝛿(𝛺 − 𝛺 𝑘)𝑑𝛺 𝛺 = ∫ 𝑑𝑡 𝑇 0 𝐶𝑗∗(𝛺 𝑘) = 𝑇𝐶𝑗∗(𝛺𝑘), (2.8.3) and similarly (𝐶𝑗∗, 𝑆̂ 𝑏)𝑏= ∫ 𝑑𝑡 𝑇 0 ∫ 𝐶𝑗∗𝛿(𝛺 𝑏− 𝛺𝑏𝑘)𝑑𝜕𝛺 𝛺 = ∫ 𝑑𝑡 𝑇 0 𝐶𝑗∗(𝛺 𝑏𝑘) = 𝑇𝐶𝑗∗(𝛺𝑏𝑘). (2.8.4)
That is, the scalar product is essentially the value that the adjoint plume has for each point Ω𝑘 in
space. The last equality holds as 𝐶𝑗∗, the Reynolds averaged concentration has no time dependence.
The adjoint duality finally reads
𝑇𝑄[𝐶𝑗∗(𝛺𝑘) + 𝐶𝑗∗(𝛺𝑏𝑘)] = 𝑇〈𝑐𝑗〉 (2.8.5)
and since Ωb𝑘 ∈ Ω𝑘
𝑄𝐶𝑗∗(𝛺
𝑘) = 〈𝑐𝑗〉. (2.8.6)
11
2.9 Source term estimator
As done by [8] denote the a posteriori probability 𝑃(𝑆̂Ω𝑘) of the source by
𝑃(𝑆̂Ω𝑘) ∝ 𝑒𝑥𝑝 (− ∑ [𝑄𝐶𝑗 ∗(𝛺 𝑘) − 〈𝑐𝑗〉 𝜎𝑗(𝑄𝐶𝑗∗(𝛺𝑘)) ] 2 𝑛 𝑗=0 ) . (2.9.1)
That is, the probability is proportional to the exponential of the negative least squares between the measured concentration and the adjoint values at each discrete point. Here the standard deviation is approximated as σj= k0QCj∗(Ωk), where k0 is a fixed constant. However first, an optimal Q = Q̂
should be estimated. This can be attained by minimizing the least square in (2.9.1) as
𝑄̂ = 𝑚𝑖𝑛 𝑄 ∑ [ 𝑄𝐶𝑗∗(𝛺𝑘) − 〈𝑐𝑗〉 𝑘0𝑄𝐶𝑗∗(𝛺 𝑘) ] 2 𝑛 𝑗=0 . (2.9.2)
The solution to this one-dimensional quadratic minimization can be solved explicitly as
𝑄̂ =∑ 〈𝑐𝑗〉 2 𝐶 𝑗∗2(Ω𝑘) ⁄ 𝑛 𝑗=0 ∑ 〈𝑐𝑗〉 𝐶𝑗∗(Ω 𝑘) ⁄ 𝑛 𝑗=0 . (2.9.3)
To find the normalized source term ŜΩk, we solve the minimization problem min𝑆̂Ω𝑘𝑃(𝑆̂Ω𝑘) for
k = 1 … n where 𝑃(𝑆̂𝛺𝑘) = 𝑒𝑥𝑝 (− ∑ [𝑄̂𝐶𝑗 ∗(𝛺 𝑘) − 〈𝑐𝑗〉 𝑘0𝑄̂𝐶𝑗∗(𝛺 𝑘) ] 2 𝑛 𝑗=0 ) 𝜉 , (2.9.4)
and the normalization factor 𝜉 is given by
𝜉 = ∑ 𝑃(𝑆̂𝛺𝑘)
𝑚 𝑘=1
. (2.9.5)
12
2.10 The Atmospheric Boundary Layer (ABL)
Due to friction that the ground exerts on the air, the velocity of the air will get lower the closer to the ground one gets, much like the boundary layer close to walls. As such a velocity profile will form which is identical in profile as the viscous boundary layers close to walls. Together with a no slip condition, i.e. The velocity at the boundary of the ground is zero. Hence it is possible to model this boundary layer utilizing the logarithmic function seen below
𝑈𝑥(𝑧) 𝑈𝑟𝑒𝑓 = 𝑢∗ 𝜅 𝑙𝑛 𝑧 + 𝑧0 𝑧0 . (2.10.1)
Here 𝑢∗is the friction velocity, a measurement for the turbulent shear stress, 𝜅 the von Karman
constant, with the value of 0.41, z the distance from the ground and 𝑧0 is the roughness length. The
roughness length is a measure for the roughness of the ground, usually some distinct length of the region in question, if the region is a forest, the average length of the trees may be used as the roughness length.
2.11 K-ε Model for the atmospheric boundary layer
In order to model the turbulent viscosity seen in the Reynolds Averaged Navier Stokes equation, Launder & Spalding [9] postulated in 1972 that it can be modelled using
𝜈𝑡 = 𝐶𝜇𝑘2
𝜖. (2.11.1)
Where k is the turbulent kinetic energy, defined as 𝑘 = 0.5 ∙ [〈𝑢1′ 2〉 + 〈𝑢2′ 2〉 + 〈𝑢3′ 2〉] and ε is the
turbulent dissipation rate, i.e. the rate at which the kinetic energy dissipates to heat. In order to model these two quantities, two transport equations are introduced for each of these quantities
𝜕𝑘 𝜕𝑡+ 𝑈𝑖 𝜕𝑘 𝜕𝑥𝑖 = 𝜕 𝜕𝑥𝑗[(𝜈 + 𝜈𝑡 𝜎𝑘) 𝜕𝑘 𝜕𝑥𝑗 ] + 𝑃𝑘− 𝜖, (2.11.2) 𝜕𝜖 𝜕𝑡+ 𝑈𝑖 𝜕𝜖 𝜕𝑥𝑖 = 𝜕 𝜕𝑥𝑗[(𝜈 + 𝜈𝑡 𝜎𝜖) 𝜕𝜖 𝜕𝑥𝑗 ] + 𝐶1 𝜖 𝑘𝑃𝑘− 𝐶2 𝜖2 𝑘. (2.11.3)
Or in words, the transport of turbulent kinetic energy and dissipation rate can be translated to 𝑟𝑎𝑡𝑒 𝑜𝑓 𝑐ℎ𝑎𝑛𝑔𝑒 + 𝑡𝑟𝑎𝑛𝑠𝑝𝑜𝑟𝑡 𝑏𝑦 𝑎𝑑𝑣𝑒𝑐𝑡𝑖𝑜𝑛 = 𝑡𝑟𝑎𝑛𝑠𝑝𝑜𝑟𝑡 𝑏𝑦 𝑑𝑖𝑓𝑓𝑢𝑠𝑖𝑜𝑛 + 𝑆𝑜𝑢𝑟𝑐𝑒 − 𝑆𝑖𝑛𝑘, Where the production term 𝑃𝑘 is modelled by
𝑃𝑘 = 2𝜈𝑡𝑆𝑖𝑗𝑆𝑖𝑗, (2.11.4) 𝑆𝑖𝑗=1 2( 𝜕𝑈𝑖 𝜕𝑥𝑗+ 𝜕𝑈𝑗 𝜕𝑥𝑖). (2.11.5)
13
Table 1: Coefficients for the modified ABL k-ε equation
σk σϵ C1 C2 Cμ
1.0 1.3 1.176 1.92 √0.09
In addition to this, a model hereafter referred to as the MMK model, [11] proposed a modification to the constant 𝐶𝜇 in order to describe the recirculation vortices behind blunt bodies which proved to
yield positive results as opposed to the standard K-ε model. The modified constant 𝐶𝜇∗ is
𝐶𝜇∗= {𝐶𝜇, Ω ̃ 𝑆̃ ⁄ > 1 Ω̃ 𝑆̃ ⁄ , Ω̃⁄ < 1,𝑆̃ (2.11.6)
Here Ω̃ is the vorticity and 𝑆̃ is the strain rate, defined respectively as
Ω̃ = (𝜕𝑈𝑖 𝜕𝑥𝑗− 𝜕𝑈𝑗 𝜕𝑥𝑖), (2.11.7) 𝑆̃ = (𝜕𝑈𝜕𝑥𝑖 𝑗+ 𝜕𝑈𝑗 𝜕𝑥𝑖). (2.11.8)
14
2.12 Boundary conditions
The inlet boundary conditions for the atmospheric boundary layer using the k-ε model should according to [12] be set as 𝑈𝑥=𝑢∗ 𝜅 𝑙𝑛 ( 𝑧 + 𝑧0 𝑧0 ), (2.12.1) 𝑘 = 𝑢∗2 √𝐶𝜇 , (2.12.2) 𝜖 = 𝑢∗3 𝜅(𝑧 + 𝑧0). (2.12.3)
However, since there are data accessible for the kinetic energy k for the wind tunnel, instead these values will be used in the simulations. Furthermore for the ground and array, the turbulent viscosity is modelled using a standard wall function but with the roughness length 𝑧0 for the ground instead of
the standard E=9.8, used for walls – as seen below. 𝑦+=𝑢∗𝑦 𝜈 , 𝐸̅ = 𝑦 + 𝑧0 𝑧0 (2.12.4) 𝜈𝑡𝑔𝑟𝑜𝑢𝑛𝑑= 𝜈 ( 𝑦 +𝜅 𝑙𝑜𝑔(𝐸̅)− 1) , 𝜈𝑡𝑤𝑎𝑙𝑙= 𝜈 ( 𝑦+𝜅 𝑙𝑜𝑔(9.8)− 1) (2.12.5) Here y is the distance to the wall, or to the ground depending on the boundary. 𝑦+ is the dimensionless distance to the wall, as previously 𝜅 is the von Karman constant and 𝑢∗ is the friction
velocity. The full list of boundary conditions applied is seen in table 2 below.
Table 2: Boundary conditions for the Adjoint Reynolds Averaged Scalar Transport with modified k-ε model.
Inlet Outlet Left / Right
wall
Top Ground Array
U (2.12.1) 𝒏 ∙ 𝛻𝑼 = 0 Slip Slip No slip No slip
P 𝒏 ∙ 𝛻𝑃 = 0 𝑃 = 0 Slip Slip 𝒏 ∙ 𝛻𝑃 = 0 𝒏 ∙ 𝛻𝑃 = 0
K Exp. values 𝒏 ∙ 𝛻𝑘 = 0 Slip Slip Wall function Wall function
ε (2.12.3) 𝒏 ∙ 𝛻𝜖 = 0 Slip Slip Wall function Wall function
𝐂𝒊∗ 𝒏 ∙ 𝛻𝐶
15
2.13 The Turbulent Schmidt number
There have been many studies revolving the turbulent Schmidt number 𝑆𝑐𝑡 = 𝜈𝑇⁄𝐷𝑇, ever since
Spalding and Launder reported values for 𝑆𝑐𝑡 = 0.7 in 1971 and 𝑆𝑐𝑡 = 0.9 in 1978 respectively.
These values have been commonly used in CFD studies for turbulent mass diffusion, even going as far as to putting these values as default in commercial CFD software [13]. Have been researching different flows in order to determine whether it is reasonable to assume the turbulent Schmidt number as either of the values specified by Launder and Sharma. His results show that it is highly flow dependent and that a value of 𝑆𝑐𝑡 = 0.3 for dispersion around a building shows the best
agreement with experimental data.
Further on, [14] has made experiments in a wind tunnel, which mimics the atmospheric boundary layer, and compiled own results with results previously found from other experiments. These results shows a strong height dependency on the Schmidt number, and constructed a power series which approximated the Schmidt number as
𝑆𝑐𝑡(𝑧) = ∑ 𝑎𝑖( 𝑧 ℎ𝑟𝑒𝑓) 𝑖 , 5 𝑖=0 (2.13.1) 𝑎 = (−0.226, 12.2, −46.2, 81.0, −67.9, 21.5).
The graph for this can be seen below
Figure 1: The power series for the turbulent Schmidt number height as described by [15]
16
2.14 Formal definition of the Adjoint Equation
Assume that the main equation is 𝜕𝐶 𝜕𝑡 + 𝑈𝑖 𝜕𝐶 𝜕𝑥𝑖− 𝜕 𝜕𝑥𝑖(𝐾 𝜕𝐶 𝜕𝑥𝑖) = 𝑆, 𝑥𝑖 ∈ Ω, t ∈ [0, 𝑇]. (2.14.1) With boundary conditions
(𝐾𝜕𝐶
𝜕𝑥𝑖) = 𝑆𝑏, 𝑥𝑖 ∈ 𝜕Ω, t ∈ [0, 𝑇], (2.14.2)
𝐶(𝑥𝑖, 0) = 0, 𝑥𝑖 ∈ Ω, (2.14.3)
for a boundary source. Otherwise the boundary conditions is 𝜕𝐶
𝜕𝑥𝑖 = 0, 𝑥𝑖 ∈ 𝜕Ω, t ∈ [0, 𝑇]. (2.14.4) Here the release of the scalar C and the measurement 〈𝐶𝑗〉 is inside of the domain Ω and is within the
time interval [0, 𝑇]. This means that S(xi, 𝑡) = S∗(xi, 𝑡) = 0 for 𝑡 ∉ [0, 𝑇]. The adjoint problem will
then together with boundary conditions be
−𝜕𝐶∗ 𝜕𝑡 − 𝑈𝑖 𝜕𝐶∗ 𝜕𝑥𝑖 − 𝜕 𝜕𝑥𝑖(𝐾 𝜕𝐶∗ 𝜕𝑥𝑖) = 𝑆∗, 𝑥𝑖 ∈ Ω, t ∈ [𝑇, 0] (2.14.5) 𝜕𝐶∗ 𝜕𝑥𝑖 = 0, 𝑥𝑖 ∈ 𝜕Ω, t ∈ [𝑇, 0] (2.14.6) 𝐶∗(𝑥 𝑖, 𝑇) = 0, 𝑥𝑖 ∈ Ω (2.14.7)
Translate the time as t1= 𝑇 − 𝑡 and the adjoint in (2.14.5) will be much more similar to (2.14.1)
𝜕𝐶∗ 𝜕𝑡1 − 𝑈𝑖 𝜕𝐶∗ 𝜕𝑥𝑖 − 𝜕 𝜕𝑥𝑖(𝐾 𝜕𝐶∗ 𝜕𝑥𝑖) = 𝑆∗, 𝑥𝑖 ∈ Ω, t1 ∈ [0, 𝑇] (2.14.8) 𝜕𝐶∗ 𝜕𝑥𝑖 = 0, 𝑥𝑖 ∈ 𝜕Ω, t1 ∈ [0, 𝑇] (2.14.9) 𝐶∗(𝑥 𝑖, 0) = 0, 𝑥𝑖 ∈ Ω (2.14.10)
The adjoint duality representation will now read, together with (2.8.1) and (2.6.2) 𝑄̂𝐶𝑗∗(Ω
17
3 Outline of numerical solution in OpenFOAM
The platform of implementation was chosen to be OpenFOAM, an open source software written in C++ which utilizes the finite volume method to solve partial differential equations. FOAM stands for Field Operation And Manipulation, and the implementation of the methods described in this work has their starting point in the already implemented kEpsilon model along with ScalarTransportFoam – the standard k-ε model and the transport equation of a passive scalar, below the only the interesting snippets of implemented code is shown.
3.1 Implementation of the Adjoint equation
The first thing one notices when working with OpenFOAM is the ease to implement the equations to be solved for. As an example, the adjoint equation in (2.14.12) would be written as a finite volume scalar matrix, which in FOAM code would be
fvScalarMatrix cstarEQ ( fvm::ddt(*c_star) +fvm::div(-phi, *c_star) -fvm::laplacian(K_, *c_star) == fvOptions(*c_star) ); cstarEQ.solve();
Where the fvm class implies that a finite volume matrix will be built, FvOptions being the method in which the point source is defined, this being documented and described further in the OpenFOAM manual [16]. The equation is then solved using the solve() class. Please note the pointer, *c_star. The reason for this as opposed to how it normally would be implemented is that since one adjoint equation has to be solved for each measurement, the code need to be dynamic such that the amount of measurements may vary. Luckily since each adjoint equation has the same boundary conditions, the code implemented reads a “dummy” file with all the boundary conditions, and then creates one scalar field c_star for each measurement supplied by the user, as seen in the code snippet below.
18 volScalarField* fields[c[0].size()];
for(int i=0; i < c[0].size(); i++) {
char fieldname[sizeof(int)*3+6];
snprintf(fieldname, sizeof fieldname, "c_star%d", i);
volScalarField* c_starttmp = new volScalarField
( IOobject ( fieldname, runTime.timeName(), mesh, IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE ), c_star ); fields[i] = c_starttmp; }; // Courtesy of H.Grahn @ FOI
Now for each adjoint equation to be solved, storing the equation from above in a separate file cstarEq.H, each adjoint concentration may be solved by looping over the pointer array fields.
for (int i = 0; i < c[0].size(); i++) {
volScalarField* c_star = fields[i];
#include "cstarEqn.H"
}
3.2 Numerical Schemes and Algorithms
OpenFOAM comes packed with different Numerical Schemes, and more often than not will it be hard to determine which to use without having to dwell deep into the code and theory behind matrix inversion. It is generally recommended to view the tutorials that come packed with OpenFOAM, find something similar to what you want to solve and use the same Schemes and Algorithms.
For the RANS model, the pressure is solved using a Generalized Algebraic Multi Grid solver, which in short means that the pressure is solved in a fine grid and a coarse grid, and then being mapped between the two. This speeds up convergence. The velocity field, k and ε are respectively solved using the smoothSolver. All of which are converged down to a residual of 10−6, at which the flow can be assumed as steady state.
Each of the adjoint equations 𝐶𝑗∗, 𝑗 ∈ 0, . . , 𝑛 is uses a preconditioned bi-conjugate gradient solver and
19
3.3 Pre-processing - Parallel computing
The ability to run in Parallel exists within OpenFOAM, which shows great scaling. By doubling the amount of processors the solution time is roughly halved, depending on the whether the grid is coarse or fine enough. If it is too coarse the scaling will decline – OpenFOAM spends too much time communicating between the processors instead of spending time to solve the equations. The solutions were performed on a dedicated computational server with 64 processors.
3.4 Post processing
Unfortunately OpenFOAM does not come with any graphical visualization tool; it is only text based and as such run via the terminal. Luckily there is open source software available for visualization. The software used for that was Paraview, which can be seen in the window below.
20
4 The domain
The wind tunnel in question is the Environmental Flow (EnFlo) Wind tunnel in the MOdelling the Dispersion of Toxic Industrial Chemicals in urban environments (MODITIC) international project. This tunnel is built to mimic the atmospheric boundary layer in order to study flow and dispersion around objects. The computational domain was chosen to be 3 meters wide, 4 meters long and 1 meter high, that is −2 ≤ 𝑥 ≤ 2, −1.75 ≤ 𝑦 ≤ 1.75 and 0 ≤ 𝑧 ≤ 1.
4.1 The simple array
No concentration measurements of the dispersed gas were made on the atmospheric boundary layer exclusively, but for a small amount of 4 blocks of length h=0.109m in a 2x2 layout. This is hereby referred to as the simple array, and is meant to be the first step of an increasingly intricate model. The figure below describes the computational domain, along with the true source position and the simple array.
Figure 3: not to scale sketch of the computational domain for the simulation, complete with coordinate system and the simple array.
4.2 Mesh
The mesh is first built via the blockMesh tool in OpenFOAM, where the domain first is created; the mesh is then altered as more features are added to the domain. First the domain is built in two parts, one upper part and one lower part, where the lower part is where everything of interest is taking place. The height of the domain is as previously stated 1 meter, the lower part of the domain was taken as when z ∈ [0,2h] ≈ [0 , 0.218] m and the upper part was taken as when z ∈ (2h, 1] m. So the upper part of the domain was discretized into 80x80x15 cells, whereas the lower part was discretized into 80x80x25 cells for a total of 320 000 cells. Further refinements was then made for y ∈ [−0.75, 0.75] throughout the domain using the refineMesh tool, refining in x and y direction by splitting the cells in half. The total number of cells before the introduction of the simple array is 571 372 cells. This is the domain for which the verifications are run, to validate the turbulence model for the wind tunnel data.
21 used by introducing a box around the simple array and enhancing the mesh by once more splitting the cells in half, this time in x, y and z. The final number of cells ended up at 1 169 050 cells and can be seen in the figures below. The reason for such a high number of cells is mainly due to the dedicated numerical scientific computing server available.
22
5 Results
First off, in order to for the inverse model to make sense, one will have to determine how well the forward model describes the dispersion of a neutral gas. Hence measurements have been performed at the EnFlo wind tunnel, to determine the neutral boundary layer, and fitting it to the logarithmic equation stated in (2.12.1). From [17] Robins have estimated the friction velocity 𝑢∗= 0.055 and the surface roughness length 𝑧0= 0.0008 with a roughness Reynolds number
𝑢∗𝑧 0
𝜈 = 2.9
This in turn tells that the kinematic viscosity of the air in the wind tunnel is 𝜈 = 1.52 ∙ 10−5[𝑚𝑠2]. Also the reference height and reference velocity are respectively ℎ𝑟𝑒𝑓 = 1[𝑚], 𝑈𝑟𝑒𝑓 = 1.068. These are
the numbers for which the simulations are based on. Furthermore the Schmidt number for air is set to 0.7.
Below is a short flow scheme of how the solution process takes place.
23
5.1 Wind profile and turbulence – mimicking EnFlo
Applying the boundary conditions stated in section 2.12, the converged solution of the modified ABL k- ε model in an empty domain yield the following velocity profile and kinetic energy profile respectively.
Figure 6: Dimensionless velocity profile for the inlet and outlet of the simulated domain, together with the measured values from the EnFlo wind tunnel.
24
25
5.2 Flow over the simple array
Under the assumption that the turbulence model describes the atmospheric boundary layer in the EnFlo wind tunnel, the following flow can be found around the simple array through different planes over the array.
Figure 9: Velocity field over the simple array for one set of the cubes, taken at y = -h
26
Figure 11: Velocity field through the center of the simple array, taken at y = 0.0.
27
Figure 13: Velocity field around the simple array, taken at z = 0.1h.
28
Figure 15: Velocity field around the simple array, taken at z = 0.5h.
29
Figure 17: Velocity field around the simple array, taken at z = 1.1h.
30
5.3 Measurement positions
The data from the tunnel has been measured at key areas of interest, as can be seen in figure 19.
Figure 19: The measurement positions from which data from EnFlo is available
As such, a forward model is tested with different turbulent Schmidt numbers, in order to determine which of value to use. The resulting graphs can be seen in A.1 and are rather interesting in that the commonly used value of 0.7 does not model the dispersion of the neutral gas in a satisfying way. Analysing the graphs, given the forward model is only made to be a quantitative ground for which Schmidt number to use, it can be seen that the value of 0.7 under predicts the turbulent mass exchange. The less commonly used value of 0.3 and the z-dependency yield roughly the same behaviour of the released gas.
With this information in mind, it is time to determine how well the adjoint, or backward model is in determining the source term. This will be performed in two steps; the first step is by applying the model to synthetic data, taken from the forward run with 𝑆𝑐 = 0.3, measurement positions that have been hypothesized as probable point measurements were the simple array to be viewed as some sort of facility which performs live data collection of the air around it. See figure 19 below.
Figure 20: Points from which a real compound could potentially have measurement stations
The synthetic data and the data from EnFlo at the corresponding points can be seen in table 3 below.
Table 3: Synthetic data generated with 𝑺𝒄𝒕= 𝟎. 𝟑 and corresponding measurements from EnFlo.
Point X[m] Y[m] Z[m] c̅[kgm−3] synthetic c̅[kgm−3] EnFlo
1 0.75 0.0 0.025 0.1163 0.0978
2 1.036 -0.267 0.025 0.0096 0.0395
3 1.036 0.0 0.025 0.0718 0.0629
4 1.036 0.267 0.025 0.0159 0.0185
31
5.4 Synthetic data test
Figure 21: Probability distribution of the source for the synthetic data given from table 3. The blue level surface plotted here is for which the probability is within 𝟏𝟎−𝟑𝟎, a number chosen in the code. Note how the shape is stretched out in the wind direction (along the x-axis) and as such the high uncertainty along the wind direction.
32
Figure 23: Predicted source strength for the synthetic data throughout the domain. Note the high variability of source strength throughout the domain that differs by an order of magnitude of 𝟏𝟎𝟏𝟕. The high values correspond to the source being so strong that it would be transported simply by diffusion to the sensors. This being unrealistic, the blue region is the interesting bit.
33
5.5 ENFLO simple array data
Before going on to the real measured data from the wind tunnel, analysing the differences between the best agreeing synthetic model data (that is, 𝑆𝑐𝑡 = 0.3) with that of the tunnel. Note the big error
primarily in point 2 (see figure 20 and table 3). One could suspect that the two points adjacent to the array in the y axis, i.e. points 2 and 4 and their error could inflict an equal error in the adjoint model. As such tests with and without for the two points will be made, mainly using the turbulent Schmidt number of 0.3.
5.5.1 Five measurements
Figure 25: Probability distribution for data from the EnFlo wind tunnel where 𝑺𝒄𝒕= 𝟎. 𝟑and five measurements. The blue
level surface plotted here is for which the probability is within 𝟏𝟎−𝟑𝟎, as chosen in the code. Notice the similar shape as in figure 21, i.e. a high certainty along the y-axis but high uncertainty along the wind direction, but that it being further upstream that the real source.
Figure 26: Probability distribution for data from the EnFlo wind tunnel where 𝑺𝒄𝒕= 𝑺𝒄𝒕(𝒛) and five measurements. The
34
5.5.2 Three measurements
Figure 27: Probability distribution for data from the EnFlo wind tunnel where 𝑺𝒄𝒕= 𝟎. 𝟑 but with three measurements.
The blue level surface plotted here is for which the probability is within 𝟏𝟎−𝟑𝟎, as chosen in the code. Notice how the shape changed, with the level surface gaining increased uncertainty in both the wind direction and the y-axis. However the ratio between them is roughly the same, with a higher uncertainty along the wind direction. It does appear to swallow both the distributions in figures 21 and 25.
Figure 28: High posteriori distribution (HPD) for which 𝑷(𝑺̂𝜴𝒌) = 𝟎. 𝟗 ∙ 𝑷(𝑺̂𝜴𝒌)𝒎𝒂𝒙 for data from the EnFlo wind tunnel
where 𝑺𝒄𝒕= 𝟎. 𝟑 and three measurements. The blue level surface plotted here is for which the probability is
35
Figure 29: Least square estimated source strength throughout the domain with three measurements and 𝑺𝒄𝒕= 𝟎. 𝟑.
Notice how the shape is still as that of the synthetic data seen in figure 23, but with a different ratio of the order of magnitude between the maximum and minimum values. The high values correspond to the source being so strong that it would be transported simply by diffusion to the sensors. This being unrealistic, the blue region is the interesting bit.
Figure 30: Source strength estimate scaled down to manageable size for three measurements and 𝑺𝒄𝒕= 𝟎. 𝟑. This shape
36
Figure 31: Probability distribution for three measurements and 𝑺𝒄𝒕= 𝑺𝒄𝒕(𝒛). The blue level surface plotted here is for
which the probability is within 𝟏𝟎−𝟑𝟎, a number chosen in the code. Notice that the shape is different as for that seen in figure 26 but resembles figure 27 in shape. The same is applies here, there is a high uncertainty along the wind direction but a high certainty in the y-axis.
Figure 32: High posteriori distribution (HPD) for which 𝑷(𝑺̂𝜴𝒌) = 𝟎. 𝟗 ∙ 𝑷(𝑺̂𝜴𝒌)𝒎𝒂𝒙for three measurements and 𝑺𝒄𝒕= 𝑺𝒄𝒕(𝒛). Unlike the synthetic data in figure 21 and for 𝑺𝒄𝒕= 𝟎. 𝟑 in figure 28 the HPD is very close to the first
37
Table 4: Comparisons between the estimated source strength for the origin (true source location) and the High posteriori distribution regions for different cases. The high posteriori distribution (HPD) is for when 𝑷(𝑺̂𝜴𝒌) = 𝟎. 𝟗 ∙ 𝑷(𝑺̂𝜴𝒌)𝒎𝒂𝒙and
the origin is where we expect to find the source and as such the values shown here is the estimated source strength where we know the source is. The real source strength is 𝑸 = 𝟏. 𝟏 ∙ 𝟏𝟎−𝟑𝒌𝒈/𝒔.
38
6 Discussion
First off ponder back to the velocity profile given in figure 6, along with the kinetic energy in figure 7. For the velocity profile the first noticeable trait is that the profile for the inlet and outlet of the domain is rather alike for all heights, with some negligible differences. But what is most pleasant is the fact that both the inlet and outlet curves are within the measured intervals from the EnFlo tunnel and that can be considered to be one reason to why the model manages to describe the boundary layer in the tunnel. Moving further to the turbulent kinetic energy, the other reason to why the model is suitable for ABL modelling is evident, the kinetic energy maintains the shape and magnitude rather well for all heights. Given a little hard to decipher, there are three values for which the kinetic energy is outside of the measured EnFlo values.
However, Figure 8 shows something that can be regarded as a problem with this specific RANS modelling in regards to the atmospheric boundary layer. Since the eddy viscosity assume isotropic turbulence, i.e. that 𝑅𝑋𝑋= 𝑅𝑌𝑌 = 𝑅𝑍𝑍 we notice that given that this is the case for our simulation, it
is by far not the case for the wind tunnel. This problem is acknowledged throughout the simulation. Being confident that the double modifications to the standard k-ε model manage to model the atmospheric boundary layer in regards to the overall flow, the simple array is then introduced into the domain, of which there are a few traits that the model manages to capture. Referring to [18] who does an experimental study two cubes, we notice first in figure 9 the recirculation zones. There is one in front of the first cube, one on top of the first cube but not on the second cube; moreover there are vortices between the cubes and also a wake behind the second cube. All of these traits are expected from experiment as previously referred. This is also witnessed from the Turbulent Kinetic Energy in figure 10. With high turbulence in front of the first cube and on top of it, connecting on to the second cube and then going on into the wake. Ponder figures 11 and 12, which show that the air passes between the cubes fairly undisturbed according to figure 11, while figure 12 shows that the turbulent region in the centre shrinks in z.
Figure 13 and 14, which are taken slightly above the ground shows the previously mentioned recirculation zones in front of the two first cubes, between the cubes downstream and finally a wake. This is further seen in figures 15 and 16, where the aforementioned regions are much more distinct. Finally figures 17 and 18 show the recirculation zones on top of the first two cubes (but none on the two other). Noticeable in figure 18 is that there is a highly turbulent region behind the two final cubes, a region which occurs at a height higher than the cubes are. All in all the adapted turbulence model behaves much as one would expect to see in such a case.
39 in the boundary layer; on the other hand the measured values are around and near buildings. as such both values are rather important, and one could suspect the measurement positions should be carefully considered. Unfortunately this could not be explored extensively in this work.
This suspicion discussed in the previous paragraph is sort of enhanced in table 3. The measured values compared to the synthetic data is rather close along x when y=0.0. But for the two points on the sides of the cubes a big difference is observed. Granted this is probably due to the point source not modelling the real source well enough, and as such the gas has not yet dispersed enough before arriving at the array. This must be kept in mind when analysing the results for the source term estimation.
Arriving on the synthetic data test, take note of the shape of the probability distribution in figure 21, of which a point source contained in this region (with varying strength) is expected to yield the same values as the ones the forward model gave. Figure 22 itself shows where the maximum posteriori probability to find the point source, astonishingly enough the maximum probability is the precise cell which was set as the point source for the forward model. As such, it is concluded that the adjoint is valid method to determine the source. Figures 23 and 24 shows the least square predicted source strength, and what is mostly interesting here is the fact that it is possible to have a source downstream from the measurement, given a very high strength of the source, which by shear diffusion would allow the gas to be transported against the wind. Happy with the model evaluation, it is time to move on to the data from the wind tunnel.
Figures 25 and 26 shows the probability densities for the two best predicting turbulent Schmidt numbers. Notice two things, the first and most obvious being how the shape of the probability density changed between these two cases. The second notable trait is that, even though the real source was not pin pointed, it is possible to explain this prediction. It is known that this sort of source, i.e. a cylindrical source, could be modelled by a point source further upstream. Hence it could be suspected that these regions represent the region of which inside a point source would be the equivalence of the true cylindrical source. This is the most noteworthy conclusion from this simulation, as methods could be constructed that could find which sources would be equivalent to a point source inside this region.
Finishing off, what would happen if data was to be removed from the model? Since the two points 2 and 4 in figure 20 show the highest error compared to a point source, what happens when these are removed is rather interesting.
40 Table 4 is mainly interesting from the synthetic data viewpoint. Given the real source strength of 1.1 ∙ 10−3km−3, knowing that for the synthetic data the true source positon, the origin, and the HPD
coincide, notice how close to the real source strength the estimate is. Note the pin point source strength estimate for the height dependent Schmidt number at the origin, but where the HPD is much closer to the array and with a much lower source strength. However, it is quite interesting how the estimated source strength for the 5 measurements, both inside the HPD but also the origin, is twice the size of the real source. Perhaps a point source of this magnitude in this region corresponds to the real source, from the synthetic data strength estimation this is not too unlikely.
41
7 Concluding words
The Source Term Estimation, or reconstructing the source term, is indeed an ill-posed problem. Ideally, one would like to model the scalar product in (2.8.1) as a union of a domain; however, it is a combinatorial nightmare trying to go through all possible unions. Under the assumption that a point source could be used to model the true source in the wind tunnel, the assumption of finding a point source could very well be valid. Since a true source such as the one found in the wind tunnel could be modelled by a point source further upstream, one could conclude that the method is indeed a success. However, more research would be preferable. One method could be to coarsen the possible point source in (2.8.1). That is to look at the neighbours of each cell and take the average of the adjoint plume and also conclude whether it would be wise to perform this on the estimated source strength also, or one or the other.
42
8 Bibliography
1. Robertson, Lennart. Extended back-trajectories by means of adjoint equations. Norrköping : Swedish Meteorological and Hydrological Institute, SMHI, 2004.
2. Celik, Ismail B. Introductory Turbulence Modeling. Morgantown, WV : West Virginia University. Mechanical & Aerospace Engineering Dept., 1999.
3. Pope, Stephen B. Turbulent Flows. Edinburgh : Cambridge University Press, 2000.
4. Warsi, Z.U.A. Fluid Dynamics - Theoretical and Computational Approaches. Mississippi : CRC Press, 1992.
5. Yaglom, Monin and. Statistical Fluid Mechanics. The Massachusetts institute of Technology : the MIT press, 1971.
6. Marchuk, Guri I. Adjoint Equations and Analysis of Complex Systems. Institute of Numerical Mathematics, Russian academy of Sciences, Moscow, Russia : Kluwer, 1994.
7. N. Brännström, L.Å. Persson. A measure theoretic approach to the linear inverse atmospheric dispersion problems. Inverse Problems. 2015, 31.
8. Yee, E. Bayesian Inversion of Concentration Data for an Unknown Number of Contaminant Sources. s.l. : Defence Research and Development Canada, 2007.
9. Launder, B. E. and Spalding, D. B. The numerical computation of turbulent flows. Computer
Methods in Applied Mechanics and Engineering. 1974, Vol. 3(2), 269-289.
10. Masson, A. El Kasmi & C. Turbulence modeling of atmospheric boundary layer flow over complex terrain: a comparison of models at wind tunnel and full scale. Wind Energy. 2010, pp. 689-704. 11. M. Tsuchiya, , S.Murakami, A. Mochida, K.Kondo, Y.Ishada. Development of a new k-e model for flow and pressure fields around bluff body. Journal of Wind Engineering and Industrial Aerodynamics. 1997, Vols. 67-68, pp.169-182.
12. Wright, D.M. Hargreaves & B.G. On the use of the k-epsilon model in commercial CFD softare to model the neutral atmospheric boundary layer. Journal of Wind Engineering and Industrial
Aerodynamics. 2007, Vol. 95, 5. pp. 355-369.
13. Yoshihide Tominaga, Ted Strathopoulos. Turbulent Schmidt number for CFD analysis with various types of flowfield. Atmospheric Environment. 2007, pp. 8091-8099.
14. Thomas K. Flesch, John H. Prueger, Jerry L. Hatfield. Turbulent Schmidt number from a tracer experiment. Agricultural and Forest Meteorology. 111, 2002, pp. 299-307.
15. The height dependence of the turbulent Schmidt number within the boundary layer. Konrad Koeltzsch. Atmospheric Environment. 2000, Vol. 34, 7, pp. 1147-1151.
43 17. Alan Robins, Paul Hayden. The netural boundary layer in the EnFlo wind tunnel. s.l. : MODITIC, 2014.
18. Havel, Brian. Experimental Investigation of the turbulent flow generated by two interfering
surface mounted cubes in a thin boundary layer. Ontario : The University of Western Ontario, 1999.
19. G. Johannesson, B. Hanley, J. Nitao. Dynamic Bayesian Models via Monte Carlo - An Introduction
with Examples. s.l. : Lawrence Livermore National Saboratory, 2004.
i
A. Appendix
A.1 Dispersion comparisons of different Schmidt numbers
ii
iii
iv
v
vi
A.2 Convergence plots
Figure 38: Convergence plot for the empty domain of the modified ABL MMK k-ε model, converged to 𝟏𝟎−𝟔
vii
Figure 40: Convergence plot for the forward dispersion of neutral gas
