Adjoint-based calibration of inlet boundary condition for atmospheric computational fluid dynamics solvers

A continuous adjoint solver is developed for calibration of the inlet velocity profile boundary condition (BC) for computational fluid dynamics (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 speed and wind direction. The solver has been implemented in the open-source CFD package OpenFOAM coupled with the local gradient-based “CONMIN-frcg” solver of the DAKOTA optimization package. The feasibility of the optimizer output is continuously monitored during the calibration process. The inlet flow profile is considered acceptable only if it can be fitted to a logarithmic or power law function with a tolerance of 3 %. Otherwise, the optimization takes the last fitted profile and asks for a new gradient evaluation. The newly developed framework has been applied in two cases, namely the Ishihara case and Kassel domain. By using the measurements over the hill in the Ishihara case, the method was able to predict the velocity profiles upstream and downstream of the hill accurately. For the Kassel domain, despite the complexity of the site, the method managed to achieve the targeted profile within a reasonable number of the solver calls.

1 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 the development of powerful computers the computational fluid dynamics (CFD) has become one of the leading tools for microscale 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 (wind speed, wind direction, 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 methods. 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 the design parameters. Then, the optimal solution can be obtained using the gradient and a relatively less cost function evaluation. Gradient computation with the finite-difference (FD) method is relatively simple. However, for a high number of design parameters it is 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 the adjoint method into compressible CFD models by Pironneau (1974), the adjoint-based optimization methods have been extensively used in shape and topology optimization (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 and (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 differentiation 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 principle, both methods can be applied to any algorithm and model which 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 and Jameson, 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 adjoint method has been well applied to different problems in atmospheric science such as wind turbine blade shape optimization (Dhert et al., 2017), wind-farm control (Goit et al., 2016; Munters and Meyers, 2017; Vali et al., 2019) and wind-farm layout optimization (King et al., 2017). Bauweraerts and Meyers (2018) have used the adjoint method and large eddy simulation (LES)-based 4D-Var data assimilation to estimate the turbulent flow field of an atmospheric boundary layer domain from lidar data. Although compared to the Reynolds Averaged Navier–Stokes (RANS) models, the LES can provide more accurate predictions of the flow field, it is still mainly used as a research tool due to its demanding computing power.

The effect of inflow boundary wind speed and its direction are significant parameters in ABL CFD simulation, but they are not often known. The focus of this paper is on the adjoint gradient-based calibration of the inlet velocity profile and inflow wind direction (WD) for a RANS-based ABL flow, which is a very common CFD solver in the wind energy industry, with only a few wind speed measurements from a metrological mast at the site. 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.

The structure of the paper is as follows: the gradient-based optimization and the theory of the adjoint method are briefly presented in Sect. 2. The derivation of the adjoint equations and its BCs from the primal flow model is explained in Sect. 3. Finally, the numerical results and the conclusions are presented in Sects. 4 and 5.

2 Gradient-based optimization and adjoint method theory

In gradient-based optimization or calibration, one needs to compute the gradient of a smooth cost function, J, with respect to the design parameters, α, at each design step,

$$\begin{array}{}\text{(1)}& {\mathit{\alpha }}^{n+\mathrm{1}}={\mathit{\alpha }}^n+\mathbb{A}\left({\mathit{\alpha }}^n,\left({\displaystyle \frac{\mathrm{d}J}{\mathrm{d}\mathit{\alpha }}}\right)\right)={\mathit{\alpha }}^n+(\mathrm{\Delta }\mathit{\alpha }{)}^n,\end{array}$$

where 𝔸 is an optimization algorithm operator that returns a perturbation Δα to the current design αⁿ. The procedure is repeated until a convergence criterion is reached. The design parameters, α, are chosen based on the optimization problem and its parametrization.

In CFD applications, computing the term $\frac{\mathrm{d}J}{\mathrm{d}\mathit{\alpha }}$ includes the differentiation of the steady-state partial differential equation (PDE) governing equation of the flow,

$$\begin{array}{}\text{(2)}& \mathit{R}(\mathit{\psi },\mathit{\alpha })=\mathrm{0}\to {\displaystyle \frac{\partial \mathit{R}}{\partial \mathit{\psi }}}{\displaystyle \frac{\partial \mathit{\psi }}{\partial \mathit{\alpha }}}+{\displaystyle \frac{\partial \mathit{R}}{\partial \mathit{\alpha }}}=\mathrm{0},\end{array}$$

where R is the residual vector of the discretized flow equations that is driven to zero and ψ stands for state variables (velocity, pressure, temperature, etc.). Eq. (2) results in a linear system,

$$\begin{array}{}\text{(3)}& {\displaystyle \frac{\partial \mathit{R}}{\partial \mathit{\psi }}}{\displaystyle \frac{\partial \mathit{\psi }}{\partial \mathit{\alpha }}}=-{\displaystyle \frac{\partial \mathit{R}}{\partial \mathit{\alpha }}},\end{array}$$

in which the term $\frac{\partial \mathit{R}}{\partial \mathit{\psi }}$ is Jacobian and $\frac{\partial \mathit{\psi }}{\partial \mathit{\alpha }}$ represents the perturbation of flow fields. Using the chain rule, the total derivative can be then computed by

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{\mathrm{d}J}{\mathrm{d}\mathit{\alpha }}}={\displaystyle \frac{\partial J}{\partial \mathit{\alpha }}}+{\displaystyle \frac{\partial J}{\partial \mathit{\psi }}}{\displaystyle \frac{\partial \mathit{\psi }}{\partial \mathit{\alpha }}}.\end{array}$$

Several methods can be used to compute the gradient from Eq. (4). For instance the complex variable technique, which overcomes the problem of choosing step width in the finite-difference method, or the forward mode (tangent linearization) application of algorithmic differentiation (AD). However, all these methods are computationally expensive when the design space is large. This is due the fact that in Eq. (4) the term $\frac{\partial \mathit{\psi }}{\partial \mathit{\alpha }}$ requires an expensive PDE simulation, which satisfies Eq. (3) for every dimension in the design space, α_i. As will be shown in the following, in the adjoint method the sensitivity can be obtained without computing this expensive term.

From Eq. (3) we can write the perturbation term as

$$\begin{array}{}\text{(5)}& {\displaystyle \frac{\partial \mathit{\psi }}{\partial \mathit{\alpha }}}=-{\left({\displaystyle \frac{\partial \mathit{R}}{\partial \mathit{\psi }}}\right)}^{-\mathrm{1}}{\displaystyle \frac{\partial \mathit{R}}{\partial \mathit{\alpha }}}\end{array}$$

leading to

$$\begin{array}{}\text{(6)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{d}J}{\mathrm{d}\mathit{\alpha }}}& ={\displaystyle \frac{\partial J}{\partial \mathit{\alpha }}}-{\displaystyle \frac{\partial J}{\partial \mathit{\psi }}}{\left({\displaystyle \frac{\partial \mathit{R}}{\partial \mathit{\psi }}}\right)}^{-\mathrm{1}}{\displaystyle \frac{\partial \mathit{R}}{\partial \mathit{\alpha }}}={\displaystyle \frac{\partial J}{\partial \mathit{\alpha }}}\\ & +\left[-{\displaystyle \frac{\partial J}{\partial \mathit{\psi }}}{\left({\displaystyle \frac{\partial \mathit{R}}{\partial \mathit{\psi }}}\right)}^{-\mathrm{1}}\right]{\displaystyle \frac{\partial \mathit{R}}{\partial \mathit{\alpha }}}.\end{array}\end{array}$$

The terms in the bracket can be identified as the adjoint system of equations from which the adjoint variable, $\stackrel{\mathrm{‾}}{\mathit{\psi }}$, can be introduced as

$$\begin{array}{}\text{(7)}& \begin{array}{rl}& {\stackrel{\mathrm{‾}}{\mathit{\psi }}}^T=-{\displaystyle \frac{\partial J}{\partial \mathit{\psi }}}{\left({\displaystyle \frac{\partial \mathit{R}}{\partial \mathit{\psi }}}\right)}^{-\mathrm{1}}\mathrm{or}\\ & {\left({\displaystyle \frac{\partial \mathit{R}}{\partial \mathit{\psi }}}\right)}^T\stackrel{\mathrm{‾}}{\mathit{\psi }}=-{\left({\displaystyle \frac{\partial J}{\partial \mathit{\psi }}}\right)}^T.\end{array}\end{array}$$

In this way, instead of solving a PDE simulation for every design variable, the adjoint system of equations needs to be solved only once. As a result, the computational cost of the gradient becomes independent of the number of design parameters (Giles et al., 2003; Mavriplis, 2007; Nielsen et al., 2010):

$$\begin{array}{}\text{(8)}& {\displaystyle \frac{\mathrm{d}J}{\mathrm{d}\mathit{\alpha }}}={\displaystyle \frac{\partial J}{\partial \mathit{\alpha }}}+{\stackrel{\mathrm{‾}}{\mathit{\psi }}}^T{\displaystyle \frac{\partial \mathit{R}}{\partial \mathit{\alpha }}}.\end{array}$$

3 Derivation of the continuous adjoint solver for ABL inflow calibration

3.1 Flow model

The ABL flow model for cases of neutral stratification consists of steady-state Reynolds Averaged Navier–Stokes (RANS) for incompressible fluid flows (Rebollo and Lewandowski, 2014), which results in the following equations for momentum and continuity:

$$\begin{array}{}\text{(9)}& {\displaystyle }& {\displaystyle }{\left(R_{\mathrm{1}},R_{\mathrm{2}},R_{\mathrm{3}}\right)}^T=(\mathit{U}\cdot \mathrm{\nabla })\mathit{U}+\mathrm{\nabla }p-\mathrm{\nabla }\cdot \left(\mathrm{2}{\mathit{\nu }}_{\mathrm{eff}}\mathbf{D}\left(\mathit{U}\right)\right),\text{(10)}& {\displaystyle }& {\displaystyle }{R}_{\mathrm{4}}=-\mathrm{\nabla }\cdot \mathit{U},\end{array}$$

where (R₁, R₂, R₃) and R₄ denote the discretized flow equations, R=(R₁, R₂, R₃, R₄)^T, in Eq. (2). The variables U and p are the state variables velocity vector and modified pressure, ν_eff stands for the sum of kinematic and turbulent viscosity, and D is the rate of strain tensor, $\mathbf{D}=\frac{\mathrm{1}}{\mathrm{2}}(\mathrm{\nabla }\mathit{U}+(\mathrm{\nabla }\mathit{U}{)}^T)$. Throughout this work, the standard k−ϵ turbulence model with the addition of a forest model is used.

3.1.1 Inflow boundary

The properties of the inflow boundary have an important effect on the solution of the interior domain. With the increasing application of LES, hybrid RANS-LES and direct numerical simulation (DNS) methods, the inflow turbulence generation has become the subject of many research works in recent decades. For a review of such methods and their application in wind energy, we refer to Wu (2017) and Stevens and Meneveau (2017). However, due to computational cost, it is still popular in the wind energy industry to use the RANS model with an inflow boundary obtained from either an analytical formula or a one-dimensional (1-D) simulation. The idea of the latter method is to solve a zero-pressure gradient equation for a 1-D domain with periodic boundary conditions in the stream- and span-wise directions (Chang et al., 2018):

$$\begin{array}{}\text{(11)}& {\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{{\mathit{U}}_x^{\prime }{\mathit{U}}_z^{\prime }}}{\partial z}}=\mathrm{0},\end{array}$$

where x (horizontal) and z (vertical) are the Cartesian coordinates. The inflow profile and its turbulent characteristics are obtained from this 1-D simulation and then are mapped to the 3-D inlet boundary.

3.1.2 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:

$$\begin{array}{}\text{(12)}& {\mathit{F}}_{\mathrm{D}}=-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{\rho }{C}_{\mathrm{d}}A\left(z\right)\left|\mathit{U}\right|\mathit{U}\end{array}$$

with density ρ, leaf-level canopy drag coefficient C_d and leaf area density A(z). The effects of the forest in the turbulence models such as k−ϵ for ABL flows have 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).

