Adjoint-based Calibration of Inlet Boundary Condition for Atmospheric CFD Solvers

A continuous adjoint solver is developed for optimization of the inlet velocity profile boundary condition for CFD simulations of the neutral atmospheric boundary layer (ABL). The adjoint solver uses interior domain wind speed observations to compute the gradient of a calibration function with respect to inlet velocity components and wind direction. The solver has been implemented in the open source CFD package OpenFOAM. The sensitivities computed by the newly developed adjoint solver are validated against the second order finite-difference method. Furthermore, the DAKOTA optimization package is 5 coupled with OpenFOAM, and a number of numerical studies are carried out including the calibration of the inlet velocity profile of a 3D complex domain.


Introduction
The wind energy industry is growing very fast and a comprehensive site assessment is a key factor in planning, installation and performance of wind farms.As a result, the interaction of the atmospheric boundary layer (ABL) flow and wind turbines is one of the most important aspects of the site assessments for wind farms.Over the last decades with development of powerful computers the Computational Fluid Dynamics (CFD) has become one of the leading tools for micro-scale simulation of the wind flow over complex terrains.However, in general, developing algorithms that capture all the physics of such complex flow regimes is an ongoing research in CFD.
In CFD simulations of atmospheric flows the correct boundary conditions (BCs) are often unknown but measurement data (e.g wind speed, wind direction and etc.) within the area of interest is available.Using the measurement data and an inverse analysis, the unknown BCs can be obtained with an optimization algorithm.This approach is called open boundary optimization and has been successfully tried in oceanography (Seiler, 1993;Chen et al., 2013) and numerical wind prediction (NWP) models (Schneiderbauer and Pirker, 2011).Another approach is to use observations and statistical analysis (Glover et al., 2011) to calibrate the CFD model parameters (e.g., inflow and turbulence model constants) .
The solution to an optimization problem can be found with different solvers.However, the methods such as genetic algorithm and evolutionary strategies (Davis, 1991;Michalewicz, 1996) require a large number of function evaluations which in CFD applications can be computationally very expensive.Alternatively, the gradient-based optimizers (Ruder, 2016) use the derivative of cost function with respect to (w.r.t.) the design parameters.Then, using the gradients direction and values and a relatively small number of cost function evaluations the optimal design can be obtained.
To compute the gradients, the finite-difference (FD) method is relatively simple but prone to error due to step-size and round-off error.Moreover, for a high number of design parameters the method becomes prohibitively expensive.In the adjoint method, the sensitivity of the objective function can be calculated independently from the number of parameters, and this considerably reduces the cost of computation.
Since the first application of adjoint method into compressible CFD models by Pironneau in 1974, the adjoint based optimisation methods have been extensively used in shape and topology optimisation (Jameson et al., 1998;Kämmerer et al., 2003;Li et al., 2006;Othmer and Grahs, 2005;Othmer et al., 2006).The adjoint gradient computation can be categorized into two main groups: 1) continuous adjoint; 2) discrete adjoint.In the continuous method (Jameson, 1988) the adjoint equations are first derived and then discretized.In the discrete adjoint (Giles and Pierce, 2000), using the chain rule, the adjoint solver is derived by line-by-line differentiating of the discretized original (primal) flow solver code.The discrete adjoint differentiation can also be automated using Algorithmic Differentiation (AD) (Griewank and Walther, 2008).
In theory, both methods can be applied to any algorithm and model that is continuously differentiable.However, the manually derived discrete adjoint differentiation of big CFD codes is laborious and error-prone.Although, the AD tools can be seen as an interesting solution to this problem their application to the codes which are written in high-level languages (e.g.C++) is still limited in terms of memory requirement and performance.A comparison of these two adjoint approaches can be found here (Giles and Pierce, 2000;Nadarajah andJameson, 2000, 2001).
The discrete adjoint version of OpenFOAM based solvers has been presented before (Towara and Naumann, 2013;Towara et al., 2015).More recently a hybrid approach has also been introduced (He et al., 2018) in which some parts of the code are differentiated by finite-difference, and better performance is reported in comparison to the pure discrete adjoint version of the code.Despite all these improvements, the continuous adjoint version of OpenFOAM solvers are still more popular.
The effect of inflow boundary wind velocity and its direction are very important parameters in an ABL CFD simulation, but they are not often known.The focus of this paper is on the calibration of the inlet velocity profile and inflow wind direction (WD) for a neutral ABL flow and cases in which wind speed measurements near the area of interest at the site are present.The available continuous adjoint solver of OpenFOAM CFD tool package (adjointShapeOptimizationFoam), which is for topology optimization of duct flows, is further developed to compute the gradients.To the knowledge of the authors, the adjoint-based calibration of the inflow boundary has not been applied before in the framework of wind resource assessment.
The structure of the paper is as follows: The primal flow with forest model is explained in Section 2. In Section 3 the gradient-based optimization and the theory of adjoint method are briefly presented.Section 4 contains the adjoint equations and its BCs derivations.The inflow wind direction and its derivative are discussed in Section 5. Finally, the numerical results and the conclusions are presented in Sections 6 and 7.

