An improved inflow turbulence generator for large eddy simulation evaluation of wind effects on tall buildings

Abstract

Accurate inflow turbulence wind profile is a key premise of large eddy simulation (LES) methodology for super tall buildings. This paper presents an improved method for narrow band synthesis random flow generation to produce a more accurate inflow turbulence of LES based on superposition of harmonics. The proposed method considers time correlation by adopting time parameters, and an accurate expression of wind speed spectrum energy by a modified spectrum integral equation is used. These improvements could generate a more accurate turbulence intensity for the atmospheric boundary layer. A simple numerical example is applied to verify the spectrum characteristic, time correlation and spatial correlation of atmospheric numerical boundary layer. Meanwhile, the effects of the introduced parameters are investigated and the method is applied to the LES of the Commonwealth Advisory Aeronautical Research Council building (CAARC). The improved method could increase the accuracy of simulated turbulence intensity. Compared with the results of CAARC model obtained from the wind tunnel test, the improved method in this study decreases the errors of root-mean-square (rms) of the based bending moment from 5.3% to 4.1% in the along-wind direction and from 5.6% to 1.3% in the cross-wind direction.

KEYWORDS:

• Atmospheric boundary layer
• time correlation
• spectrum energy modification
• super tall building
• wind load

1. Introduction

LES technology has been widely applied to the turbulence flow simulation in the numerical wind tunnel test of tall buildings, owing to the computer performance enhancement and based on the development of a computational fluid dynamics method (Alminhana et al., 2018; Huang et al., 2007; Nozu et al., 2008; Xie et al., 2021; Zhang et al., 2022). The turbulence inflow characteristics, including the correlation of time and space and wind power spectrum, significantly affect the flow behavior in the computational domain for LES. Inappropriate inflow parameters could lead to longer development length required to obtain a fully developed turbulence flow. Setting appropriate parameters to satisfy the turbulence flow characteristics of atmospheric boundary layer (ABL) and generating accurate and efficient inflow turbulence wind profile are key problems of LES simulation. These concepts have also been the research focus in computational wind engineering in recent years.

Keating et al. (2004) classified turbulent flow generation methods into three categories: (1) precursor database method, (2) recycling method, and (3) synthetic turbulence method. The first two methods mainly divide the entire simulation of flow into two steps. In the generation step, the statistically steady wind profile is obtained, and the result is saved in a database for the calculating step (Liu & Pletcher, 2006; Lund et al., 1998). Both methods have the disadvantage of being computationally intensive, requiring large-scale servers and being computationally inefficient. The third method, namely, synthetic turbulence method is more commonly used for LES simulation. Huang et al. (2010) divided the synthetic turbulence method into three categories: harmonic synthesis (WAWS), vortex method, and random turbulence generation (RFG). The vortex method uses a random 2D vortex at the inlet and perturbations to a specified mean velocity profile to produce a spatiotemporally correlated velocity field that provides a reasonable mean turbulent kinetic energy distribution and average inflow dissipation rate. However, the target spectrum and statistical characteristics of turbulence, as well as the non-uniformity and anisotropy of turbulence are not well defined in the turbulence generation. The WAWS yields a turbulent velocity field that satisfies the target power spectrum and the cross-spectrum (Hoshiya, 1972; Iwatani, 1982; Kondo et al., 1997; Maruyama & Morikawa, 1994). However, the turbulence field obtained by this method is not dependent on the computational grid, and therefore, does not satisfy the continuity condition. The RFG was first proposed by Smirnov et al. (2001) and can produce turbulence fields that satisfy the Gaussian spectrum of the divergence-free turbulence field. However, this method does not consider spatial correlation. Huang et al. (2010) proposed a discrete synthetic turbulence with the spatial correlation method (DSRFG). Aboshosha et al. (2015) improved the target power spectrum based on the properties of the DSRFG method and spatial correlations, and then proposed a coherent discrete synthetic turbulence method (CDRFG). Yu et al. (2018) proposed a narrowband synthesis random inflow generation method (NSRFG) to address the inadequacies of the CDRFG and DSRFG methods in terms of low accuracy and computational efficiency of results. This method considers the spatial correlation of the turbulent field with more concise fluctuating wind field expressions and significantly improves computational efficiency. However, NSRFG also has some shortcomings, follow: the time dependence is not considered and the fluctuating wind speed field expressions are derived with insufficient rigorousness. The advantages and disadvantages of various stochastic turbulence methods are summarized in Table 1. In addition, Lund et al. (1998) first proposed the recycling rescaling method, which can only simulate flat terrain and can be applied to rough terrain after being improved by Nozawa and Tamura (2002).

Table 1. Random inflow turbulent generation method.