3.2 Adjoint model

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:

$$\begin{array}{}\text{(13)}& \begin{array}{rl}\mathrm{minimize}& J\left({\mathit{U}}_{\mathrm{M}},{\mathit{U}}_{\mathrm{S}},\mathit{\alpha }\right)=\sum _{i=\mathrm{1}}^k{\left[{\mathit{U}}_{{\mathrm{M}}_i}-{\mathit{U}}_{{\mathrm{S}}_i}\right]}^{\mathrm{2}};\\ & \mathrm{subject}\mathrm{to}\mathit{R}(\mathit{U},p,\mathit{\alpha })=\mathrm{0},\end{array}\end{array}$$

where R stands for the spatial residual of the flow equations with U and p the discretized velocity and pressure, respectively. ${\mathit{U}}_{{\mathrm{M}}_i}$ and ${\mathit{U}}_{{\mathrm{S}}_i}$ are the measured and simulated wind velocities at the same location in the domain, respectively. The variable α represents the design variables which are considered to be the velocity at inlet faces of the CFD mesh through this work.

The derivation of adjoint equations as in the preceding section is arguably the most straightforward way to introduce the adjoint equations and understand their advantages. However, the first developments for using the adjoint equations in CFD applications were done using a Lagrange multiplier argument (Hinze et al., 2008). From this point of view, the inner product of the PDE of flow equations and a new set of variables vanishes the variations of state variables, $\frac{\partial \mathit{U}}{\partial \mathit{\alpha }}$ and $\frac{\partial p}{\partial \mathit{\alpha }}$. By introducing the adjoint variables $\stackrel{\mathrm{‾}}{\mathit{U}}$ and $\stackrel{\mathrm{‾}}{p}$ for adjoint velocity and adjoint pressure, respectively, ${\stackrel{\mathrm{‾}}{\mathit{\psi }}}^T=(\stackrel{\mathrm{‾}}{\mathit{U}}$, $\stackrel{\mathrm{‾}}{p})$, the cost function can be reformulated to a Lagrange function as