Flow Model
The ABL flow model for cases of neutral stratification consists of steady-state Reynolds Averaged Navier-Stokes (RANS) equations for incompressible fluid flows (Rebollo and Lewandowski, 2014): (1) where R is the spatial residual of the equations that is driven to zero.The variables V and p are the state variables velocity and pressure, ν e f f stands for the sum of kinematic and turbulent viscosity and D is the rate of strain tensor, D = 1 2 (∇V + (∇V ) T ).

Inflow Boundary
Mainly two different methods are used to set the inflow boundary condition for ABL flow simulations: The first method is to derive a mathematical formula, either theoretically or empirically, which describes the vertical extrapolation of wind speeds V (z).At altitudes relevant to wind-power production (10-200 [m]), the theoretical formula is often based on the Monin-Obukhov similarity theory (MOST) (Foken, 2006), with height z, roughness length z 0 , Monin-Obukhov length L, friction velocity V * and von Karman constant κ = 0.40.The function Φ is a stability correction and for the neutral condition zero.On the other hand, the empirical power law gives where n is the friction coefficient or Hellmann exponent (Bañuelos-Ruedas et al., 2011).
The second method which is common in CFD RANS simulations is to solve a zero-pressure gradient equation for a onedimensional (1D) domain with periodic boundary conditions in stream-and span-wise directions (Chang et al., 2018): The inflow profile and its turbulent characteristics are obtained from this 1D simulation and then are mapped to the 3D inlet boundary.

Forest Effect
Forest canopies modify the available free volume of the terrain domain and introduce an explicit drag term to the momentum equation as below: with density ρ, leaf-level canopy drag coefficient C d and leaf area per unit volume A(z).The effects of the forest in the turbulence models such as k − ε for ABL flows has been extensively discussed in the literature and several formulas are presented to make the turbulence model consistent with the modified momentum equation.For more details, the reader may be referred to (Lopes da Costa, 2007).

Gradient-Based Calibration
Calibration algorithms seek to maximize the agreement between simulation outputs and measurements.In the context of ABL based model calibration the data are often wind speed and direction at one or more locations of a potential wind farm site.The CFD-based calibration can be formulated as a constrained optimization problem with a scalar objective function: where R stands for the spatial residual of the flow equations with V and p the discretized velocity and pressure, respectively.
V M i and V S i are the measured and simulated wind velocities at the same location in the domain.The variable α represents the design variables which are considered to be the velocity at inlet faces of the CFD mesh through this work.
In gradient-based optimisation, one needs to compute the gradient (dJ/dα) of this cost function with respect to the design parameters at each design step, where A is an optimization algorithm operator that returns a perturbation δ α to the current design α n .The procedure is repeated until a convergence criterion is reached.

Gradient Computation
Computing the term dJ dα includes the differentiation of the steady-state PDE of flow equations, R = 0, where ψ stands for state variables V and p.The Eq. ( 9) results in a linear system, after which the total derivative can be computed by The calculation of gradient by Eq. ( 11) is called the tangent-linear approach.Comparing to the finite difference (FD) method, which is error-prone, the derivative computed by the tangent-linear method is exact.However, every change in the design variables, α, requires the solution to the linear system of Eq. ( 10).This makes the tangent-linear approach, similar to the FD, computationally expensive when the number of design parameters is not small.