Random turbulence generation method	Advantages	Disadvantages
RFG	Time correlation	Only generates the wind speed field satisfaction. Does not consider spatial correlation
DSRFG	Any target spectrum	Low spatial correlation Low calculating efficiency Does not consider time correlation
CDRFG	Any target spectrum Good spatial correlation	Low calculating efficiency Not considering time correlation
NSRFG	Any target spectrum Good spatial correlation High calculating efficiency	Not considering time correlation

An improved inflow turbulence generation method, namely, the modified narrowband synthesis random inflow generation method (MNSRFG) is proposed based on the NSRFG method. A fluctuating wind field expression for the MNSRFG method is proposed, and a new expression introduces time parameters to consider the time dependence of the simulated wind speed. Moreover, a spectral energy correction factor in the spectral integration is introduced to accurately calculate the energy of the wind speed power spectrum. The improved method generates a fluctuating wind speed field that theoretically satisfies the ABL turbulence characteristics (divergence-free condition, time correlation, spatial correlation, and power spectral characteristics). The feasibility of the method is verified by numerical examples of typical ABL inflow turbulence. Then, the method was applied to a super high-rise building scale model (CAARC) for LES, and the results were consistent with those of the wind tunnel test (Yeo, 2011). The method can be applied successfully to the wind load simulations of supertall buildings.

2. Review of the NSRFG method

In the NSRFG method, the fluctuating wind speed field is given by the following simple harmonic function: (1) $u_i(x_j,t,f_n)=p_{i,n}\sin (2\pi f_{n}t+{\varphi }_n)$(1) where $u_i$ represents velocity in three directions (i = 1, 2, 3 for along-wind, across-wind, and vertical, respectively); j = 1, 2, 3 denotes x, y, and z directions, respectively; t is time; $f_n$ is frequency, and $p_{i,n}$ and ${\varphi }_n$ are the corresponding amplitude and phase of the trigonometric function, respectively. ${\varphi }_n\sim U(0,2\pi )$ indicates that ${\varphi }_n$ follows uniform distribution at the range of 0 to $2\pi$. n is the number of discrete points of the wind speed spectrum. The mean square value of the fluctuating velocity can be calculated by the following equation: (2) $\begin{array}{rl}{u}_{rms,i}^2(x_j,t,f_n)& =\underset{T\to \mathrm{\infty }}{lim}\frac{1}{T}{\int }_0^T{[p_{i,n}\sin (2\pi f_{n}t+{\varphi }_n)]}^{2}dt\\ & =\frac{p_{i,n}^2}{2}\end{array}$(2)

The power spectrum of the fluctuating wind speed in ABL satisfies the Karman spectrum, which is given by the following: (3) $u_{rms,i}^2(x_j,t,f_n)=S_{u,i}(f_n)\mathrm{\Delta }f$(3)