$$\begin{array}{}\text{(14)}& \begin{array}{rl}L& :=J+\underset{\mathrm{\Omega }}{\int }\left({\stackrel{\mathrm{‾}}{U}}_x{R}_{\mathrm{1}}^T+{\stackrel{\mathrm{‾}}{U}}_y{R}_{\mathrm{2}}^T+{\stackrel{\mathrm{‾}}{U}}_z{R}_{\mathrm{3}}^T+\stackrel{\mathrm{‾}}{p}{R}_{\mathrm{4}}^T\right)d\mathrm{\Omega }\\ & =J+\underset{\mathrm{\Omega }}{\int }(\stackrel{\mathrm{‾}}{\mathit{U}},\stackrel{\mathrm{‾}}{p})\mathit{R}d\mathrm{\Omega },\end{array}\end{array}$$

where ${\stackrel{\mathrm{‾}}{U}}_x$, ${\stackrel{\mathrm{‾}}{U}}_y$ and ${\stackrel{\mathrm{‾}}{U}}_z$ are the adjoint velocity components and Ω is the flow domain. For the sensitivities of the cost function with respect to the design parameters, we have to compute the total variation of L:

$$\begin{array}{}\text{(15)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{d}L}{\mathrm{d}\mathit{\alpha }}}=& {\displaystyle \frac{\partial J}{\partial \mathit{\alpha }}}+{\displaystyle \frac{\partial J}{\partial \mathit{U}}}{\displaystyle \frac{\partial \mathit{U}}{\partial \mathit{\alpha }}}+{\displaystyle \frac{\partial J}{\partial p}}{\displaystyle \frac{\partial p}{\partial \mathit{\alpha }}}+\underset{\mathrm{\Omega }}{\int }\\ & \left[(\stackrel{\mathrm{‾}}{\mathit{U}},\stackrel{\mathrm{‾}}{p})\left({\displaystyle \frac{\partial \mathit{R}}{\partial \mathit{U}}}{\displaystyle \frac{\partial \mathit{U}}{\partial \mathit{\alpha }}}+{\displaystyle \frac{\partial \mathit{R}}{\partial p}}{\displaystyle \frac{\partial p}{\partial \mathit{\alpha }}}\right)+(\stackrel{\mathrm{‾}}{\mathit{U}},\stackrel{\mathrm{‾}}{p}){\displaystyle \frac{\partial \mathit{R}}{\partial \mathit{\alpha }}}\right]d\mathrm{\Omega }.\end{array}\end{array}$$

Choosing the Lagrange multipliers such that the variation with respect to the state variables vanishes, leads to

$$\begin{array}{}\text{(16)}& {\displaystyle \frac{\partial J}{\partial \mathit{U}}}{\displaystyle \frac{\partial \mathit{U}}{\partial \mathit{\alpha }}}+{\displaystyle \frac{\partial J}{\partial p}}{\displaystyle \frac{\partial p}{\partial \mathit{\alpha }}}+\underset{\mathrm{\Omega }}{\int }\left[(\stackrel{\mathrm{‾}}{\mathit{U}},\stackrel{\mathrm{‾}}{p})\left({\displaystyle \frac{\partial \mathit{R}}{\partial \mathit{U}}}{\displaystyle \frac{\partial \mathit{U}}{\partial \mathit{\alpha }}}+{\displaystyle \frac{\partial \mathit{R}}{\partial p}}{\displaystyle \frac{\partial p}{\partial \mathit{\alpha }}}\right)\right]d\mathrm{\Omega }=\mathrm{0}.\end{array}$$

Then the sensitivity of the cost function can be given by

$$\begin{array}{}\text{(17)}& {\displaystyle \frac{\mathrm{d}L}{\mathrm{d}\mathit{\alpha }}}={\displaystyle \frac{\partial J}{\partial \mathit{\alpha }}}+\underset{\mathrm{\Omega }}{\int }(\stackrel{\mathrm{‾}}{\mathit{U}},\stackrel{\mathrm{‾}}{p}){\displaystyle \frac{\partial \mathit{R}}{\partial \mathit{\alpha }}}d\mathrm{\Omega },\end{array}$$

which excludes the state variable sensitivities.

The theory presented here is based on the work of Othmer (2008), which derives an adjoint solver for topology optimization of duct flows to reduce the pressure loss between inlet and outlet boundaries. A more detailed derivation of the model is given by Hinterberger and Olesen (2011). By neglecting the turbulent viscosity variation, assuming the “frozen turbulence” hypothesis; replacing the derivative of Eq. (9) with the forest source term; and inserting Eq. (10) into Eq. (16), this gives

$$\begin{array}{}\text{(18)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\partial J}{\partial \mathit{U}}}& {\displaystyle \frac{\partial \mathit{U}}{\partial \mathit{\alpha }}}+{\displaystyle \frac{\partial J}{\partial p}}{\displaystyle \frac{\partial p}{\partial \mathit{\alpha }}}+\underset{\mathrm{\Omega }}{\int }\stackrel{\mathrm{‾}}{\mathit{U}}\cdot \left[\left({\displaystyle \frac{\partial \mathit{U}}{\partial \mathit{\alpha }}}\cdot \mathrm{\nabla }\right)\mathit{U}+(\mathit{U}\cdot \mathrm{\nabla }){\displaystyle \frac{\partial \mathit{U}}{\partial \mathit{\alpha }}}-\mathrm{\nabla }\right.\\ & \left.\cdot \left(\mathrm{2}{\mathit{\nu }}_{\mathrm{eff}}\mathbf{D}\left({\displaystyle \frac{\partial \mathit{U}}{\partial \mathit{\alpha }}}\right)\right)+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{C}_{D}A\left(\left|{\displaystyle \frac{\partial \mathit{U}}{\partial \mathit{\alpha }}}\right|\mathit{U}+\left|\mathit{U}\right|{\displaystyle \frac{\partial \mathit{U}}{\partial \mathit{\alpha }}}\right)\right]d\mathrm{\Omega }\\ & -\underset{\mathrm{\Omega }}{\int }\stackrel{\mathrm{‾}}{p}\mathrm{\nabla }\cdot {\displaystyle \frac{\partial \mathit{U}}{\partial \mathit{\alpha }}}d\mathrm{\Omega }+\underset{\mathrm{\Omega }}{\int }\stackrel{\mathrm{‾}}{\mathit{U}}\cdot \mathrm{\nabla }{\displaystyle \frac{\partial p}{\partial \mathit{\alpha }}}d\mathrm{\Omega }=\mathrm{0}.\end{array}\end{array}$$