Adjoint Method
The main idea of the adjoint method is to multiply the PDE of flow equations with new variables such that the variations of state variables, ∂ ψ ∂ α , vanish.By introducing the adjoint variables U and q, the cost function can be reformulated to a Lagrange function as where Ω is the flow domain.Then the sensitivity of the cost function can be given by which excludes the state variables sensitivities.Consequently, the computation of sensitivity becomes independent of the number of design parameters and, compared to other methods, relatively cheap .However, a new set of equations namely adjoint equations needs to be solved.

Adjoint Solver for ABL Inflow Calibration
The theory presented here is based on the work of C. Othmer 2008 which derives an adjoint solver for topology optimization of duct flows to reduce the pressure loss between inlet and outlet boundaries.However, there are some fundamental differences between the current work and that adjoint solver that can be summarized as followings: a) in the duct flow optimization the cost function is pressure loss, while here the objective function is defined for calibration as expressed in Eq. ( 7

Adjoint Equations
To zero out the variation of state variables w.r.t.design parameters and by neglecting the turbulent viscosity variation, assuming the "frozen-turbulent" hypothesis, the following expression should hold: or equivalently Decomposition of parts into interior domain, Ω, and its boundaries, Γ, leads to reformulation of Eq. ( 15) as follows Due to the definition of the cost function Eq. ( 7), its direct variation comes only from the interior domain.Moreover, it does not have any derivative w.r.t. the pressure field.The corresponding terms are zeroed out in Eq. ( 16).The only cost function derivative that remains is w.r.t. the inflow and velocity in the interior of the domain and at the locations where the measurements are available.From the latter we have In addition, as it was explained before, the forest model adds an explicit source term, Eq. ( 6), to the momentum equation.After division by density, the derivative of the forest source term gives Using Eqs.(16, 17 and 18) the adjoint equations can be derived as where ω j is the cell volume.In comparison to derived equations in (Othmer, 2008), the design variables, α, are not anymore in the adjoint momentum equation.Instead, the difference between measurement and simulation at observation locations and the differentiated term of forest model are present as source.

Boundary Conditions
The boundary integrals of Eq. ( 16) can be mathematically re-formulated and reduced to where n is the unit normal vector from the boundary faces.The adjoint BCs should be chosen such that the above equations are held meaning that after any small perturbations of flow fields, V + δV and p + δ p, the primal equations are still satisfied.
Except for the inlet, which is the design parameter.
Generally, for an ABL CFD domain no-slip wall (zero fixed velocity) and zero pressure gradient conditions are imposed on the ground.For a wall type of boundary in which δV is zero the first integral of Eqn. ( 22) is cancelled.Then, the only way to satisfy the following conditions is to apply a no slip (U = 0) condition on the ground.No BC can be derived on the ground for the adjoint pressure but consistent with the primal a zero gradient condition is applied.
For the top and outlet boundaries of the domain a zero gradient velocity ((n • ∇)δV = 0) and zero fixed pressure (p = 0) are the common conditions for the primal system.These conditions fulfil the Eq. ( 21) and cancel the second integral of Eq. ( 22).
The only term that remains is the first term of the Eq. ( 22) which needs to be zeroed out, After decomposition into tangent and normal components it can be shown that the relations below should hold where underscripts n and t represent the normal and in-plane components respectively.The adjoint BCs can be summarized as ground (wall): It is worth mentioning that the last term of the adjoint pressure which includes the kinematic viscosity, in implementation is often neglected (Nilsson et al., 2014).Moreover, the adjoint variables at the inlet should not be chosen to zero out the inlet velocity perturbations because the design variables are the inlet velocities, δV inlet = δ α.Instead, the zero gradient condition is imposed on the inlet for both adjoint velocity and adjoint pressure to have a well-posed system.Finally, from the integral over the boundary term in Eq. ( 16) it is clear one needs to evaluate the following expression to compute the sensitivity.
As it was mentioned before, in ABL CFD simulations it is common to simulate first a 1D domain with periodic boundary to obtain the inflow boundary condition.The mesh of the 1D precursor run follows the vertical structure of the 3D domain as it is shown in Figure 1.This implies that the number of faces in each column of the 3D inflow boundary is the same as the number of cells in the 1D mesh.Moreover, and ideally, the face center heights in the 3D mesh are equal to their counterpart cell height in 1D.The inflow wind direction (WD) effect can be expressed by a rotation angle, θ , which rotates the inflow from its default west to east (WD=270 • ) direction, The differentiation of Eqs.(32 and 33) gives The adjoint solver which was explained in the previous section computes the derivative of the cost function w.r.t.3D inflow velocities at each face of the boundary: Using the chain rule, one can compute the sensitivity w.r.t each cell of the 1D inflow velocity as follows The Eq. ( 38) means the x and y gradient components of the 3D inflow faces at the same column j are accumulated and then multiplied by cos(θ ) and sin(θ ) respectively before being summed.For the sake of clarity, the Eq. ( 38) can be re-written as Using the same analogy, it can be shown that the sensitivity w.r.t rotation angle can be obtained by where the dot sign stands for the inner product of the two vectors.

Numerical Results
The adjoint solver and Eq. ( 30) are implemented based on the "simpleFoam" incompressible CFD solver of OpenFOAM-4.1.
In this section, first the accuracy of the gradients obtained by the developed solver is tested against the second-order FD method in a simple 3D domain.Then, the inflow calibration of a real complex terrain is presented.For all the simulations in this work the general roughness length of the domain is z 0 = 0.05 [m] and the turbulent eddies are modelled by standard k-ε model with canopy model of Liu et. al 1996.Moreover, an ABL wall function is used to apply the roughness-related logarithmic law near to the ground, which is consistent with Monin-Obukhov similarity theory (Chang et al., 2018).For the primal CFD simulation of a circular domain the standard "inletOutlet" BC is used which checks if the flow is flowing into domain or out of it and switches between fixed value and zero gradient respectively.This BC is further developed to apply the derived adjoint BCs in a similar way and based on the flow direction on the boundary.Table 1.Comparison of the wind direction gradients by finite differences (FD) and via the adjoint approach.

Gradient Verification
For gradient evaluation, the 1D velocity profile inflow is rotated by 30 • (WD=240 • ).Then, using the simulation result of this set-up and the reference data of the first simulation the sensitivities are computed.The gradients obtained by the developed adjoint solver are plotted against the second order FD gradients in Figure 4.In general, the trends of the sensitivity profiles are similar.Moreover, the gradients are in excellent agreement except only for the heights between 50 [m] and 100 [m] in which the maximum relative error is ε rel = 0.10.This difference can be traced back to two main reasons: a) even though the second order FD can be regarded as a relatively good reference it is still not exact; b) in the derivation of the adjoint solver equations and BCs there are some simplifications such as the assumption of the frozen turbulent.
The wind direction gradients are tabulated in Table 1.The derivative of calibration cost function w.r.t. the change in wind direction is much higher that w.r.t. the inflow boundary.Here also, there is a good agreement between the FD and the adjoint gradients and the relative error of wind direction sensitivity is close to that of inflow gradients.This is of course not surprising solvers provide the required information whenever it is needed and this process continues until a certain convergence criterion is satisfied.
As it was explained in section 2.1, the inflow boundary of an ABL domain can be in general represented by a logarithmic or power-law relationship.During optimization, the optimizer may ask for a cost function evaluation with a new inflow boundary which is highly unrealistic for an ABL domain.This may result in a poor numerical convergence or even divergence of the solver.To avoid such scenarios, the curve fit capability from the Scipy library of Python is used to smooth and fit the new inflow to a boundary which has a logarithmic/power law characteristic.This part is called smoothing in the optimization flow chart.If the new inflow boundary is severely skewed such that it could not be fitted to either of those types of functions then the optimizer restart the optimization from the last fitted inflow boundary and asks for the gradient evaluation to continue.
The turbulent properties of the inflow boundary are not part of the design parameters, but instead, in each new flow solver call of the optimization the turbulent parameters of the k-ε model inside of the domain are initialized with the last converged solution.In this way the turbulence model parameters are also gradually updated toward the end the of optimization when the inlet boundary velocities have reached their optimum value.
The convergence criterion is defined such that the optimization stops when the absolute difference between simulated and measured value at all heights is below 0.2 [ms −1 ]. Figure 7 shows the history of optimization.The optimization has called for a total number of 41 primal solver runs for cost function evaluation.In addition 8 adjoint gradient evaluations were required.
The optimization convergence graph shows that there is a spike both in the cost function and the wind direction.This can be explained by the fact that in early iterations the derivative of cost function w.r.t the WD is much bigger than w.r.The initial and optimal results are compared in Figure 8.The optimized inlet velocity is almost the same as the reference profile.Indeed, the output of the optimizer could be the exact reference profile if the convergence criterion was stricter.The velocity profile near to MM200 mast, is in good agreement with the pseudo measurements at all five selected heights.

