Large‐eddy simulation of plume dispersion in a turbulent boundary layer flow generated by a dynamically controlled recycling method

When conducting large‐eddy simulations (LESs) of plume dispersion in the atmosphere, crucial issue is to prescribe time‐dependent turbulent inflow data. Therefore, several techniques for driving LESs have been proposed. For example, in the original recycling (OR) method developed by Kataoka and Mizuno (Wind and Structures, 2002, 5, 379–392), a mean wind profile is prescribed at the inlet boundary, the only fluctuating components extracted at the downstream position are recycled to the inlet boundary. Although the basic turbulence characteristics are reproduced with a short development section, it is difficult to generate target turbulent fluctuations consistent with realistic atmospheric turbulence. In this study, we proposed a dynamically controlled recycling (DCR) method that is a simple extension of the OR procedure. In this method, the magnitude of turbulent fluctuations is dynamically controlled to match with the target turbulent boundary layer (TBL) flow using a turbulence enhancement coefficient based on the ratio of the target turbulence statistics to the computed ones. When compared to the recommended data of Engineering Science Data Unit (ESDU) 85020, the turbulence characteristics generated by our proposed method were quantitatively reproduced well. Furthermore, the spanwise and vertical plume spreads were also simulated well. It is concluded that the DCR method successfully simulates plume dispersion in neutral TBL flows.


| INTRODUCTION
In local-scale atmospheric dispersion at short distances from an emission source, an important issue is the accurate prediction of airborne radionuclide materials released from nuclear facilities for safety and consequence assessments.For analyzing atmospheric dispersion of a plume released from a stack over complex terrains, wind tunnel experiment is used as a rational tool for incorporating complicated effects from local topography and/or building effects.Many experiments have been conducted on plume dispersion in turbulent boundary layer (TBL) flows from a long time ago.For example, Robins (1978) extended the wind tunnel technique of simulating the atmospheric boundary layer flow to the dispersion of a plume released from a point source and demonstrated that the plume spreads obtained from a wind tunnel are similar in magnitude to those in the atmosphere.Fackrell and Robins (1982) conducted wind tunnel experiments of concentration fluctuations of a plume by using a high-response hydrocarbon analyzer and investigated the characteristics of mean and fluctuating concentrations, peak concentrations, and turbulent fluxes of concentrations.
On the other hand, with rapid development of computational technology, computational simulations using large-eddy simulation (LES) has been evolving as an alternative to wind tunnel experiments and widely used in various fields for fundamental and application research purposes.When applying LES techniques to plume dispersion in the atmosphere, inflow turbulence conditions have to be carefully prescribed in order to simulate plume dispersion in the atmosphere.
The inflow turbulence methods have been proposed by many researchers over the past decades.There are typically two methods: synthetic method and recycling method.These were reviewed by Tabor and Baba-Ahmadi (2010), Wu (2017), andZhong et al. (2021).In the former, Okaze and Mochida (2017) and Xie and Castro (2008) developed the synthetic method based on the Cholesky decompositions of turbulence flux tensor.The main advantage is that the target TBL flows can be obtained by prescribing the turbulence statistics such as turbulent fluxes, spatial and time correlations.
In the latter, Lund et al. (1998) developed the rescaling-recycling method.They produced a spatiallydeveloping TBL flow by rescaling the vertical profile of the wind velocities at a downstream station and recycling it to the inlet.In this method, the length of the development section is usually 10δ.In our previous study (Nakayama et al., 2008), we conducted LESs of plume dispersion using Lund's method and showed that the simulation results of mean concentration profiles are reproduced well in comparison to the wind tunnel experimental data of Fackrell and Robins. Kataoka and Mizuno (2002) simplified it and proposed the recycling method.In this original recycling (OR) method, only the fluctuating wind components extracted at a downstream station are recycled to the mean flow at the inlet.This approach has an advantage in easily generating basic TBL flows with a short development section of usually 2δ.However, it is difficult to quantitatively reproduce the target turbulence intensities because the turbulence statistics cannot be pre-specified at the inlet condition.Therefore, we attempted to generate target TBL flows by using the recycling method and various roughness obstacles such as a tripping fence and roughness blocks over a long fetch as well as wind tunnel experiments (Nakayama et al., 2013(Nakayama et al., , 2014)).Although this approach generated neutral TBL flows well, it requires a long development section for quantitively obtaining target turbulence intensities.
In this study, we incorporated a turbulence enhancement coefficient into the OR method to efficiently reproduce target turbulent fluctuations.Our objective is to conduct LESs of plume dispersion in TBL flows using the dynamically controlled recycling (DCR) method that is a simple extension of the OR procedure and examine the performance.