The following equation can be obtained from Equations (2) and (3): (4) $p_{i,n}=\sqrt{2S_{u,i}(f_n)\mathrm{\Delta }f}$(4) where $S_{u,i}(f_n)$ is the spectral value of the Karman spectrum at frequency $f_n$, and $\mathrm{\Delta }f$ is the bandwidth. Substituting Equation (4) into Equation (1), and considering spatial correlation, the fluctuating wind field expression for LES is as follows: (5) $u_i(x,t)=\sum _{n=1}^N\sqrt{2S_i(f_n)\mathrm{\Delta }f}\sin (k_{j,n}\tilde{x}+2\pi f_{n}t+{\varphi }_n)$(5) where N is the discrete number of power spectra with $k_{j,n}=\{k_{1,n}(k_{2,n})k_{3,n}\}$; to maintain the divergence-free condition of the inflow turbulent field, Equation (6) should be satisfied, as follows: (6) $\{\begin{array}{l}{P}_{1,n}{\displaystyle \frac{k_{1,n}}{L_{1,n}}}+P_{2,n}{\displaystyle \frac{k_{2,n}}{L_{2,n}}}+P_{3,n}{\displaystyle \frac{k_{3,n}}{L_{3,n}}}=0|{\mathit{k}}_n|=1\end{array}$(6)

Equation (6) shows that $k_{j,n}$ follows a uniform distribution on a spatial circular curve, the value of which can be calculated from the parametric equation of the spatial circular curve. In Equation (5), ${\tilde{x}}_{j,n}=x_j/L_{j,n}$, and $L_{j,n}$ is the turbulence spatial scale parameter; its expression is as follows: (7) $L_{j,n}=\frac{U_{av}}{f_n{c}_j{\gamma }_j}$(7) where $c_j$ and ${\gamma }_j$ are the attenuation and adjustment coefficients of the target spatial correlation in direction j, respectively, and ${\gamma }_j$ enables the resulting inflow turbulent field to satisfy the target spatial correlation condition. $U_{av}$ is the average value of $u_i(x,t)$.

3. Improving the NSRFG method

The spatial dependence of the wind speed expression of the NSRFG method is considered through Equation (7), but the temporal dependence is not considered. Therefore, the improved method introduces the time parameter ${\tau }_0$ and the wind speed expression is as follows: (8) $u_i(x,t)=\sum _{n=1}^N{p}_{i{,}_n}\sin (k_{j,n}\tilde{x}+2\pi f_n\frac{t}{{\tau }_0}+{\varphi }_n)$(8)

Meanwhile, the spectral energy expression of Equation (3) in the NSRFG method suffers from the following problem: the integration uses a summation approach, which is an approximate solution with some errors. To address this problem, the improved method introduces a spectral energy correction factor ${\beta }_i$ defined by the following: (9a) $\begin{array}{rl}{u^2}_{rms,1}& =(I_u{U}_{av}{)}^2={\beta }_1{\int }_0^{\mathrm{\infty }}{S}_u(f_n)df\Rightarrow {\beta }_1\\ & =\frac{{(I_u{U}_{av})}^2}{{\int }_0^{\mathrm{\infty }}{S}_u(f_n)df}\end{array}$(9a) (9b) $\begin{array}{rl}{u^2}_{rms,2}& =(I_v{U}_{av}{)}^2={\beta }_2{\int }_0^{\mathrm{\infty }}{S}_v(f_n)df\Rightarrow {\beta }_2\\ & =\frac{{(I_v{U}_{av})}^2}{{\int }_0^{\mathrm{\infty }}{S}_v(f_n)df}\end{array}$(9b) (9c) $\begin{array}{rl}{u^2}_{rms,3}& =(I_w{U}_{av}{)}^2={\beta }_3{\int }_0^{\mathrm{\infty }}{S}_w(f_n)df\Rightarrow {\beta }_3\\ & =\frac{{(I_w{U}_{av})}^2}{{\int }_0^{\mathrm{\infty }}{S}_w(f_n)df}\end{array}$(9c) where, ${\beta }_{\text{1}}$, ${\beta }_2$, and ${\beta }_3$ denote the spectral energy correction factor in along-wind direction, across-wind direction, and vertical-wind direction, respectively. $I_u$, $I_v$, and $I_w$ denote the turbulence intensity in along-wind direction, across-wind direction, and vertical-wind direction, respectively. Then, the expression in Equation (3) is transformed, as follows: (10) $\sum _{i=1}^3{u^2}_{rms,i}\text{(}x\text{, }t\text{, }{f}_n\text{)}=\sum _{i=1}^3\sum _{n=1}^N{\beta }_i{S}_i(f_n)\Delta f$(10) The expression of new $p_{i,n}$ is given by the following: (11) $p_{i,n}=\sqrt{2{\beta }_i{S}_{u,i}(f_n)\mathrm{\Delta }f}$(11)

Thus, the velocity expression for the LES fluctuating wind field of Equation (5) can be rewritten as follows: (12) $u_i(x,t)=\sum _{n=1}^N\sqrt{2{\beta }_i{S}_{u,i}(f_n)\mathrm{\Delta }f}\sin (k_{j,n}\tilde{x}+2\pi f_n\frac{t}{{\tau }_0}+{\varphi }_n)$(12) Substituting Equation (11) into Equation (6) yields the expressions for $k_{j,n}$, as follows: (13) $\{\begin{array}{l}{\text{k}}_{1,n}=-{\displaystyle \frac{q_{2,n}^2+q_{3,n}^3}{A_n}}\sin \theta {\text{k}}_{2,n}={\displaystyle \frac{q_{1,n}{q}_{2,n}}{A_n}}\sin \theta +{\displaystyle \frac{q_{3,n}}{B_n}}\cos \theta {\text{k}}_{3,n}={\displaystyle \frac{q_{1,n}{q}_{3,n}}{A_n}}\sin \theta -{\displaystyle \frac{q_{2,n}}{B_n}}\cos \theta \end{array}$(13) where $q_{i,n}=p_{i,n}/L_{i,n}$, $A_n=\sqrt{{(q_{2,n}^2+q_{3,n}^2)}^2+q_{1,n}^2{q}_{2,n}^2+q_{1,n}^2{q}_{3,n}^2}$, $B_n=\sqrt{q_{2,n}^2+q_{3,n}^2}$, and $\theta$ is uniformly distributed with a mean of 0 and a variance of $2\pi$. The parameter $k_{j,n}$ can be calculated directly from Equation (13) without solving the nonlinear set of equations similar to the CDRFG methods (Aboshosha et al., 2015), thereby improving the computational efficiency of the inflow turbulence.

4. Validation of the MNSRFG method

4.1. Wind speed power spectrum characteristic

The wind profile at the boundary layer of the inlet atmosphere is selected in accordance to the urban terrain described in the Load Code for The Design of Building Structures: GB50009-2012 (2012) with the roughness of the landform $\alpha =0.22$. No specification on winds is available in the incoming cross direction (u) and the vertical direction (w) of the turbulence degree. Thus, the following expressions are provided based on the international ESDU specification: ${\sigma }_v/{\sigma }_u=1-0.22{\cos }^4(\pi /2\cdot z/h)$ and ${\sigma }_w/{\sigma }_u=1-0.44{\cos }^4(\pi /2\cdot z/h)$ (Characteristics of Atmospheric Turbulence Near the Ground, Part II. Single Point Data for Strong Winds (Neutral Atmosphere), 1993). The detailed parameters of ABL are shown in Table 2. Figure 1 (a) shows the time series of fluctuating velocity in three directions at a typical point at 1 m height. Figures 1(b)–(d) show the power spectrum in three directions at the same position, consistent with the Karman spectrum expressed in equations (14), indicating that the power spectral characteristics of the fluctuating wind speed time range of the turbulent field generated by the proposed method satisfy the requirements of ABL turbulence. $L_u$, $L_v$, and $L_w$ denote the turbulence integral scale in along-wind direction, across-wind direction, and vertical-wind direction, respectively. (14a) $\begin{array}{rl}{S}_u(f)& =\frac{4{(I_u{U}_{av})}^2(L_u/U_{av})}{{[1+70.8{(f{L}_u/U_{av})}^2]}^{5/6}}\end{array}$(14a) (14b) $\begin{array}{rl}{S}_v(f)& =\frac{4{(I_v{U}_{av})}^2(L_v/U_{av})[1+188.4{(2f{L}_v/U_{av})}^2]}{{[1+70.8{(2f{L}_v/U_{av})}^2]}^{11/6}}\end{array}$(14b) (14c) $\begin{array}{rl}{S}_w(f)& =\frac{4{(I_w{U}_{av})}^2(L_w/U_{avg})[1+188.4{(2f{L}_w/U_{av})}^2]}{{[1+70.8{(2f{L}_w/U_{av})}^2]}^{11/6}}\end{array}$(14c)

Figure 1. Simulated wind speed time history and power spectrum, (a) Simulated wind speed time history, (b) u, (c) v, (d) w.

Table 2. Parameters of ABL.

Parameters	Equations
Wind speed profile	$U_{av}=U_{\text{ref}}(z/H_{\text{ref}}{)}^{\alpha }$, where $U_{\text{ref}}$ is the reference wind speed at $H_{\text{ref}}$ height and z is the height
Turbulence intensity profile	$I_u(z)=I_{10}(z/10{)}^{-\alpha }$, $I_v\text{(}z\text{) = }{I}_u(z){\sigma }_v/{\sigma }_u$,$I_w\text{(}z\text{) = }{I}_u(z){\sigma }_w/{\sigma }_u$
Turbulence integral scale	$L_u(z)=300(z/300{)}^{0.46+0.074\ln z_0},L_v(z)=0.5({\sigma }_v/{\sigma }_u{)}^3{L}_u(z),L_w(z)=0.5({\sigma }_w/{\sigma }_u{)}^3{L}_u(z)$ where $z_0$ is equal to 0.7 m

4.2. Spatial correlation

Satisfying the spatial correlation is one of the basic requirements for simulating the turbulence characteristics of the atmospheric boundary layer. The spatial correlation of the simulation-generated turbulence field requires the three velocity components u, v, and w to be validated in the x, y, and z directions, respectively. To test the spatial correlation of the generated turbulent field, intervals were selected in each of the x, y, and z directions at a height of 1 m with total space of 2 m, consisting of 20 points with interval of 0.1 m. The monitoring points are shown in Figure 2. Parameters in Equation (7) are set as follows: C₁= 5, C₂= 8, and C₃= 10 and γ₁= 3.2, γ₂= 1.6, and γ₃= 1.4. Figure 3 shows the spatial correlation of u, v, and w components in three directions, and the objective function is adopted from the spatial correlation function proposed by Hémon and Santi (2007). Figure 3 shows that the spatial correlation of u, v, and w velocities in three directions is generally consistent with the objective function.

Figure 2. Monitoring points in x, y, and z directions.

Figure 3. Spatial correlation of u, v, and w components in three directions.

4.3. Time correlation

The time correlation is defined as follows: (15) $R(m\delta \tau )=\frac{1}{M-m}\sum _{j=0}^{M-m}u(j\delta \tau )u[(j+m)\delta \tau ]$(15) where m is an integer such that ${\tau }_m=m\delta \tau$ and $0\le m<M$, $\delta \tau$ is the time step, and M is the length of the vector ${\tau }_m$. The target function of the time correlation is defined as follows: (16) $R_i(\tau )=e^{-|\tau |/T_i}$(16)

The time scale is defined as follows: (17) $T_i={\int }_0^{\mathrm{\infty }}{R}_i(\tau )d\tau \equiv \sum _{j=0}^{M_0}{R}_i(j\delta \tau )\delta \tau$(17) where M₀ is less than M. As the time lag tends to infinity, low-frequency fluctuations result in time-dependent close to zero. In the absence of a sufficient upper limit, the time scale cannot be approximated by Equation (17) (Castro & Paz, 2013). In this study, the time scale is proposed to be calculated by setting M₀ corresponding to the first zero value point of the time dependence function. Figure 4 shows the time dependence of the velocity in three directions. The target function of time is defined in Equation (17). The enlarge plots at the range from 0 to the value of $\tau$, where correlation function is equal to zero, are shown to compare the discrepancies. T_u, T_v, and T_w are equal to 0.0545, 0.0174, and 0.0087, respectively. Figure 4 shows that the coincidence is excellent in the direction of the incoming flow u, which has the greatest influence on the wind load, and the time dependence is almost the same as the objective function. Although the v and w directions have some discrepancies, the general trend is consistent. Thus, they are considered to satisfy the time-dependent characteristics of atmospheric boundary layer turbulence.

Figure 4. Time correlation (the right plots are the enlarge plots at the range from 0 to the value of $\tau$, where the correlation function is equal to zero). (a) u direction, (b) v direction, (c) w direction.

To further analyze the effect of time parameter ${\tau }_0$ on the time scale $T_i$, the objective function is estimated according to Taylor’s assumptions, i.e. $T_i={\text{L}}_i/U_{avg}(i=u,v,w)$. ${\tau }_0$ is a non-dimensional time parameter and denote the value corresponding to best time correlation at the given computed range of $\tau$ in Figure 5. The figure presents the mean and standard deviation values of the time scales on a sample of 70 velocity time series. The time parameter ${\tau }_0$ can be modified using different values in the MNSRFG method, i.e. the MNSRFG method can be varied to obtain a range of time scales from which the appropriate values can be selected, and then the temporal dependence is better satisfied and closer to the actual atmospheric turbulence. Within the analyzed range of 0.75–1.5, selecting the time scale values closer to the target object without harming other physical characteristics (target turbulence velocity spectrum, spatial correlation, and pulsed RMS values of the simulated turbulent flow) is possible.

Figure 5. Statistical analysis of the influence of ${\tau }_0$ on the time scale.

4.4. Modified coefficient of wind speed power spectrum

The effect of the size of the frequency interval on the RMS value of the velocity time series is analyzed to validate the spectral energy correction factor proposed in Equation (10). Theoretically, as the spectrum becomes finer discretized, the energy content in each frequency should be included in the velocity time series. Table 3 shows the results of the analysis, indicating that as $\mathrm{\Delta }f$ becomes progressively smaller, the MNSRFG method converges to the target value with smaller errors than the NSRFG method. In addition, the comparison of the turbulence intensities generated by the two methods in the three directions (Figure 6) indicates that the MNSRFG method generates turbulence intensities that are more consistent with the simulation, and the reduced values of turbulence intensities are less. The target values of turbulence intensities are defined in Table 3. Thus, the proposed MNSRFG method for simulating the atmospheric boundary layer inflow turbulence can produce velocity time series that better match the target values of the physical problem into closer to the actual atmospheric boundary layer turbulence.

Figure 6. Comparison of turbulence intensity in three directions.

Table 3. Comparison of the standard deviation of the velocity components of the two methods.

	$\mathrm{\Delta }f$	10	5	2	1	0.2	Target
σu	NSRFG	1.08	1.25	1.31	1.31	1.31	1.35
	Relative error	−20.0%	−7.5%	−2.8%	−2.8%	−2.8%
	MNSRFG	1.34	1.34	1.35	1.35	1.35
	Relative error	0.5%	0.4%	0.1%	0.1%	0.1%
σv	NSRFG	1.00	0.99	0.99	0.99	0.99	1.05
	Relative error	4.6%	5.6%	6.0%	6.0%	6.0%
	MNSRFG	1.04	1.06	1.06	1.05	1.05
	Relative error	1.2%	0.7%	0.6%	0.02%	0.02%
σw	NSRFG	0.65	0.65	0.65	0.65	0.65	0.74
	Relative error	12.3%	12.4%	12.4%	12.4%	12.4%
	MNSRFG	0.72	0.73	0.74	0.74	0.74
	Relative error	3.3%	1.0%	0.3%	0.3%	0.3%

5. Numerical example

The MNSRFG method is applied to the LES of the wind flow around a CAARC standard supertall building model. The results are compared with the wind tunnel test results, with the results of the NSRFG method.

5.1. CAARC model and wind tunnel test

The CAARC model, which is the standard model of supertall buildings commonly used in the wind engineering industry, was used as the object of the analysis. The parameters of the CAARC model are shown in Table 4. The dimensions of the considered model are 45.72 m × 30.48 m × 182.88 m, as shown in Figure 7 (a). The model consists of 45 layers, and the height of each layer is 4.0 m with a mass of 223016 kg. The first three orders of the model frequency were 0.20, 0.27, and 0.40 Hz, and the damping ratio is 5% for the wind displacement calculation. In accordance with the Chinese Building Code, the Structural Load Code (Load Code for The Design of Building Structures: GB50009-2012, 2012), a rigid model with a model ratio of 1:300 to the CAARC standard model was tested in the pressure measurement wind tunnel. It requires the simulation of a Class C geomorphic wind field in the model area. The model is arranged with seven layers of measurement points, with 20 measurement points in each layer, as shown in Figure 7(b). Figure 7(c) shows the test model in the wind tunnel.

Figure 7. Parameters of the CAARC model. (a) 3D model, (b) Pressure taps, (c) Wind tunnel test of CAARC model.

Table 4. Parameters of the CAARC model.

Model dimensions (m)		Layer number	Height of each layer (m)	Weight of each layer (kg)	Frequency (Hz)		Damping ratio
Length	45.72	45	4	223016	First	0.02	5%
Breadth	30.48				Second	0.27
Height	182.88				Third	0.40

5.2. Numerical model

5.2.1. Parameter setting

Large eddy flow simulation of the CAARC standard model was performed using UDF programming linked to the commercial software Fluent 15.0, considering the computer’s computing power and minimum grid size for simulation at the model scale. The computational domain size and boundary conditions are set with reference to Aboshosha et al. (2015) and Tominaga et al. (2008), as shown in Figure 8(a). The model has a blocking rate of 1.7% that is less than 3%, and the watershed is set to satisfy the blocking rate requirement (Kim et al., 2013). The pressure-implicit with splitting of operators (PISO) algorithm was adopted for the numerical iterative solution. Time discretization is in the second-order implicit format, and spatial discretization is in the bounded center differential format. In addition, the wall adaptive local eddy-viscosity (WALE) model is used as the sub-grid model of the LES (Yu et al., 2018). The expression of Reynolds number is as follow: $Re=\rho VD/\nu$. Where $\rho$ is the air density, V is the freestream velocity, D is the length scale, and $\nu$ is the coefficient of dynamic viscosity of air. In this study, the value of the freestream velocity used in the Reynolds number is 11.1 m/s. The length scale used in the Reynolds number is the width of building model and the value of it is equal to 0.1016 m. The parameters in the CFD match that in the experiment. Therefore, the Reynolds number in the CFD that is 7.7 × 10⁴ is equal to that in the experiment. The time step is selected as 0.002 s and the maximum Courant–Friedrichs–Lewy number is approximately 2, satisfying the requirements of computational stability and convergence. The total duration of the numerical simulation calculation is 12 s, which is sufficient to ensure that the wind flow field is adequately developed. The last 10 s of data for the analysis of the results was obtained and converted to 621 s to satisfy the 10 min requirement. Numerical simulations were performed on a high-performance computer server with 12-core parallel computing, the convergence criterion for iterative computing is defined as all variables, and the individual massless residuals are less than 10⁻⁵.

Figure 8. Calculation of watershed size and corresponding boundary conditions and computational grids. (a) Calculation of watershed size and corresponding boundary conditions. (b) Computational grids for empty domain. (c) Computational grids for wind tunnel test of CAARC model.

The time correlation parameter of the MNSRFG method is set to 0.85, and the spatial correlation parameter is assumed to be C₁= 5, C₂= 8, and C₃ = 10 and γ₁= 3.2, γ₂ = 1.6, and γ₃ = 1.4. The calculation efficiency of the MNSRFG method is insignificantly different from the NSRFG method. The grid independence test is performed using grid calculations with different densities, and the minimum cell size is set to 0.03 times of the width of the windward side of the building. Two grid schemes, G1 and G2, were adopted to achieve the grid independency for the computational results. The smallest grid sizes near the model surface were selected as 0.03B for grid G1 and 0.06B for grid G2 to balance the computational resources available and the numerical efficiency. Approximately 6.0 million 3D structural grid cells for grid G1 and 3.5 million cells for grid G2 were generated for the simulation of empty domain and the CAARC model, as shown in Figure 8 (only grid scheme G1 is shown in this figure). The simulation results from two grid schemes were very close. The value of the y + parameter of the grid used is at the range of 20–140. The cell size growth factor is 1.05, the total number of grid elements is approximately 6 million, and the grid division is shown in Figures 8(b) and (c).

5.2.2. Equilibrium wind field validation

The self-sustainability of the flow field simulated by the proposed MNSRFG method is verified using an empty watershed consistent with the original building model with the grid size unchanged. Sixteen monitoring points were placed vertically at intervals of 0.1 m from the bottom to the top of the empty basin, in line with the center of placement of the original building model. Those points are set every 30 m at prototype scale corresponding to 0.1 m at the model scale along the height. The highest point is at the height of 480 m that is higher than the gradient wind height. Velocity characteristics simulated under an empty basin with velocity timescale obtained from the last 10 s of the simulated timescale data are used to derive the mean wind profile and turbulence intensity profile, as shown in Figure 9. The figure shows that the wind profile simulated by the proposed MNSRFG method has better self-retaining flow field and turbulence intensity, and the results are closer to the predicted target values than the results of the conventional NSRFG method, with a 3% improvement in turbulence simulation accuracy at the top of the building location.

Figure 9. Comparison of the mean wind speed and turbulence intensity profile. (a) Mean wind speed and (b) Turbulence intensity.

5.2.3. Wind speed field and wind pressure coefficient distribution

Figure 10 shows the results of the instantaneous velocity field at 8 s for the CAARC standard model computed from LES. The results include the instantaneous velocity contour plots and traces of the turbulent fields in the horizontal section (at the 1/2H building height) and the mid-perpendicular (at the 1/2B building width), as shown in Figure 10. The figure shows that when an incoming wind load acts on the surface of a building, it separates on both sides of the building, creating a separation vortex in the incoming wind. This is mechanism of the cross-wind response of supertall buildings. It also shows an inconstant vortex shedding on the sides and top. Figure 11 also shows that the traces clearly depict the details of the complex turbulent structures at different scales in the vicinity of the building.

Figure 10. Flow field: overall contour plots of instantaneous velocity in the 8 s. (a) Plan sectional view at 1/2H building height. (b) Vertical section view at 1/2B building width.

Figure 11. Flow field: a partially enlarged contour plots of the instantaneous velocity in the 8s. (a) Enlarged plan sectional view at 1/2H building height. (b) Enlarged vertical section view at 1/2B building width.

Figure 12. Mean wind pressure coefficient. (a) Wind tunnel test, (b) NSRFG and (c) MNSRFG.

Figure 13. Fluctuating wind pressure coefficient. (a) Wind tunnel test, (b) NSRFG, (c) MNSRFG.

Figure 14. Mean wind pressure coefficient and RMS of wind pressure coefficient at 2/3 H of the building model. (a) Mean wind pressure coefficient. (b) RMS of wind pressure coefficient.

Figures 12 and 13 compare the mean wind pressure and fluctuating wind pressure coefficient distribution diagrams of the CAARC standard model obtained from the MNSRFG method with those obtained from the NSRFG method and the wind tunnel test. Overall, the wind pressure distributions of NSRFG and MNSRFG methods are in agreement with the data of the wind tunnel test except that those near the base of the building. On the windward side and leeward side near the base of building, the mean pressure and fluctuating wind pressure are all slightly less than those of wind tunnel test. For the sidewall, the values of mean wind pressure are negative. The absolute values of two LES method are also less than that of wind tunnel test. The same conclusion can be drawn for the fluctuating wind pressure distribution on the sidewall. The discrepancies are possibly related to the grid space near the ground, the wall model, and the formation of a horseshoe vortex around the base of the building. The flow near the base of building is complex. For the detailed analysis, a layer of pressure coefficient data for measurement points at 2/3 height of the CAARC model (Figure 7) is shown in Figure 14 for comparisons with the results of the wind tunnel test and the NSRFG method. The mean wind pressure results of the MNSRFG method and the NSRFG method agree well with the wind tunnel test results in terms of size, but the mean wind pressure distributions at the leeward side of both methods are slightly smaller than those of the wind tunnel test. The accuracy of the improved MNSRFG method is higher than that of the conventional NSRFG. However, for the fluctuating wind pressure, Figure 14 shows the differences between the numerical simulation results of the two methods and the wind tunnel test results. Fluctuations of flow field pressure due to non-constant flow can lead to fluctuating wind pressure producing numerical errors when the synthetic turbulence is used in the LES for flow conditions (Kim et al., 2013). In addition, the mechanism of fluctuating wind pressure action is much more complex than the mean wind pressure, and the numerical errors are related to the simulated turbulence characteristics, sub-grid models, wall functions, setting time step, numerical algorithm, discrete scheme, and convergence criteria. Figure 14 also shows that for fluctuating wind pressure, the results of the MNSRFG method are much closer to the wind tunnel test results than those of the NSRFG method on the windward side and sidewall. However, it should be noted that the accuracy of the MNSRFG method is less than that of the NSRFG method on the leeward side. That will result in the reduction of the accuracy of MNSRFG method on the base moment. Fortunately, the large fluctuating values do not occur at leeward side. In particular, the error of the mean square value of the fluctuating wind pressure coefficient for the sidewall measurement points with large fluctuating values was reduced from 18%–29% (NSRFG) to 6%–12% (MNSRFG).

5.2.4 Base moments and wind-induced response

Along-wind base moment coefficient C_My, cross-wind base moment coefficient C_Mx, and base torsional coefficient C_Mz are defined as follows: (18) $\begin{array}{rl}& C_{M_y}=\frac{M_y}{\frac{1}{2}\rho V_h^{2}B{H}^2},C_{M_x}=\frac{M_x}{\frac{1}{2}\rho V_h^{2}D{H}^2},\\ & C_{M_z}=\frac{M_z}{\frac{1}{2}\rho V_h^{2}BDH}\end{array}$(18) where V_h is the reference wind speed, B and D are the widths in the along-wind and cross-wind directions of the building model, respectively, and H is the height of the building model. Figure 15 shows the power spectrum of the basal bending moment coefficient, from which the main characteristics of the power spectrum of the base bending moment in the along-wind and cross-wind directions simulated by MNSRFG and NSRFG methods, respectively, agree well with the wind tunnel test. The trend of the power spectrum obtained from numerical simulation is also in good agreement with the wind tunnel test results. However, compared with the NSRFG method, the power spectrum of the base moment obtained by the MNSRFG method combined with the LES agrees better with the wind tunnel test. In addition, the frequency corresponding to the highest peak of the power spectrum of the cross-wind-to-base moment work factor is the frequency of the building’s vortex shedding, which is an important parameter for evaluating the cross-wind aerodynamic characteristics of high-rise structures. Figure 15 shows that the power spectral peaks of the cross-wind base moment from the LES using the MNSRFG method have good coincidence with the wind tunnel test results, indicating that the MNSRFG method can predict the vortex shedding frequency of building models more accurately for structural design. In the torsional direction, the simulation results are relatively different from the wind tunnel test results and not as good as those in the two other directions. Table 5 compares the RMS values, the relative error of the two methods, and the mean value with those of the wind tunnel test. The table also shows that the improved MNSRFG method has a better performance in the along-wind direction, cross-wind direction, and torsional direction than that of the NSRFG method. The error of the base moment coefficient decreased from 5.3%, 5.6%, and 11.1% (NSRFG) to 4.1%, 1.3%, and 6.6% (MNSRFG), respectively. The mean values of the MNSRFG and NSRFG method are approximate to those of the wind tunnel test.

Figure 15. Comparison of power spectrum of base moment coefficient. (a) Along-wind direction, (b) Cross-wind direction, (c) Torsion.

Table 5. Comparison of values of the base moment coefficient.

	CMx			CMy			CMz
Method	RMS	error	Mean	RMS	error	Mean	RMS	error	Mean
Wind tunnel test	0.1206	/	0.0156	0.1351	/	0.6372	0.0422	/	0.0073
MNSRFG	0.1191	1.3%	0.0035	0.1407	4.1%	0.6426	0.0394	6.6%	0.0079
NSRFG	0.1138	5.6%	0.0020	0.1279	5.3%	0.6370	0.0375	11.1%	0.0052

Figure 16 shows the maximum displacement response of the top of the structure at different wind speeds. It indicates that in the along-wind and cross-wind directions, the results of the simulation calculations by the MNSRFG method are closer to the wind tunnel test than those of the NSRFG method. The improvement in the along-wind direction is better than that in the cross-wind and torsional directions. When the incoming wind speed is 35 m/s, the maximum displacement response errors in the along-wind, cross-wind, and torsional direction reduce from 21.1%, 19.2%, and 19.3% (NSRFG) to 7.6%, 17.3%, and 14.1% (MNSRFG). In general, the trends of the results in the three directions are more consistent, and the results of the wind response calculations are generally consistent with those of the wind tunnel tests.

Figure 16. Wind-induced vibration response. (a) Along-wind direction, (b) Across-wind direction, (c) Torsion.

6. Conclusion

Simulation of ABL turbulence by the improved narrowband synthesis leads to the following conclusion:

1. The MNSRFG method, which considers time dependence and spectral energy correction, is capable of producing simulations that better match the target value of the physical problem. The velocity time series of the atmospheric boundary layer is used to more accurately simulate the turbulence characteristics. The generated turbulence field satisfies divergence-free condition, time correlation, spatial correlation, and power spectral characteristics.

2. The MNSRFG method improves the accuracy of the wind pressure simulation on the CAARC scalar surface, similar to the NSRFG method. Compared with the fluctuating wind pressure generated by NSRFG, the MNSRFG method reduces the error from 18%–29% to 6%–12%.

3. The MNSRFG method can significantly improve the basal bending moment simulation over the previous NSRFG method in the along-wind, transverse, and torsional directions. The basal moment coefficient errors decreased from 5.3%, 5.6%, and 11.1% to 4.1%, 1.3%, and 6.6%, respectively. The cross-wind vortex shedding frequency obtained by the MNSRFG method agrees better with the wind tunnel test than that of the NSRFG method.

4. The analysis of the structural wind vibration response shows that the improvement in the along-wind response is better relative to the cross-wind and torsional response. Compared with the NSRFG method, the error in the displacement response can be reduced from 21.1% to 7.6% at a wind speed of 35 m/s.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

The authors would like to acknowledge the support from the National Natural Science Foundation of China [grant number 51925802 and 51808153], the Natural Science Foundation of Guangdong Province [grant number 2021A1515011824], the 111 Project [grant number D21021], and the Guangzhou Municipal Science and Technology Project [grant number 20212200004].