Conclusions
In this paper, it has been shown that the ABL CFD solvers can be calibrated via adjoint-based inlet boundary optimization.
Based on the frozen-turbulence hypothesis, the adjoint equations and its boundary conditions for such problems are derived.
The developed solver has been coupled with the DAKOTA optimization package and applied to the 3D Kassel terrain for a 14 Wind Energ.Sci.Discuss., https://doi.org/10.5194/wes-2019-6Manuscript under review for journal Wind Energ.Sci. Discussion started: 11 February 2019 c Author(s) 2019.CC BY 4.0 License.
); b) instead of the volume porosity distribution here the design parameters are the velocity components at inflow boundary faces; c) in this study the adjoint equations and boundary conditions are derived based on the definition of the calibration function and the characteristics of the ABL domain including the forest model.

Figure 1 .
Figure 1.The results from 1D precursor run are copied to the 3D inflow boundary of the main domain which has the similar number of cells in vertical direction.

37) 8
Wind Energ.Sci.Discuss., https://doi.org/10.5194/wes-2019-6Manuscript under review for journal Wind Energ.Sci. Discussion started: 11 February 2019 c Author(s) 2019.CC BY 4.0 License.where i and j represent the row-and column-wise position of a face on the boundary and the total number of the faces in the 3D circular boundary is N = n × m.