| THE DCR METHOD
The OR method proposed by Kataoka and Mizuno (2002) is formulated as follows: Here, u and U are the instantaneous and mean winds of the streamwise component x ð Þ.The suffixes of y, z, t, inlt, and recy are spanwise and vertical directions, time, the inlet boundary, and the recycle station, respectively.ϕ z ð Þ and u recy are a damping function and the spanwiseaveraged winds at the recycle station.The turbulent fluctuations for each component obtained by the difference between u recy and u recy are recycled to the inlet.In this study, we incorporate a turbulence enhancement coefficient into the OR method as follows: Here, the β is the turbulence enhancement coefficient for dynamically controlling turbulent fluctuations to match with the target turbulence statistics.σ u_tgt is the target standard deviation of wind velocities for the streamwise component given from wind tunnel or the ESDU (Engineering Science Data Unit, 1985) recommended data, etc. σ u_recy z,t ð Þ is the spanwise-averaged standard deviation of the wind velocity at the recycle station.The turbulent fluctuations for each component are dynamically scaled to the target standard deviation by driving the coefficient β.

| SIMULATION SETTINGS
We conduct LESs of plume dispersion in the neutral TBL flows generated by the two different methods.These results are compared with the existing wind tunnel experimental data.

| Wind tunnel experimental dataset
To evaluate our proposed method for predicting atmospheric dispersion of a plume released from a tall stack, we compared to mainly the experimental data of Fackrell and Robins (1982).The experiments were carried out in the large boundary layer wind tunnel.This wind tunnel has a 24 m long test section, a cross section with 9.1 m in width and 2.7 m in height.In the experiment, the neutral TBL flow was simulated by a barrier and surface roughness.The scale of the modeled boundary layer was 1:550.
Two passive plumes were released from an elevated and a ground-level source.The elevated source height was 0.23 m that corresponds to about 125 m in full scale.The streamwise variations of vertical and lateral plume spreads at the release heights were measured.

| Numerical model
LESs were conducted using LOHDIM-LES (Nakayama & Takemi, 2020;Nakayama et al., 2013Nakayama et al., , 2014)).The coupling algorithm of the velocity and pressure fields is based on the MAC method with the Adams-Bashforth scheme for time integration.The Poisson equation is solved by the SOR method.For the spatial discretization, a second-order-accurate central difference is used.The subgrid-scale turbulent effect is represented by the standard Smagorinsky model with a constant of 0.1 (Smagorinsky, 1963).

| Computational model
The computational domain size and grid resolution are 3600 m Â 800 m Â 600 m and 10 m Â 10 m Â (5-30) m in the streamwise, spanwise, and vertical directions, respectively.
The boundary conditions are: time-dependent turbulent inflow data at the inlet; a zero-gradient condition at the outflow; a free-slip condition for the horizontal components and a zero-speed condition for the vertical component at the top; the Monin-Obukhov similarity theory (Monin & Obukhov, 1954) at the bottom; and a periodic condition at the spanwise boundaries.

| Computational conditions
The vertical profile of the mean inflow U z ð Þ is given by the power law exponent α of 0.14 estimated from the roughness length z 0 by the following equation (Counihan, 1975).
Here, z 0 and U z ð Þ are 0.05 m and has wind speed of 10 m s À1 at 10 m height, respectively.The height of the lower part of the simulated TBL δ is 300 m.The recycle station is set to 2δ downstream position from the inlet.The vertical profiles of the target standard deviations of wind velocities estimated from the ESDU data are shown in Figure 1.These are given to σ tgt z ð Þ in Equation (3).The time step interval was 0.02 s.The total length of the simulation run was 40 min.The β z,t ð Þ for each component at each height were set to 1.5 to effectively generate turbulent fluctuations for the first 10 min.The plume was continuously released 20 min after the simulation start.The statistics values of the flow and dispersion fields were computed over the last 10 min.Because the