Decomposition of parts into interior domain, Ω, and its boundaries, Γ, leads to reformulation of Eq. (18) as follows

$$\begin{array}{}\text{(19)}& {\displaystyle }\begin{array}{rl}\underset{\mathrm{\Gamma }}{\int }& \left[\stackrel{\mathrm{‾}}{\mathit{U}}\cdot \mathit{n}+{\displaystyle \frac{\partial J_{\mathrm{\Gamma }}}{\partial p}}\right]{\displaystyle \frac{\partial p}{\partial \mathit{\alpha }}}d\mathrm{\Gamma }+\underset{\mathrm{\Gamma }}{\int }\left[\mathit{n}\left(\stackrel{\mathrm{‾}}{\mathit{U}}\cdot \mathit{U}\right)+\stackrel{\mathrm{‾}}{\mathit{U}}(\mathit{U}\cdot \mathit{n})\right.\\ & \left.+\mathrm{2}{\mathit{\nu }}_{\mathrm{eff}}\mathit{n}\cdot \mathbf{D}\left(\stackrel{\mathrm{‾}}{\mathit{U}}\right)-\stackrel{\mathrm{‾}}{p}\mathit{n}+{\displaystyle \frac{\partial J_{\mathrm{\Gamma }}}{\partial \mathit{U}}}\right]{\displaystyle \frac{\partial \mathit{U}}{\partial \mathit{\alpha }}}d\mathrm{\Gamma }\\ & +\underset{\mathrm{\Gamma }}{\int }\left[-\mathrm{2}{\mathit{\nu }}_{\mathrm{eff}}\mathit{n}\cdot \mathbf{D}\left({\displaystyle \frac{\partial \mathit{U}}{\partial \mathit{\alpha }}}\right)\cdot \stackrel{\mathrm{‾}}{\mathit{U}}\right]d\mathrm{\Gamma }\\ & +\underset{\mathrm{\Omega }}{\int }\left[-\mathrm{\nabla }\cdot \stackrel{\mathrm{‾}}{\mathit{U}}+{\displaystyle \frac{\partial J_{\mathrm{\Omega }}}{\partial p}}\right]{\displaystyle \frac{\partial p}{\partial \mathit{\alpha }}}d\mathrm{\Omega }+\underset{\mathrm{\Omega }}{\int }\left[-\mathrm{\nabla }\stackrel{\mathrm{‾}}{\mathit{U}}\cdot \mathit{U}\right.\\ & -(\mathit{U}\cdot \mathrm{\nabla })\stackrel{\mathrm{‾}}{\mathit{U}}-\mathrm{\nabla }\cdot \left(\mathrm{2}{\mathit{\nu }}_{\mathrm{eff}}\mathbf{D}\left(\stackrel{\mathrm{‾}}{\mathit{U}}\right)\right)+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{C}_{D}A\left|\mathit{U}\right|\stackrel{\mathrm{‾}}{\mathit{U}}\\ & \left.+\mathrm{\nabla }\stackrel{\mathrm{‾}}{p}+{\displaystyle \frac{\partial J_{\mathrm{\Omega }}}{\partial \mathit{U}}}\right]{\displaystyle \frac{\partial \mathit{U}}{\partial \mathit{\alpha }}}d\mathrm{\Omega }=\mathrm{0}.\end{array}\end{array}$$

J_Γ and J_Ω, respectively, stand for the part of the cost function which is dependent on the flow state values on the boundary and volume of the domain. Due to the definition of the cost function (Eq. 13), its direct variation comes only from the interior domain $\left(\frac{\partial J_{\mathrm{\Gamma }}}{\partial \mathit{U}}=\mathrm{0},\frac{\partial J_{\mathrm{\Gamma }}}{\partial p}=\mathrm{0}\right)$. Moreover, it does not have any derivative of the pressure field $\left(\frac{\partial J_{\mathrm{\Omega }}}{\partial p}=\mathrm{0}\right)$. The only derivative of the cost function is with respect to the inflow and velocity in the interior of the domain and at the locations where the measurements are available. From the latter we have

$$\begin{array}{}\text{(20)}& {\displaystyle \frac{\partial J_{\mathrm{\Omega }}}{\partial \mathit{U}}}=-\mathrm{2}\left({\mathit{U}}_{{\mathrm{M}}_i}-{\mathit{U}}_{{\mathrm{S}}_i}\right)i=\mathrm{1},\mathrm{2},\mathrm{\dots }.\end{array}$$

Using Eq. (19) and (20) the adjoint equations can be derived as

$$\begin{array}{}\text{(21)}& {\displaystyle }\begin{array}{rl}-\mathrm{\nabla }\stackrel{\mathrm{‾}}{\mathit{U}}\cdot \mathit{U}& -(\mathit{U}\cdot \mathrm{\nabla })\stackrel{\mathrm{‾}}{\mathit{U}}=-\mathrm{\nabla }\stackrel{\mathrm{‾}}{p}+\mathrm{\nabla }\cdot \left(\mathrm{2}{\mathit{\nu }}_{\mathrm{eff}}\mathbf{D}\left(\stackrel{\mathrm{‾}}{\mathit{U}}\right)\right)\\ & +\left({\displaystyle \frac{\mathrm{2}}{{\mathit{\omega }}_i}}\right)\left({\mathit{U}}_{{\mathrm{M}}_i}-{\mathit{U}}_{{\mathrm{S}}_i}\right)-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{C}_{\mathrm{D}}A\left|\mathit{U}\right|\stackrel{\mathrm{‾}}{\mathit{U}},\end{array}\text{(22)}& {\displaystyle }\mathrm{\nabla }\cdot \stackrel{\mathrm{‾}}{\mathit{U}}=\mathrm{0},\end{array}$$

where ω_i is the volume of the cell in which the measurement is located.

3.2.1 Boundary conditions

The boundary integrals of Eq. (19) can be mathematically reformulated and reduced to