A
cylindrical domain with 1000 [m] radius and 300 [m] height is chosen.The mesh of domain has 209K hexahedral cells in which the inflow/outflow boundary has 49 rows and 172 columns and a total number of 8428 faces (N = 49 × 172 = 8428).The topoSet utility of OpenFOAM is used to select a number of cells as forest in a box size of 300 [m] in x and y and 40 [m] in z direction.The canopy drag and leaf area density for all forest cells are C d = 0.3 and A = 0.0033 [m −1 ].The operational Reynolds number based on the free-stream velocity, forest area height and air kinematic viscosity is Re h = 4.8 × 10 5 .

Figure 2 .
Figure 2. The reference simulation (WD=270 • ) of the circular domain with forest: 1D inflow boundary (left) and the CFD simulated velocity at z = 45 [m] above the ground (right).

Figure 3 .
Figure 3.The reference simulation (WD=270 • ) of the circular domain with forest: velocity profile and selected reference speeds (left) and the cross sectional (x-z) view of velocity field at the center of domain, y = 0, (right).The green grid cells represent the forest region.

Figure 4 .
Figure 4. Comparison of the inflow boundary gradients by finite differences (FD) and via the adjoint approach.

Figure 6 .
Figure 6.The calibration flowchart of the inflow boundary using DAKOTA optimization package and the developed adjoint solver.
Wind Energ.Sci.Discuss., https://doi.org/10.5194/wes-2019-6Manuscript under review for journal Wind Energ.Sci. Discussion started: 11 February 2019 c Author(s) 2019.CC BY 4.0 License.neutral stratification condition.Using the reference wind speeds and its direction over a hill within the terrain, the optimal inlet velocities have been found which are in good agreement with the inflow boundary profile of the reference simulation.The presented adjoint solver has the potential to be further developed by including Coriolis force, turbulence model and thermal stratification.