| Flow field
Figure 2 shows time series of the β z, t ð Þ at heights of 50, 125, and 200 m.Those parameters largely vary for about 5 min just after the dynamically calculating start.Then those become nearly constant 20 min after the calculating start, which indicates that turbulent equilibrium states are nearly achieved.In particular, the β z, t ð Þ for v-component at 50 m height is much larger than the target value even in the turbulent equilibrium state because the spanwise turbulent intensity generated by the OR method is smaller than the target one.This quantitative difference between the calculated and target values produces the large β z,t ð Þ.We discuss this matter below.Figure 3 shows instantaneous fields of the TBL flows generated by the OR and DCR methods.It is found that the turbulent motions simulated by the DCR method are more active than the OR method.As Simens et al. (2009) pointed out, the eddy persistence can induce artificial periodicity if the reference plane is located too close to the inlet boundary in the case of using recycling type method.Such artificial periodicity can be clearly seen in both instantaneous velocity fields because of a short development section of 2δ.Figures 4 and 5  Both mean wind profiles are matched with the power law profile.The turbulence intensities for each component generated by the OR method are much smaller than the curve of the ESDU data for smooth surface case although those show the same profiles at each downstream position.On the other hand, the streamwise turbulence intensities by the DCR method are distributed between the curves for the smooth and moderately rough surfaces while the spanwise turbulence intensities are distributed along the one for the smooth surface.The vertical turbulence intensities are larger than the curve for the moderately rough case at the source position of x=δ ¼ 2. However, those become to be distributed between the curves for the smooth and moderately rough surfaces at the downstream position from the source.The Reynolds stresses of the OR method is also much smaller than the experimental data of Fackrell and Robins (1982) as shown in Figure 6.On the other hand, those of the DCR method are largely improved although those are still underestimated a little.
Figure 7 shows normalized autocorrelation for the streamwise wind velocities at the recycle station.In particular, it is clearly observed that the autocorrelation at 200 m height in the DCR method show positive sharp peaks with a time period of nearly 2δ=U z ð Þ. Figure 8 shows power spectral density of the streamwise wind velocities S. The results in both methods show sharp peaks at frequency f corresponding to U z ð Þ=2δ especially in the lower frequency range.These tendencies are due to the artificial periodicity induced by a recycling type method with a short development section.The energy content of the spectrum in the DCR method is a little larger than that in the OR method in the higher frequency range.
One of the primary drawbacks of the recycling method is the emergence of artificial correlations.However, the DCR method quantitatively reproduces the target turbulence of the TBL well for simulating plume dispersion in the atmosphere by dynamically controlling the turbulence enhancement coefficients.

| Dispersion field
Figure 9 shows the instantaneous fields of plume dispersion in the two TBLs.Since the turbulence intensities generated by the DCR method are quantitatively reproduced well, the plume spreads become larger than those in the OR method up to large downstream distances from the point source.
Figure 10 shows a comparison of the spanwise and vertical plume spreads of the LES with the Pasquill-Gifford (P-G) chart (Turner, 1970) and the wind tunnel experimental data.The most common stability classification scheme is the P-G curves.The curves of atmospheric stability classes A-F indicate highly unstable, moderately unstable, slightly unstable, neutral, moderately stable, and extremely stable, respectively.The spanwise plume spread of Fackrell and Robins (1982) has a tendency to gradually become smaller with a downstream distance from the point source while the vertical one is distributed between the curves of the classes C and D. On the other hand, both plume spreads of Michioka et al. (2011) are distributed between the two curves up to large downstream distances.
Since the turbulence intensities by the OR method are considerably underestimated, both plume spreads are much smaller than the curve of the stability class D. On the other hand, the spanwise plume spread obtained by the DCR method is a little larger than the curve D near the point source and become smaller than that with a downstream distance.The vertical one is overestimated and is nearly distributed along the curve C near the downstream position of the point source.This is due to the overestimation of the vertical turbulence intensities at the source position.However, that becomes to be distributed between the curves of the classes C and D with a downstream distance from the source.Michioka et al. (2011) generated large-scale turbulent motions using an active grid installed at the front of the test section in the wind tunnel and concluded that the spanwise plume spread is enhanced by the lateral large-scale turbulent motions.It seems from their experiments that the underestimation of the spanwise plume spread is due to the technical problem of representing the lateral large-scale turbulent motions in the LESs.In future work, we plan to incorporate the technique for generating the lateral largescale turbulent motions into the DCR method.However, both plume spreads are considerably improved in comparison to the OR method case.Furthermore, the distribution patterns of the plume spreads are quantitatively the same as the wind tunnel experimental data of Fackrell and Robins (1982).It is concluded from these results that our proposed DCR method practically simulates plume dispersion in target TBL turbulences with a short development section of 2δ.