$$\begin{array}{}\text{(23)}& {\displaystyle }\underset{\mathrm{\Gamma }}{\int }[\stackrel{\mathrm{‾}}{\mathit{U}}\cdot \mathit{n}]{\displaystyle \frac{\partial p}{\partial \mathit{\alpha }}}d\mathrm{\Gamma }=\mathrm{0},\text{(24)}& {\displaystyle }\begin{array}{rl}\underset{\mathrm{\Gamma }}{\int }& \left[\mathit{n}\right(\stackrel{\mathrm{‾}}{\mathit{U}}\cdot \mathit{U})+{\mathit{\nu }}_{\mathrm{eff}}(\mathit{n}\cdot \mathrm{\nabla })\stackrel{\mathrm{‾}}{\mathit{U}}-\stackrel{\mathrm{‾}}{p}\mathit{n}]\cdot {\displaystyle \frac{\partial \mathit{U}}{\partial \mathit{\alpha }}}d\mathrm{\Gamma }\\ & -\underset{\mathrm{\Gamma }}{\int }\left[{\mathit{\nu }}_{\mathrm{eff}}(\mathit{n}\cdot \mathrm{\nabla }){\displaystyle \frac{\partial \mathit{U}}{\partial \mathit{\alpha }}}\cdot \stackrel{\mathrm{‾}}{\mathit{U}}\right]d\mathrm{\Gamma }=\mathrm{0},\end{array}\end{array}$$

where n is the unit normal vector from the boundary faces. Except for the inlet, which is the design space, the adjoint BCs should be chosen such that the above equations are held.

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 $\frac{\partial \mathit{U}}{\partial \mathit{\alpha }}$ is zero, the first integral of Eq. (24) is canceled. Then, the only way to satisfy the conditions,

$$\begin{array}{}\text{(25)}& {\displaystyle }& {\displaystyle }\stackrel{\mathrm{‾}}{\mathit{U}}\cdot \mathit{n}=\mathrm{0},\text{(26)}& {\displaystyle }& {\displaystyle }(\mathit{n}\cdot \mathrm{\nabla }){\displaystyle \frac{\partial \mathit{U}}{\partial \mathit{\alpha }}}\cdot \stackrel{\mathrm{‾}}{\mathit{U}}=\mathrm{0},\end{array}$$