| CONCLUSIONS
We developed the DCR method using the turbulence enhancement coefficient and applied to LESs of plume dispersion in the simulated TBL flows.In this method, the magnitude of turbulent fluctuations is dynamically controlled to match with the target turbulence statistics using the coefficient based on the ratio of the target values to the computed ones.
First, the LES results of the TBL flows generated by the OR and DCR methods were compared with the recommended data of ESDU 85020.In the former method, the turbulence intensities are much smaller than the ESDU data.In the latter method, the target turbulence intensities are quantitatively reproduced well.The Reynolds stress profiles of the DCR method are quantitatively reproduced well in comparison to the wind tunnel experimental data.Although one of the primary drawbacks of the recycling method is the emergence of artificial correlations, it is shown that the DCR method effectively generates the target TBL turbulences with a short development section of 2δ.
Then, we applied to LESs of plume dispersion in the two TBL flows.In the former method, the spanwise and vertical plume spreads were considerably smaller than the P-G chart.In the latter method, those are improved and are matched with the wind tunnel experimental data well although the spanwise plume spread is underestimated especially at a large downstream distance.This is due to the technical problem of representing the lateral large-scale turbulent motions in the LESs.From these results, it is concluded that our proposed method successfully simulates plume dispersion in the target TBL flows.
For predicting plume dispersion under real meteorological conditions, it is promising to couple with a mesoscale meteorological simulation (MMS) model used for prescribing the input conditions.Presently, many researchers have been studying on various coupling techniques using inflow turbulence generation methods.Because there is a dependency of a development section length on generated turbulence scales, the crucial issue is to fill the gap between large-scale meteorological disturbances simulated by MMS models and small-scale turbulences generated by LES models (Nakayama et al., 2012).Typical coupling approaches are to determine the inflow boundaries by the MMS outputs linearly interpolated on the LES model grids.These approaches reasonably reproduce mean wind flow fields under realistic meteorological conditions.However, the turbulence kinetic energy (TKE) in the MMS models is not fully utilized to quantitatively maintain the magnitude of turbulent fluctuations in LES models coupled to them.Our proposed method has an advantage in quantitatively reproducing realistic turbulent fluctuations by driving the turbulence enhancement coefficient based on the TKE in MMS models.

F
I G U R E 1 Target of the standard deviation of wind velocities for each component.stack height in nuclear facility is usually greater than 100 m, we set to the release height of 125 m and compared with the experimental data for the elevated source case.
compare mean flows and turbulence characteristics obtained by the two different methods.The ESDU provides comprehensive turbulence characteristics of a neutrally stratified atmospheric boundary layer based on independent field measurements ranging from the ground surface to 300 m.It recommends vertical profiles of turbulence intensities for each wind component and their relationship depending on surface roughness.The z 0 of 0.01, 0.1, and 0.3 m shown in Figure 5 corresponds to the curves for smooth, moderately rough, and rough surface cases, respectively.

F
I G U R E 2 Time series of the turbulence enhancement coefficients.

F
I G U R E 3 Instantaneous fields of the turbulent boundary layer flows generated by the two different methods.F I G U R E 4 Vertical profiles of mean wind velocities.

F
I G U R E 6 Vertical profiles of Reynolds stress.F I G U R E 7 Normalized autocorrelation of the streamwise wind velocities at the recycle station.F I G U R E 8 Power spectral density of the streamwise wind velocities at the recycle station.

F
I G U R E 9 Instantaneous field of the plume dispersion simulated by the two methods.F I G U R E 1 0 Streamwise variation of spanwise and vertical plume spreads at the stack height.