is to apply a no-slip ($\stackrel{\mathrm{‾}}{\mathit{U}}=\mathrm{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 (($\mathit{n}\cdot \mathrm{\nabla })\frac{\partial \mathit{U}}{\partial \mathit{\alpha }}=\mathrm{0}$) and zero fixed pressure (p=0) are the common conditions for the primal system. These conditions fulfill Eq. (23) and cancel the second integral of Eq. (24). The only term that remains is the first term of Eq. (24), which needs to be zeroed out,

$$\begin{array}{}\text{(27)}& \left[\mathit{n}(\stackrel{\mathrm{‾}}{\mathit{U}}\cdot \mathit{U})+{\mathit{\nu }}_{\mathrm{eff}}(\mathit{n}\cdot \mathrm{\nabla })\stackrel{\mathrm{‾}}{\mathit{U}}-\stackrel{\mathrm{‾}}{p}\mathit{n}\right]\cdot {\displaystyle \frac{\partial \mathit{U}}{\partial \mathit{\alpha }}}=\mathrm{0}.\end{array}$$

After decomposition into tangent and normal components it can be shown that the relations below should hold

$$\begin{array}{}\text{(28)}& {\displaystyle }& {\displaystyle }\stackrel{\mathrm{‾}}{p}=\stackrel{\mathrm{‾}}{\mathit{U}}\cdot \mathit{U}+U_n{\stackrel{\mathrm{‾}}{U}}_n+{\mathit{\nu }}_{\mathrm{eff}}(\mathit{n}\cdot \mathrm{\nabla }){\stackrel{\mathrm{‾}}{U}}_n,\text{(29)}& {\displaystyle }& {\displaystyle }\mathrm{0}=U_n{\stackrel{\mathrm{‾}}{U}}_t+{\mathit{\nu }}_{\mathrm{eff}}(\mathit{n}\cdot \mathrm{\nabla }){\stackrel{\mathrm{‾}}{U}}_t,\end{array}$$

where subscripts n and t represent the normal and in-plane components, respectively. The adjoint BCs can be summarized as

$$\begin{array}{}\text{(30)}& {\displaystyle }\mathrm{ground}\left(\mathrm{wall}\right):\stackrel{\mathrm{‾}}{\mathit{U}}=\mathrm{0},\mathit{n}\cdot \mathrm{\nabla }\stackrel{\mathrm{‾}}{p}=\mathrm{0};\text{(31)}& {\displaystyle }\begin{array}{rl}\mathrm{top}/\mathrm{outlet}:& \stackrel{\mathrm{‾}}{p}=\stackrel{\mathrm{‾}}{\mathit{U}}\cdot \mathit{U}+U_n{\stackrel{\mathrm{‾}}{U}}_n+{\mathit{\nu }}_{\mathrm{eff}}(\mathit{n}\cdot \mathrm{\nabla }){\stackrel{\mathrm{‾}}{U}}_n,\\ & {\stackrel{\mathrm{‾}}{U}}_t=\mathrm{0}.\end{array}\end{array}$$

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. 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. (19) it is clear one needs to evaluate the following expression,

$$\begin{array}{}\text{(32)}& \begin{array}{rl}{\displaystyle \frac{\partial J}{\partial \mathit{\alpha }}}& ={\displaystyle \frac{\partial J}{\partial {\mathit{U}}_{\mathrm{inlet}}}}=\mathit{n}\left({\stackrel{\mathrm{‾}}{\mathit{U}}}_{\mathrm{inlet}}\cdot {\mathit{U}}_{\mathrm{inlet}}\right)\\ & +{\stackrel{\mathrm{‾}}{\mathit{U}}}_{\mathrm{inlet}}\left({\mathit{U}}_{\mathrm{inlet}}\cdot \mathit{n}\right)+\mathrm{2}{\mathit{\nu }}_{\mathrm{eff}}\mathit{n}\cdot \mathbf{D}\left({\stackrel{\mathrm{‾}}{\mathit{U}}}_{\mathrm{inlet}}\right),\end{array}\end{array}$$

to compute the sensitivity.

3.2.2 Wind direction effect

As it was mentioned before, in ABL CFD simulations it is common to simulate first a 1-D domain with a periodic boundary to obtain the inflow boundary condition. Then the cell center velocity of the 1-D run is copied directly to its counterpart boundary face in the 3-D domain. As a requirement, the number of cells in the 1-D mesh and faces in the vertical direction of the 3-D inflow boundary should be the same (see Fig. 1). Moreover, and ideally, the face center heights in the 3-D mesh are equal to their counterpart cell height in 1-D. Although, in current work, a circular inflow–outflow boundary is considered, with some small modification in the code the method can also be applied to other boundary shapes.

Figure 1 The inflow velocities of each cell from 1-D precursor run (a) are copied to the boundary of the 3-D domain (b), which has a similar number of cells in the vertical direction. Ideally the height of each face in the 3-D domain boundary is exactly the same as its counterpart cell in the 1-D mesh.

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,

$$\begin{array}{}\text{(33)}& {\displaystyle }& {\displaystyle }\mathrm{WD}=\mathrm{270}{}^{\circ }-\mathit{\theta },\text{(34)}& {\displaystyle }& {\displaystyle }{\left({\mathit{U}}_{\mathrm{inlet}}\right)}_x={\mathit{U}}_{\mathrm{1}\text{-}\mathrm{D}}\cos \left(\mathit{\theta }\right),\text{(35)}& {\displaystyle }& {\displaystyle }{\left({\mathit{U}}_{\mathrm{inlet}}\right)}_y={\mathit{U}}_{\mathrm{1}\text{-}\mathrm{D}}\sin \left(\mathit{\theta }\right),\text{(36)}& {\displaystyle }& {\displaystyle }{\left({\mathit{U}}_{\mathrm{inlet}}\right)}_z=\mathrm{0}.\end{array}$$

The differentiation of Eqs. (34) and (35) gives

$$\begin{array}{}\text{(37)}& {\displaystyle }{\displaystyle \frac{\partial {\left({\mathit{U}}_{\mathrm{inlet}}\right)}_x}{\partial {\mathit{U}}_{\mathrm{1}\text{-}\mathrm{D}}}}=\cos \left(\mathit{\theta }\right);{\displaystyle \frac{\partial {\left({\mathit{U}}_{\mathrm{inlet}}\right)}_x}{\partial \mathit{\theta }}}=-{\mathit{U}}_{\mathrm{1}\text{-}\mathrm{D}}\sin \left(\mathit{\theta }\right);\end{array}$$

$$\begin{array}{}\text{(38)}& {\displaystyle }{\displaystyle \frac{\partial {\left({\mathit{U}}_{\mathrm{inlet}}\right)}_y}{\partial {\mathit{U}}_{\mathrm{1}\text{-}\mathrm{D}}}}=\sin \left(\mathit{\theta }\right);{\displaystyle \frac{\partial {\left({\mathit{U}}_{\mathrm{inlet}}\right)}_y}{\partial \mathit{\theta }}}={\mathit{U}}_{\mathrm{1}\text{-}\mathrm{D}}\cos \left(\mathit{\theta }\right).\end{array}$$

The adjoint solver which was explained in the previous section computes the derivative of the cost function with respect to 3-D inflow velocities at each face of the boundary:

$$\begin{array}{}\text{(39)}& {\displaystyle \frac{\partial J}{\partial {\left({\mathit{U}}_{\mathrm{inlet}}\right)}_x}}{|}_{i,j};{\displaystyle \frac{\partial J}{\partial {\left({\mathit{U}}_{\mathrm{inlet}}\right)}_y}}{|}_{i,j};i=\mathrm{1},\mathrm{\dots },n;j=\mathrm{1},\mathrm{\dots },m,\end{array}$$

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 3-D circular boundary is $N=n\times m$.

Using the chain rule, one can compute the sensitivity with respect to each cell of the 1-D inflow velocity as follows

$$\begin{array}{}\text{(40)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\partial J}{\partial {\mathit{U}}_{\mathrm{1}\text{-}\mathrm{D}}}}{|}_i& =\left(\sum _{j=\mathrm{1}}^m{\displaystyle \frac{\partial J}{\partial {\left({\mathit{U}}_{\mathrm{inlet}}\right)}_x}}{|}_{i,j}\right){\displaystyle \frac{\partial {\left({\mathit{U}}_{\mathrm{inlet}}\right)}_x}{\partial {\mathit{U}}_{\mathrm{1}\text{-}\mathrm{D}}}}\\ & +\left(\sum _{j=\mathrm{1}}^m{\displaystyle \frac{\partial J}{\partial {\left({\mathit{U}}_{\mathrm{inlet}}\right)}_y}}{|}_{i,j}\right){\displaystyle \frac{\partial {\left({\mathit{U}}_{\mathrm{inlet}}\right)}_y}{\partial {\mathit{U}}_{\mathrm{1}\text{-}\mathrm{D}}}}\\ & =\left(\sum _{j=\mathrm{1}}^m{\displaystyle \frac{\partial J}{\partial {\left({\mathit{U}}_{\mathrm{inlet}}\right)}_x}}{|}_{i,j}\right)\cos \left(\mathit{\theta }\right)\\ & +\left(\sum _{j=\mathrm{1}}^m{\displaystyle \frac{\partial J}{\partial {\left({\mathit{U}}_{\mathrm{inlet}}\right)}_y}}{|}_{i,j}\right)\sin \left(\mathit{\theta }\right).\end{array}\end{array}$$

Eq. (40) means the x and y gradient components of the 3-D 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, Eq. (40) can be rewritten as

$$\begin{array}{}\text{(41)}& {\displaystyle \frac{\partial J}{\partial {\mathit{U}}_{\mathrm{1}\text{-}\mathrm{D}}}}={\displaystyle \frac{\partial J}{\partial {\left({\mathit{U}}_{\mathrm{3}\text{-}\mathrm{D}}\right)}_x}}\cos \left(\mathit{\theta }\right)+{\displaystyle \frac{\partial J}{\partial {\left({\mathit{U}}_{\mathrm{3}\text{-}\mathrm{D}}\right)}_y}}\sin \left(\mathit{\theta }\right).\end{array}$$

Using the same analogy, it can be shown that the sensitivity with respect to rotation angle can be obtained by

$$\begin{array}{}\text{(42)}& \begin{array}{rl}{\displaystyle \frac{\partial J}{\partial \mathit{\theta }}}& =-\left[{\displaystyle \frac{\partial J}{\partial {\left({\mathit{U}}_{\mathrm{3}\text{-}\mathrm{D}}\right)}_x}}\cdot {\mathit{U}}_{\mathrm{1}\text{-}\mathrm{D}}\right]\sin \left(\mathit{\theta }\right)\\ & +\left[{\displaystyle \frac{\partial J}{\partial {\left({\mathit{U}}_{\mathrm{3}\text{-}\mathrm{D}}\right)}_y}}\cdot {\mathit{U}}_{\mathrm{1}\text{-}\mathrm{D}}\right]\cos \left(\mathit{\theta }\right),\end{array}\end{array}$$

where the dot sign stands for the inner product of the two vectors.

4 Numerical results

The adjoint solver and Eq. (32) 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 2nd-order FD method in a simple 3-D 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.05 (m) and the turbulent eddies are modeled by standard k−ϵ model with canopy model by 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).

4.1 Gradient verification

A cylindrical domain with 1000 m radius and 300 m height is chosen. The mesh of the domain has 209 000 hexahedral cells in which the inflow–outflow boundary has 49 rows and 172 columns and a total number of 8428 faces ($N=\mathrm{49}\times \mathrm{172}=\mathrm{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⁻¹, respectively. The operational Reynolds number based on the free-stream velocity, forest area height and air kinematic viscosity is $R{e}_{\mathrm{h}}=\mathrm{4.8}\times {\mathrm{10}}^{\mathrm{5}}$.

Figure 2 The CFD-simulated velocity obtained from the reference simulation (WD = 270^∘): the general view of the domain including the cubic forest and the velocity field on a plane at the center of the domain, y=0, (a). The red line on the plane represents the location of desired target velocity profile (see Fig. 3). The velocity field on a plane at z=45 m above the ground (b).

For the primal CFD simulation of a circular domain the standard “inletOutlet” BC is used which checks if the flow is flowing into the 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.

The gradient evaluation is carried out for a simulation in which the inflow wind direction is WD = 240^∘. To have some reference wind speeds the terrain is simulated assuming the wind blows from west to east (WD = 270^∘). The domain and the velocity field on the planes y=0 m and z=45 m are shown in Fig. 2. The 1-D inflow boundary and the target velocity profiles in the domain are plotted in Fig. 3. The forest effect can be seen in the flow field and on the wind shear of the profile. Please note that the 30^∘ difference between these two simulations is not the step size for the finite-difference computation. The finite-difference step size for gradient validation of wind direction is 0.3^∘.

Figure 3 The 1-D inflow boundary (a) and the velocity profile and its selected target speeds in the domain (b).

The gradients obtained by the developed adjoint solver are plotted against the 2nd-order FD gradients in Fig. 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 and 100 m in which the maximum relative error is ϵ_rel=0.10. This difference can be traced back to grid resolution and the derivation of the adjoint equations and BCs which includes some simplifications such as the assumption of the frozen turbulence.

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

The wind direction gradients are tabulated in Table 1. The derivative of calibration cost function with respect to the change in wind direction is much higher than the change in 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 because it was shown in the previous section the sensitivity with respect to θ is obtained by mathematical operations only after the adjoint gradients are available. That said, the accuracy of the gradients computed by the adjoint solver and their signs show that they can be used for the purpose of the gradient-based calibration.

Figure 5 The calibration flowchart of the inflow boundary using the DAKOTA optimization package and the developed adjoint solver.

Table 1 Comparison of the wind direction gradients by finite differences (FD) and via the adjoint approach.

4.2 Inflow calibration

For the cases, studied in this section, the in-house terrainMesher software of the Fraunhofer IWES is used to generate the mesh. The primal flow field is simulated with an in-house CFD solver (Chang et al., 2018). The solver is a customized ABL-based version of the simpleFoam solver in the OpenFOAM package with a modified k−ϵ turbulence model which behaves like a standard model for the neutral condition.

To optimize the inlet velocity profile, the primal and adjoint solvers are coupled with the DAKOTA optimization package (Adams et al., 2017). The local gradient-based “CONMIN-frcg” solver of DAKOTA is used, which is based on the conjugate-gradient algorithm by Fletcher and Reeves (Reeves and Fletcher, 1964; Hager and Zhang, 2006). Starting from an initial guess, the algorithm updates the design parameters, αⁿ, using the recurrence of Eq. (1) in which

$$\begin{array}{}\text{(43)}& (\mathrm{\Delta }\mathit{\alpha }{)}^n=s^n{\mathbb{D}}^n\end{array}$$

and the positive step size sⁿ is obtained by a line search, and the directions 𝔻 are generated by the following rules:

$$\begin{array}{}\text{(44)}& \begin{array}{rl}{\mathbb{D}}^{n+\mathrm{1}}=& -g^{n+\mathrm{1}}+{\mathit{\beta }}^n{\mathbb{D}}^n;g^n={\left[\mathrm{\nabla }J\left({\mathit{\alpha }}^n\right)\right]}^T;\\ & {\mathit{\beta }}^n={\displaystyle \frac{{\left|g^{n+\mathrm{1}}\right|}^{\mathrm{2}}}{{\left|g^n\right|}^{\mathrm{2}}}};{\mathbb{D}}^{\mathrm{0}}=-g^{\mathrm{0}}.\end{array}\end{array}$$

Figure 5 shows the flow chart of the calibration, in which all the steps are followed sequentially. The optimizer starts with an initial guess (both inlet velocity and WD) and repeatedly asks either for the cost function value or the gradients. The primal and the adjoint solvers provide the required information, whenever it is needed, and this process continues until a certain convergence criterion is satisfied.

As a simplistic method, the inflow boundary of a RANS ABL domain can be represented by an empirical power law function or an analytically obtained logarithmic function, which is based on the Monin–Obukhov similarity theory (MOST) (Foken, 2006). During optimization, the optimizer may ask for a cost function evaluation with a new inflow boundary, which is highly unrealistic for an ABL domain, leading to poor numerical stability or even divergence. One may assume that the inflow boundary is an analytical empirical function and, instead of the inflow velocities, calibrate the parameters of that function. Having the gradient via adjoint solver and using the chain rule, the gradient of the cost function with respect to these parameters can be easily obtained. However, the inflow boundary of a real 3-D complex terrain is neither a power law nor a logarithmic function and such parametrization may fail.

Figure 6 Main dimensions of the test case with 1000 m length in the y direction. The velocity flow field on the plane at the center of the domain, y=0.

Figure 7 Optimization history (a) and inlet velocity profile comparison (b).

A simple approach is used in this study. This part is called the feasibility check in the optimization flow chart and is a Python script. First of all, the optimizer output is smoothed to avoid having any spikes in the inflow profile. Then it is checked whether a logarithmic, $f_{\mathrm{1}}\left(x\right)=A\ln (Bx+C)$, or a power law, $f_{\mathrm{2}}\left(x\right)=A(\frac{x}{B}{)}^C$, function could be fitted into it. If either of these functions is fitted and its coefficient of determination is above 0.96, the smoothed inflow from optimizer (not the fitted profile!) is accepted for the CFD solver. Otherwise, the optimization takes the last fitted profile and asks for a new gradient evaluation. In this way, the inflow boundary is not necessarily a logarithmic or power law profile, and, moreover, it is not so unrealistic to be problematic for the solver. As an alternative, constraints or a penalization term can be added to the objective function. This will be explored in future works when for instance the inflow turbulence properties are also considered design parameters.

4.2.1 Ishihara case

As a case study, an ABL domain with a 3-D hill at the center is considered (see Fig. 6). The hill has the shape $z=h{\cos }^{\mathrm{2}}(\sqrt{x^{\mathrm{2}}+y^{\mathrm{2}}}/\mathrm{2}L)$ with h=40 m and L=100 m. The scaled wind tunnel study of the case has been presented by Ishihara et al. (1999). The domain is meshed with 2×10⁶ hexahedral elements. The roughness value of the domain is set to be z₀=0.04 m. The operational Reynolds number based on the free-stream velocity, hill height and air kinematic viscosity is $R{e}_{\mathrm{h}}=\mathrm{1.5}\times {\mathrm{10}}^{\mathrm{4}}$.

Using Ishihara et al. (1999) wind tunnel measurements, the x component of the velocity over the hill (U_x) is used for the inflow velocity calibration. The velocity flow field at the center of the domain and the wake behind the hill can be seen in Fig. 6. Both the primal and the adjoint solvers have parallel scalability of OpenFOAM toolbox. The runtime of the adjoint solver is 60 %–70 % of the primal flow simulation, which has 30 min wall-clock time run with 24 CPU cores. It is worth noting that in the adjoint solver there is no adjoint turbulence equation to be solved.

The stopping criterion is such that the absolute error between the measurements and the simulated velocities over the hill is ϵ_abs<0.1 m s⁻¹. Figure 7 shows the history of optimization and the comparison of inlet velocity profiles. The optimization has converged with 14 primal and 12 adjoint calls. The optimal velocity profile is in good agreement with the experiment. There is a small deviation starting from the height $\frac{z}{h}>\mathrm{2}$. However, the comparison of the normalized velocity profiles over the hill, shown in Fig. 8, confirms that the optimized inlet boundary is able to reproduce the experiment velocity profiles.

Figure 8 Normalized vertical profiles of longitudinal velocity component on the central plane of the hill. U_h is the velocity at the hill height in the undisturbed region of the domain between inlet and hill.

Figure 9 The cylindrical domain of the Kassel terrain (a) and the velocity field at the site at 40 m perpendicular distance from the terrain's surface points, which is simulated with WD = 213^∘ (b). This simulation was used as the reference for the calibration.

4.2.2 Kassel case

For the second case study, the neutral condition of “Kassel Experiment” is considered. The domain, located near Kassel in Germany, is one of the cases of the New European Wind Atlas (NEWA) project (EU-ERA-NET(Accessed September 9, 2018). There are two meteorological masts at the site including a 200 m high mast (MM200). The wind rose of the site, provided by NEWA, indicates that most of the time the wind blows from the southwest (SW) to northeast (NE).

The site is represented by a cylindrical domain with 15 km radius and 4 km height (see Fig. 9). The structured mesh is generated with the Fraunhofer IWES in-house software “terrainMesher” with 80 cells in the vertical direction. A mesh independence study is conducted to verify the suitability of a mesh of 7×10⁶ hexahedral elements. The vicinity of the hill (z_hill=428 m) where the MM200 mast is installed consists of forested area. The leaf area density of the site is obtained from the airborne lidar data and is provided by the NEWA project (Dörenkämper, 2018).

Although the wind speed measurements of the MM200 mast from the site are available, there is not enough information for wind direction at certain heights. This information is necessary for the adjoint solver in which the difference between the components of the measured and simulated velocities is a force term on the right-hand side of the adjoint momentum equation (Eq. 21). Moreover, in the calibration of the solver with real measurements, it would become difficult to discuss the source of error when there is a discrepancy between the target and the calibrated profiles. This is because, aside from the inflow calibration process, which is the aim of the current study, many other parameters (e.g., turbulence model and the accuracy of forest and ground roughness map) are involved. Instead of using the real measurements of the mast, the velocity profile near to the mast from a reference simulation is considered. The selected wind speeds can be regarded as some pseudo-measurements.

Figure 10 Inflow boundary calibration for the Kassel domain; cost function convergence (a), wind direction history (b) and gradients history (c).

Figure 11 The target, initial and optimized 1-D inflow velocity (a). The velocity profile over the hill from which the pseudo-measurements are chosen (b).

The initial guess WD is defined as 270^∘, meaning wind blows from west to east. 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 of the optimization when the inlet boundary velocities have reached their optimum value. The adjoint solver runtime for this case is also 60 %–70 % of the primal solver wall-clock time, which is 33 min with 120 CPU cores.

The convergence criterion is defined such that the optimization stops when the absolute difference between simulated and measured value, ϵ_abs, at all heights is below a certain value. Figure 10 shows the history of optimization and gradients. The optimization has called for a total number of 41 primal solver runs for cost function evaluation. In addition eight 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 with respect to the WD is much bigger than with respect to inflow. The optimizer updates the WD based on the first gradient computation. This continues for a few iterations until the cost function value, instead of decreasing, increases. Then the optimizer calls for a new gradient computation. From that point onward both the inlet velocities and the wind direction are gradually updated. The optimum WD is found as 216^∘.

The initial and optimal results are compared in Fig. 11. The calibration is stopped based on a criterion, which is defined as ϵ_abs<0.2 m s⁻¹. Although there is a very small deviation between target and optimized inflow profiles, the velocity profile near to the MM200 mast is in good agreement with the pseudo-measurements at all five selected heights and its error is well below the accepted threshold. Here a couple of points should also be noted. Firstly, a tighter error criterion would increase the calibration iterations and subsequently the number of CFD solver calls. Secondly, the sensitivity of wind speed at a certain point in the domain with respect to a very small change in the inflow boundary is dependent on many parameters such as terrain complexity, wind direction, CFD model, etc., and cannot be easily generalized.

5 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 and using the reference wind speeds at certain heights inside of the domain, 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 two 3-D terrains for a neutral stratification condition: (1) the Ishihara test case and (2) the Kassel case. For the Ishihara test case, the optimal inflow was reached with 14 primal and 12 adjoint solver calls. The calibration of the inflow and its directions for the complex terrain of Kassel was achieved after 41 primal and 8 adjoint solver calls. In both cases, an absolute error between the measurements and the simulated velocities was used as a stopping criterion.

The main conclusion remarks of the study can be summarized as follows. (a) The developed calibration framework and the adjoint solver can be successfully applied to even complex domains. (b) The feasibility check of the optimizer output is crucial. Otherwise, at some point in the calibration process, the requested inflow leads to either poor performance or even complete failure (i.e., divergence) of the CFD solver. (c) The convergence criterion can have a big impact on the total number of solver calls. One possibility is to associate the criterion with the uncertainty of the measurements, which can be explored in future work. (d) The process can be further improved to reduce the number of CFD solver calls. For instance, a quasi-Newton optimizer (e.g., BFGS, the Broyden–Fletcher–Goldfarb–Shanno algorithm) or a better parametrization could be used. (e) The presented adjoint solver has the potential to be further developed by including Coriolis force, turbulence model and thermal stratification.