CFD Model Verification and Aerodynamic Analysis in Large-Scaled Venlo Greenhouse for Tomato Cultivation

Abstract

To address the challenges of climate change and food security, the establishment of smart farm complexes is necessary. While there have been numerous studies on the productivity and environmental control of individual greenhouses, research on greenhouse complexes is considerably limited. Conducting environmental studies during the design phase of these complexes poses financial constraints and practical limitations in terms of on-site experiments. To identify potential issues that may arise when developing large-scale greenhouse complexes, it is possible to utilize modeling techniques using computational fluid dynamics (CFD) to assess environmental concerns and location issues before constructing the facilities. Consequently, simulating large-scale CFD models that incorporate multiple greenhouses and atmospheric conditions simultaneously presents significant numerical challenges. The objective of this study was to design and verify the 3D CFD model for a large-scale Venlo greenhouse, where acquiring field data before construction is not feasible for designing a greenhouse complex. The verification of the CFD models was conducted using the improved grid independence test (GIT) and wall Y⁺ approaches. The findings revealed that a grid resolution of 0.8 m and a first-layer height of 0.04 m were suitable for developing large-scale greenhouse models, resulting in a low Root Mean Square Error (RMSE) of 3.9% and a high coefficient of determination (R²) of 0.968. This process led to a significant reduction of 38% in the number of grid cells. Subsequently, the aerodynamic characteristics and regional ventilation efficiency were analyzed in a 3D greenhouse model for developing a new large-scale greenhouse complex.

1. Introduction

The recent issue of climate change has led to an increase in the annual global-average air temperature and unpredictable precipitation patterns [1]. As a result, open-field cultivation has become cost-inefficient, unreliable, and contributes to an inconsistent food production system [2,3]. The rise in population has further exacerbated food insecurity challenges in numerous countries worldwide [4]. To tackle these challenges, greenhouse cultivation has emerged as an effective solution [5]. Under the unfavorable climate conditions, it is necessary to establish large-scale protected agricultural facilities to meet the increasing population pressure on food production systems. The Korean government recently initiated a project to develop large-scale greenhouse complexes on reclaimed land aimed at cultivating and exporting high-quality horticultural products year-round [6]. Additionally, research has been conducted on atmospheric conditions and greenhouse construction on Korean reclaimed lands [7,8].

To develop a large-scale greenhouse complex for commercial purposes, it is crucial to identify and understand environmental factors. The microclimate factors inside the greenhouse have a significant impact on crop productivity, and it is important to ensure the stability and uniformity of the microclimate, which can be controlled by the airflow pattern based on the ventilation system. The ventilation system in greenhouse designs aims to optimize the internal microclimate in a greenhouse such as the temperature, humidity, light, and carbon dioxide levels while also considering energy efficiency [9,10,11,12]. Common greenhouses used in Korea utilize natural ventilation. Natural ventilation greenhouses are sensitive to external environmental factors, such as wind direction, wind speed, temperature, rainfall, etc. To study the microclimate inside the greenhouse, it is essential to research the ventilation efficiency caused by the airflow pattern from the external environment in the naturally ventilated structure. Ventilation efficiency is influenced by the greenhouse’s structure and orientation, external environmental conditions, and the presence and type of crops. As a result, many studies have been conducted to evaluate the effectiveness of natural ventilation on the internal microclimate of greenhouses [13,14,15,16,17,18,19,20].

Greenhouse design studies have utilized both field experiments and computation fluid dynamics (CFD) modeling and simulations to investigate the greenhouse microclimate and structural safety in relation to environmental problems [21,22,23,24]. In the greenhouse structural design phase, evaluating ventilation efficiency by installing facilities on-site is time-consuming and costly, making it practically impossible. Especially when setting up large-scale facilities, it is necessary to apply simulation technologies such as CFD to overcome the challenges of on-site monitoring. Modeling and simulating a large greenhouse system offers an alternative approach to address challenges associated with field experiments [25,26]. However, validating these large-scale models poses difficulties due to the unavailability of field data from such structures before construction. Moreover, in order to achieve acceptable model accuracy, the computational grid must be sufficiently fine to minimize computational errors [27]. Consequently, modeling large-scale structures with a fine computational grid results in a large volume of computational grid cells, presenting significant limitations in computational capability and time.

The CFD modeling research process includes modeling, discretization, grid verification, field validation processes, data acquisition, selection of simulation regime, turbulence models, and reporting of results [28]. To utilize the CFD model, verification by checking grid design fidelity and validation by comparing field-monitored data were crucial. Many studies have carried out the verification of the CFD model. Specifically for greenhouses, they are typically verified by comparing the quantitative and qualitative patterns of simulation model predictions with on-site measured environmental data, mainly the temperature [26,29,30]. However, previous models have mainly focused on small-scale greenhouses of less than 0.1 hectare, making them less suitable for developing design and control strategies for large-scale greenhouses and greenhouse complexes considering the internal microenvironment of greenhouses. And procedural guidelines are not enough for specifically addressing the design and validation of large-scale CFD models for agricultural structures, where the scale is a limiting factor for both modeling and data acquisition processes.

The main objective of this study was to aerodynamically analyze a large-scale greenhouse before construction. A CFD model was designed by an accuracy verification procedure using the grid independence test and wall boundary conditions. By utilizing the verified CFD model, a three-dimensional, 2-hectare Venlo greenhouse was designed and simulated to investigate crucial factors such as internal microclimate distributions and ventilation efficiency based on airflow patterns corresponding to the external environment.

2. Materials and Methods

2.1. Numerical Simulation

The computational fluid dynamics (CFD) method enables the explicit calculation of airflow patterns through the numerical solution of the Reynolds-averaged form of the Navier–Stokes equations [19,31]. The CFD model was used to analyze the microclimate and airflow pattern in the greenhouse by aerodynamic approaches. The equations governing the mass, momentum, energy, and concentration in the transient state are expressed as follows:

$$\left(\frac{𝜕\varnothing }{𝜕\mathrm{t}}\right)+\frac{𝜕}{𝜕{\mathrm{x}}_{\mathrm{j}}}\left(u_j\varnothing \right)=\frac{𝜕}{𝜕{\mathrm{x}}_{\mathrm{j}}}\left({\Gamma }_{\varphi }\frac{𝜕\varnothing }{𝜕{\mathrm{x}}_{\mathrm{j}}}\right)+s_{\varnothing },$$

(1)

where $\varnothing$ is the variable of interest, the three velocity components $u_j\left(\mathrm{m}\times {\mathrm{s}}^{-1}\right)$, the temperature T (K), and the specific humidity ${\Gamma }_{\varphi }$ and $s_{\varnothing }$ represent the diffusion coefficient and source term of $\varnothing$. The details of the equations were described in [32]. A previous study carried out the determination of the most suitable turbulence model for greenhouse CFD modeling [7]. They reported the Renormalization-group (RNG) k-epsilon turbulent model to yield the highest level of accuracy, with an R² value of 0.99. In this study, the RNG k-epsilon model in conjunction with enhanced wall functions was selected which was expressed using Equations (2) and (3):

$$\mathsf{\rho }\left(\frac{\mathrm{d}\mathsf{\kappa }}{\mathrm{d}\mathrm{t}}\right)=\frac{𝜕}{𝜕{\mathrm{x}}_{\mathrm{i}}}\left[{\mathsf{\alpha }}_{\mathrm{k}}{\mathsf{\mu }}_{\mathrm{e}\mathrm{f}\mathrm{f}.}\frac{𝜕\mathrm{k}}{𝜕{\mathrm{x}}_{\mathrm{i}}}\right]+{\mathrm{G}}_{\mathrm{k}}+{\mathrm{G}}_{\mathrm{b}}-\mathsf{\rho }\mathsf{\epsilon }-{\mathrm{Y}}_{\mathrm{M}},$$

(2)

$$\mathsf{\rho }\left(\frac{\mathrm{d}\mathsf{\kappa }}{\mathrm{d}\mathrm{t}}\right)=\frac{𝜕}{𝜕{\mathrm{x}}_{\mathrm{i}}}\left[{\mathsf{\alpha }}_{\mathrm{k}}{\mathsf{\mu }}_{\mathrm{e}\mathrm{f}\mathrm{f}.}\frac{𝜕\mathrm{k}}{𝜕{\mathrm{x}}_{\mathrm{i}}}\right]+{\mathrm{G}}_{\mathrm{k}}+{\mathrm{G}}_{\mathrm{b}}-\mathsf{\rho }\mathsf{\epsilon }-{\mathrm{Y}}_{\mathrm{M}},$$

(3)

The Boussinesq model was employed to account for the effect of gravity. This was achieved by including the buoyancy force resulting from air density variations as a source term in the momentum equation [33]. The simulation software used for this study was Ansys Fluent (ver. 2022 R2, Ansys Inc., Rochester, New York, NY, USA) which utilizes a pressure-based finite-volume discretization code. The crops canopy in greenhouses has a great impact on natural ventilation efficiency and the distribution of the microenvironment. The effect of the crop canopy was modeled by considering the crop domain as porous media with the momentum loss effect. The Darcy–Forchheimer equation was used to introduce the drag force induced by the crops canopy into the CFD model as a source term, shown in Equations (4)–(7) [23,26,31,34,35]. For flow characteristics through the plant canopy, the nonlinear quadratic term of Equation (4) was considered. The momentum source term can also be expressed in terms of the unit volume of the canopy as in Equation (5). Equating Equations (4) and (5) and considering the nonlinear momentum loss term, the momentum source term was reduced, as in Equation 6. The porous media characteristics of the crop canopy were described by the inertial resistance factor ($c_2$) given by Equation (7) [36].

$${\mathrm{S}}_{\varnothing }=-\left(\frac{\mathsf{\mu }}{{\mathrm{K}}_{\mathrm{p}}}\ast \mathrm{u}+\frac{{\mathrm{C}}_{\mathrm{F}}}{\sqrt{{\mathrm{K}}_{\mathrm{P}}}}\ast \mathsf{\rho }\ast {\mathrm{u}}^2\right),$$

(4)

$${\mathrm{S}}_{\varnothing }=-{\mathrm{I}}_{\mathrm{L}\mathrm{A}\mathrm{v}}\mathsf{\rho }{\mathrm{C}}_{\mathrm{D}}{\mathrm{u}}^2=-\frac{{\mathrm{C}}_{\mathrm{F}}}{\sqrt{{\mathrm{K}}_{\mathrm{P}}}}\ast \mathsf{\rho }\ast {\mathrm{u}}^2,$$

(5)

$${\mathrm{S}}_{\varnothing }=-{\mathrm{I}}_{\mathrm{L}\mathrm{A}\mathrm{v}}{\mathrm{C}}_{\mathrm{D}}{\mathrm{u}}^2,$$

(6)

$${\mathrm{c}}_2=\frac{2\ast {\mathrm{C}}_{\mathrm{F}}}{{\left({\mathrm{K}}_{\mathrm{p}}\right)}^{0.5}},$$

(7)

where $\mathsf{\rho } \left(\mathrm{k}\mathrm{g}\times {\mathrm{m}}^{-3}\right)$ is the fluid density, $\mathsf{\mu } \left(\mathrm{k}\mathrm{g}\times {\mathrm{s}}^{-1}{\times \mathrm{m}}^{-1}\right)$ is the dynamic viscosity, ${\mathrm{K}}_{\mathrm{p}}\left({\mathrm{m}}^2\right)$ is the permeability of the porous media, and ${\mathrm{C}}_{\mathrm{F}}$ is the nonlinear momentum loss term.

2.2. Target Greenhouse Design

The target area for this study was the Saemangeum-reclaimed land (35°46′28.34″ N, 126°40′2.66″ E) located between Gunsan (E 126°45′, N 36°00′) and Buan (E 126°42′, N 35°43′) in Korea (Figure 1). In this region, among the crops mainly produced in large greenhouse facilities, we designed the facility targeting tomatoes considering the crop canopy, ventilation requirement, and greenhouse design. The target greenhouse type was multi-span Venlo, accounting for 80% of the recent smartfarm systems in Korea. Also, Venlo greenhouses have a relatively larger eave height of about 6 m which can accommodate a large range of crop heights, such as tomato and paprika. The weather conditions of the target area used in this study were the average weather conditions for Buan and Gunsan.

2.3. Atmospheric Boundary-Layer (ABL) Design

The impact of the external environment on the internal microenvironment of the greenhouse was considered by designing the wind characteristics in the immediate vicinity of the greenhouse. The size of the computational domain representing the external environment was determined based on recommendations from [35,37,38]. For the CFD models, the overall height of the computational domain was set to 10 times the height (H) of the model greenhouse, while the horizontal extent was set to 15 times. The authors of [7] determined the aerodynamic surface characteristics of reclaimed land in Korea. These findings, combined with the 11-year average weather data for the summer, were used to develop wind and turbulence profiles using the log law method (Equation (8)) described by [39], and the specific equations involved in this method are represented by Equations (8)–(10).

$$\mathrm{U}\left(\mathrm{z}\right)=\frac{{\mathrm{u}}^{*}}{\mathrm{k}}\ln \frac{\mathrm{z}+{\mathrm{z}}_0}{{\mathrm{z}}_0},$$

(8)

$$\mathrm{k}\left(\mathrm{z}\right)=\frac{{\mathrm{u}}^{*2}}{\sqrt{{\mathrm{C}}_{\mathsf{\mu }}}.},$$

(9)

$$\mathsf{\epsilon }\left(\mathrm{z}\right)=\frac{{\mathrm{u}}^{*3}}{\mathrm{k}\left(\mathrm{z}+{\mathrm{z}}_0\right)},$$

(10)

where U(z) ($\mathrm{m}\times {\mathrm{s}}^{-1}$) is the vertical wind speed, z (m) is the reference height, $u^{*}$ is the shear velocity $(\mathrm{m}\times {\mathrm{s}}^{-1})$, k is the Karman constant, $z_0$ is the roughness height (m), and $C_{\mu }$ is the empirical constant. The turbulent kinetic energy $\epsilon \left(z\right) \left({\mathrm{m}}^2\times {\mathrm{s}}^{-3}\right)$, wind profile, and turbulent kinetic energy dissipation rate (k) for reclaimed lands were designed using the surface roughness factor. The designed profiles were applied as inlet boundary conditions to the model using compiled User-Defined Functions (UDFs) in the main CFD code, Ansys Fluent software. The parameters and values used in the design are summarized in

Table 1. Wind profile turbulent kinetic energy and dissipation rate design parameters.

Parameter	Value
Karman constant (k)	0.42
Roughness length (Z0)	0.03 m
Empirical constant	0.09
Reference height (${\mathrm{z}}_{\mathrm{R}}$)	10 m
Summer average wind speed	3.08 m·s−1
Frictional velocity (u*)	0.227 m·s−1

.

2.4. Model Verification Theories

To enhance the accuracy of the CFD model, we conducted grid verification using the GIT (grid-independent test) and yY⁺ theory. The GIT is a crucial step in validating CFD simulations. It involves the systematic refinement of the numerical grid used to discretize the fluid domain and comparing the resulting solution with a reference solution or experimental data. The goal is to determine the minimum grid resolution required to obtain a converged and accurate solution. This is to ensure that the numerical solution is accurate and reliable while minimizing computational costs [28,40]. The need for a grid independence test arises from the fact that the numerical solution of the governing equations of fluid flow is based on discretizing the domain into a finite set of grid points. The accuracy of the solution depends on the resolution of the grid. A coarse grid may lead to numerical errors and inaccuracies in the solution, while a fine grid may result in excessive computational costs and unnecessary detail [41]. Several methods have been used to perform a grid independence test: the Richardson extrapolation method, the grid convergence index (GCI), and the grid refinement study are the commonest methods. In turbulence models, viscous fluid flow at larger Reynolds numbers, the fluid’s inertia overcomes the viscous stresses, and the flow becomes unsteady, termed as turbulent flow. The CFD model was designed in the range of low Reynolds numbers, specifically between 2750 < Re < 4970. All flows of practical engineering interest are virtually turbulent [42]. Several studies have also reported the flow inside greenhouses to be typically turbulent. The turbulent flow velocity profile near the wall has three distinct layers: the viscous sublayer, the log layer, and the defect layer. The log layer is the portion of the boundary layer where the sublayer and defect layer merge and the law of the wall accurately represents the velocity. The turbulence vanishes near the walls due to the no-slip boundary condition for the velocity as well as the blocking effect caused by the wall. In the viscous sublayer, an adequate numerical resolution of a solution is required which calls for a very fine mesh. This is because of the thinness of the sublayer and high solution gradients. This makes calculations time-consuming and may be impractical for large-scale industrial application [43]. Fine numerical mesh and wall functions approaches are commonly employed in the CFD calculation of the turbulent flow near walls to resolve the thin near-wall sublayer. The wall functions approach requires less computational effort and is thus strongly favored by industrial calculations. However, their performance is often poor because of inappropriate implementation and partly because the schemes themselves have inherent limitations [44].

Wall Y⁺ is a dimensionless number that determines the flow region where the turbulence model resolves the boundary layer. It depends on the height of the first-layer cell from the wall boundaries. Each turbulence model operates in a specific range of Y⁺ to accurately resolve the boundary layer. For the RNG model, the wall Y⁺ value should range between 30 and 300 [45]. This range was imperatively used to determine the corresponding first-layer height (Y_P) for the model using methods in Equation (11).

$${\mathrm{Y}}^{+}=\frac{{\mathrm{u}}_{\tau \mathrm{a}\mathrm{i}\mathrm{r}}{\mathrm{Y}}_{\mathrm{P}}}{{\mathrm{V}}_{\mathrm{a}\mathrm{i}\mathrm{r}}}; {\mathrm{u}}_{\tau \mathrm{a}\mathrm{i}\mathrm{r}}=\sqrt{\frac{{\tau }_{\mathrm{w}\mathrm{a}\mathrm{l}\mathrm{l}}}{\mathsf{\rho }}}, \mathrm{R}\mathrm{e}=\frac{\mathsf{\rho }{\mathrm{U}}_{\mathrm{\infty }}\mathrm{L}}{{\mathrm{V}}_{\mathrm{a}\mathrm{i}\mathrm{r}}}, {\mathrm{C}}_{\mathrm{f}}=\frac{0.026}{\mathrm{R}{\mathrm{e}}^{\frac{1}{7}}}, {\tau }_{\mathrm{w}\mathrm{a}\mathrm{l}\mathrm{l}}=\frac{{\mathrm{c}}_{\mathrm{f}}\mathsf{\rho }{\mathrm{U}}_{\mathrm{\infty }}^2}{2},$$

(11)

where Re is the Reynolds number, $\mathsf{\rho }$ is the air density $\left(\mathrm{kg}{ \times \mathrm{m}}^{-3}\right)$, L is the reference length (m), ${\mathrm{C}}_{\mathrm{f}}$ is an empirical constant, ${\mathrm{U}}_{\mathrm{\infty }}$ is the free-stream velocity of air ($\mathrm{m}\times {\mathrm{s}}^{-1}$), and ${\tau }_{\mathrm{w}\mathrm{a}\mathrm{l}\mathrm{l}}$ is the wall shear force ($\mathrm{kg}{\mathrm{m}}^{-1}\times {\mathrm{s}}^{-2})$.

In the current study, the grid refinement approach was used to perform the grid-independent study for establishing the most optimal grid conditions. Wall functions and Y⁺ approaches were used to optimize model accuracy. For verification, statistical indices R² and RMSE were used as judgment criteria for grid optimal conditions efficient for designing large-scale greenhouse models. The grid condition that satisfies the 5% requirement or lower for the RMSE and the value of 0.95 or higher for R² were defined as the optimal conditions for model verification.

2.5. Ventilation Efficiency Theories

To evaluate the natural ventilation, the tracer gas decay method (TGD) was used considering CO₂ gas [46,47]. TGD is a qualitative method that uses a tracer gas introduced in the structure to estimate the ventilation rates. It is based on the analysis of the decay rate of the tracer gas inside the study structure. TGD considers the greenhouse internal airflow characteristics, ventilator configuration, and the external environment; therefore, it has been widely used to calculate local ventilation rates inside structures. For this study, the TGD approach was used to estimate the regional ventilation rates inside the greenhouse. The air exchange rates were calculated using the equation as follows:

$$A{E}_{TGD}=\frac{\ln \left(\frac{{\mathrm{C}}_0}{{\mathrm{C}}_{\mathrm{t}}}\right)}{t-t_0}$$

(12)

where C₀ is the concentration of the tracer gas (ppm) when time is 0 s. C_t is the concentration at time t (s), and AE_TGD is the air exchange rate (AER, min⁻¹).

The ventilation requirement was used to evaluate the performance of the natural ventilation of the greenhouse. The ventilation requirement of a structure can be described as the amount of the ventilation rate to maintain an optimal environment, such as temperature, relative humidity, CO₂, etc., inside the structure for proper crop growth. In this study, the ventilation requirement for temperature control was computed by the solar radiation expected on-site (Equation (8)). The parameters used for the evaluation are shown in

Table 2. Ventilation requirement evaluation parameters and values following expected and designed solar radiation range.

Parameter	Symbol	Value
Volumetric specific heat of air (W·m−3·°C−1)	CV	0.3
Correction rate of heated area	α	1.2
Solar transmissivity of greenhouse cover	τ	0.7
Evapotranspiration coefficient (Tomato)	f	0.5
Overall heat transfer coefficient (Wm2·min−1·°C−1)	K	0.08
Greenhouse flow area (m2)	Af	19,200
Greenhouse cover area (m2)	Ac	24,946
Greenhouse volume (m3)	V	124,800

.

$$\mathrm{V}\mathrm{R}=\frac{1}{{\mathrm{C}}_{\mathrm{V}}}\left\{\frac{\mathsf{\alpha }\tau \mathrm{S}\left(1-\mathrm{f}\right){\mathrm{A}}_{\mathrm{f}}}{\nabla \mathrm{T}}-\mathrm{k}{\mathrm{A}}_{\mathrm{c}}\right\}\ast \frac{1}{\mathrm{V}}$$

(13)

2.6. Study Procedure

In this study, we analyzed the fluid dynamics and ventilation efficiency of the naturally ventilated Venlo greenhouse during the design phase of a large-scale facility. Based on meteorological data, we set the boundary conditions for the simulation and verified the model using the GIT and Y⁺ theory. Using the verified model, we conducted a fluid dynamics analysis of the airflow pattern into the naturally ventilated greenhouse and analyzed the regional ventilation efficiency using the TGD theory. Twenty years of weather data for the target area was analyzed, and average summer data were used to design the CFD model boundary conditions, vertical wind profile, and turbulence intensity. The most frequent wind speed range at both weather stations was 0.5–3.3 m·s⁻¹. The annual summer wind speed was then determined as 3.08 m·s⁻¹ and was used as the speed at the reference height velocity.

2.6.1. Greenhouse CFD Modeling and Validation

Field monitoring was conducted in the standard Venlo greenhouse, designed for 0.1 ha, 8 span, and 19 tomato crop rows following the structural design standard of 1999 from the Ministry of Food, Agricultural, Forestry and Fisheries of Korea. The standard greenhouse horizontal dimensions (32 m by 32 m) were upscaled by a factor of 3.5 in width and 5 in length to form a 30-span, 2 ha large-scale model as a target greenhouse in this study. The CFD model was designed by the previously verified CFD model by comparing the field-monitored data with the temperature deviation of 1.12 ± 0.64 °C which showed good agreement qualitatively and quantitatively [26]. To have more than one grid cell across the plant canopy rows during model discretization, 3 plant rows in the standard greenhouse were simplified into 1 row in the model greenhouse. Natural ventilation operation conditions were considered, with only roof-top ventilation. The roof vent cover opening angle was 60, to represent the average operation conditions during the summer season. The modeling process was implemented in the Ansys Space claim (Ver. 2022 R2, USA) software to make the 3D CAD model for the 2 ha greenhouse.

Figure 2 shows a schematic representation of the small-scale standard greenhouse and the upscaled greenhouse model with a 2 ha area. For simulation purposes, the greenhouse external environment was developed following the methods described in the Materials and Methods, Section 2.3. The overall computation domain was cylindrical with a radius of 105 m and a height of 70 m.

To apply a natural ventilation model, consideration of the external environment of the greenhouse is necessary. A CFD model takes into account that both the greenhouse and its exterior have a very wide domain resulting in the number of grids, which crucially affects the computation time and accuracy. A CFD model was used to assess the grid independence and optimization of the grid size. To make the CFD model, the model was sectioned along the center of the greenhouse across the crop canopy in the vertical plane. The cross-section model domain was discretized into computational grids using Ansys Meshing software. For reference data generation, the finest grid of the 0.2 m resolution model was developed. For the GIT study, several CFD models with coarser grid resolutions of 0.3 m, 0.4 m, 0.5 m, 0.6 m, 0.7 m, 0.8 m, 1.0 m, 2.0 m, and 3.0 m were considered for finding the optimal grid determination. The CFD models were simulated in Ansys Fluent software where the model boundary conditions and crop canopy models were implemented. The CFD model steady-state regime and simulations were run until the preset convergence criteria were met. The 0.2 m grid boundary conditions are illustrated in Figure 3, and the Ansys Fluent CFD settings, operating conditions, and models used are shown in

Table 3. Boundary conditions for 3D CFD model of the large-scale Venlo greenhouse with 2 ha area.

Variable	CFD Solver Setting
Solver	Pressure-based solver
Numerical algorithm	SIMPLE algorithm
Discretization	Second order
Operating pressure	101,325 Pa
Gravity	9.81 (m·s−2)
Air density	Variable (ideal gas)
Air viscosity	1.79 (kg·m−1·s−1)
Turbulence	RNG κ-ε
Empirical wall function	Enhanced wall functions
External wind velocity	3.1 (m·s−2)
Crop canopy	Porous media (0.8 m width, 3.5 m height)

. In the three-dimensional CFD model, the outer boundary region was divided into 16 parts to perform simulations considering different wind directions and speeds. The model was designed so that the windward side, where the air comes in, and the leeward side, where the air goes out, could be set using a single model. To increase the stability of the flow in the upper part of the greenhouse in the model, it is necessary to secure sufficient height. At the same time, by setting the ceiling as symmetrical, it was more flexible and able to maintain linearity than setting it as a wall [6,8]. The Venlo-type greenhouse has a floor area of 2 ha, and tomatoes were assumed as a porous media and reflected in the design (Figure 4). The model was initially simulated in the steady state until convergence and then switched to the transient regime for aerodynamic approaches with TGD. The transient simulation flow time was 10 min with 1 s time steps, including carbon-dioxide trace gas. The CFD model was computed on a computing system with an i7 CPU, 32 GB RAM, and a 1600 Ti GPU. After initializing in the steady-state condition, transient modeling was conducted for computing the ventilation efficiency with the tracer gas concentrations. The simulations were conducted for a total of 180 s, with a time step of 1 s, resulting in a total of 12~60 h to complete the CFD model according to the grid resolution from 0.2 to 3.0 m.

Model verification was performed using the GIT and Y⁺ methods. Based on the reviewed literature about greenhouse CFD modeling, a grid resolution of 0.2 m and the RNG k-ε turbulence model have been reported as optimal for designing the greenhouse CFD model through validation with field-monitored data. In this study, a 0.2 m grid resolution and the RNG k-ε turbulence model were used to design the CFD model for reference data that were considered as the baseline for model verification judgment.

The GIT was used to determine the minimum grid resolution required to obtain a converged and accurate solution. This was to ensure that the numerical solution is accurate and reliable, while minimizing computational costs. Wall Y⁺ was used to resolve the near-wall boundary layer resulting from the high solution gradient at the model walls. Model verification was performed by comparing air velocity results from the reference data which were validated by previous research. The model was considered verified when the following conditions were satisfied: R² ≥ 0.95, RMSE ≤ 5%, and 30 ≤ Y⁺ ≤ 300. The procedure used in model verification is shown in Figure 5. It involved 5 steps described below:

• Identifying a recommended grid size. The 0.2 m grid resolution and RNG k-ε turbulence model were adapted from the literature and used in the 3D CFD model to generate a reference set of data.
• Selection of grid resolution range for GIT study. With reference to the reviewed literature, a range between 0.3 m and 3 m for the resolution was selected for optimal resolution determination.
• GIT study. The average air velocity data for all the coarse grid resolutions were compared to the reference data. The grid resolution above, which the results started to diverge from the reference data, was selected as the optimal grid and tested in the 3D CFD model.
• Y⁺ study. The selected optimal grid was used to check the Y⁺ test with the averaged air velocity extracted at the study heights in the CFD model.
• Verification. The 3D results were statistically compared to the reference data by evaluating the R² and RMSE and checking the Y⁺. Near-wall mesh refinement using the first-layer height of 0.04 m was used to iteratively adjust the wall Y⁺ value into the target range. When the verification procedure was all satisfied, the natural ventilation, microenvironment distribution, and flow pattern were analyzed to assess the performance of the designed model.

Figure 5. Model verification study procedure.

Figure 5. Model verification study procedure.

2.6.2. Aerodynamic Analysis of CFD Greenhouse Model

The flow of air entering the naturally ventilated greenhouse was analyzed aerodynamically. Using the verified CFD model, the flow field was qualitatively and quantitatively analyzed according to horizontal and vertical cross sections. To analyze the airflow pattern and the local ventilation efficiency, the greenhouse was subdivided into 20 analysis regions (5 by 4), and each analysis region was identified based on their row and column positions in the greenhouse (Figure 6). The parameters were also analyzed according to height, at 0.7 m, 2.2 m, and 3.5 m heights within the crop canopy from the ground. In the model accuracy verification stage, the simulated data were extracted from 1000 points within the model at the three study heights from the center of the greenhouse for the reference data in the CFD verification model (Figure 7). The velocity distribution in the CFD model was also analyzed at 30 m and 60 m from the center of the greenhouse to both sides along the length of the greenhouse at all 3 analysis heights.

3. Results and Discussion

3.1. Greenhouse CFD Modeling and Validation

The CFD model was verified by the GIT and Y⁺ values. For the grid-independent test, we set the results from the CFD model calculation, which used the 0.2 m grid resolution and the same model design conditions as the validated model through previous research, as the reference model. For the reference model, the average air velocity in both the greenhouse and plant domains was 0.654 m·s⁻¹ (Figure 8). This value is close to 0.58 m·s⁻¹ as reported by [30] in a 1 ha greenhouse with a tomato canopy. Generally, average greenhouse air velocities have been reported to range between 0 and 1.5 m·s⁻¹ for the natural ventilation regime. The air velocity results based on larger grid resolutions between 0.3 m and 0.8 m showed no significant difference compared to the 0.654 m·s⁻¹ obtained from the simulation with the reference grid size (Figure 8). For grid sizes larger than 0.8 m, the average air velocity diverged from the reference value. The 0.8 m grid size was hence selected as the optimal grid resolution in the CFD model.

The maximum initial wall Y⁺ value was 1800 on the windward side of the greenhouse cover material (Figure 9). The average wall Y⁺ value was determined to be 600, which was far outside the target range. After evaluation, it was determined that a first-layer height of 0.04 m would be suitable to reduce the Y⁺ value within the required range for an accurate resolution of the boundary layer. This adjustment resulted in a maximum wall Y⁺ value of 270 in the CFD model, which satisfied the requirement for the RNG k-epsilon turbulence model with enhanced wall functions. To refine the near-wall zone, all the cells were refined using four inflation layers, with the first layer having a height of 0.04 m and a growth rate of 1.2 toward the external domain. This refinement only increased the total number of cells by 120,000. As a result, the overall 3D computational domain consisted of a grid with 5.5 million cells.

After a statistical comparison of the air velocity results between the designed CFD model with the reference model, the R² values were 0.984, 0.985, and 0.934 while the RMSE values were 5.61, 3.96, and 2.20 at the 0.7 m, 2.2 m, and 3.5 m heights, respectively (

Table 4. Analysis height-based statistical indices from comparing 0.8 m grid size average air velocity results to the reference data set generated from 0.2 m grid size.

Height	R-Square	RMSE
0.7	0.984	5.61
2.2	0.985	3.96
3.5	0.934	2.20
Average	0.968	3.92

). The highest R² value (0.984) and highest RMSE value (5.61) were observed at the 0.7 m height. The lowest R² and RMSE values (0.93 and 2.2) were observed at the 3.5 m height (Figure 10). This could be due to the effect of the direct influence of the external wind to the local airflow at 3.5 m in the greenhouse through the roof window vents (Figure 11). The average R² and RMSE for all three analysis heights were 0.968 and 3.923%, respectively, which both meet the verification conditions. Although there was some deviation in areas with low velocity, in areas with velocity above 0.5 m/s, which form the main flow, the distribution and pattern were very similar to those of the reference model.

3.2. Aerodynamic Analysis of CFD Greenhouse Model

3.2.1. Internal Airflow Pattern

The analysis of the greenhouse internal and external airflow patterns revealed that only a small fraction of the external air streams entering the greenhouse through the roof windows reached the interior. Due to the structural characteristics of the Venlo greenhouse, there are only windows on the roof. Thus, the external air entering from the side walls flows upwards following the structure of the greenhouse. This creates a negative pressure which causes the air to be drawn out from the windward side, and the air is then introduced through the inlet on the leeward side (Figure 12). When analyzing the internal flow distribution in the computed CFD results, the air entering near the roof window on the windward side tends to be exhausted almost immediately. This is believed to be due to the backflow formed inside. Elsewhere, due to the negative pressure, air is entering in the opposite direction from the leeward side. The inflowing air is concentrated and flows toward the center of the greenhouse. As a result, a large stagnation zone is formed on the leeward side (Figure 13). Most of the streams that entered the greenhouse were observed in the central region, flowing toward the leeward side where the pressure was lower. This flow distribution can be attributed to the insufficient direct entry of air streams from the outside into the greenhouse.

The air velocity profiles from the windward side of the greenhouse were analyzed vertically at the center and at locations 30 m and 60 m away from the center in both the left and right directions, illustrated in Figure 14. The air velocity profile at 30 m from the center, both to the right and left sides, closely resembled that at the center, with the highest velocity occurring at a height of 0.7 m and the lowest velocity at a height of 3.5 m. However, at 60 m from the center, significant variations in the velocity profiles were observed, with the velocities sharply dropping to 0 m/s at certain points. This result shows that the distribution of the wind velocity inside a 2 ha greenhouse varies with the horizontal distance from the center. Similar velocity distribution profiles and trends have been reported in several previous studies.

3.2.2. Regional Ventilation Efficiency

There are limitations to quantitatively evaluating the ventilation structure using only the airflow pattern and velocity distribution. To quantitatively analyze the uniformity of local ventilation efficiency, the TGD method using a CFD model was applied. The greenhouse CFD model was divided into a total of 60 detailed regions by height and locations to analyze the overall and regional ventilation efficiency. The ventilation efficiency was assessed by computing the local distributions of tracer gas concentrations within the greenhouse. The time-dependent variation in tracer gas concentration at a height of 2.2 m within the greenhouse and plant canopy, at flow times of 2, 4, 6, and 10 min, is illustrated in Figure 14. After filling the tracer gas uniformly inside the greenhouse CFD model, the concentration of the gas inside becomes diluted over time from the beginning of the ventilation. The change in concentration of the tracer gas over 10 min was applied to the TGD method to analyze the regional ventilation efficiency. As analyzed in the airflow pattern, due to the negative pressure formed outside, air backflows into the greenhouse from the leeward side, so the change in concentration shows the influence of this backflow.

As a result of analyzing the ventilation efficiency by height, it was observed that the efficiency increased with height (Figure 15). The lowest average air exchange rate (AER) was found to be 0.2 min⁻¹ at a height of 0.7 m, with the lowest standard deviation. Interestingly, despite the highest air velocities being observed at a height of 0.7 m, the ventilation efficiency was the lowest at this height. Conversely, the reverse trend was observed at a height of 3.5 m. These results suggest that the air at the bottom of the greenhouse had a higher mean age of air (MAA) due to the induced internal pressure distribution resulting from the crop canopy effect and the flow pattern at the roof window. The distribution of ventilation efficiency among the analysis regions revealed three distinct clusters, as indicated by the dotted lines in Figure 15. The regions that were furthest from the windward side exhibited the lowest air exchange rates, with an average of 0.2 min⁻¹. Generally, the ventilation efficiency exhibited a decreasing trend as the distance from the windward side increased.

Figure 16 represents the regional ventilation efficiency according to the height from the floor. The ventilation efficiency displayed an increasing trend with an increase in height. This is due to the flow of air entering the greenhouse from the roof and moving downwards resulting in the concentration of air upwards. As the air travels down, its movement becomes more stagnant, or it fails to reach certain areas due to the influence of crops. In zones such as R11–12 and R51–52, where the air flows from the center toward the sides and exhausts, the ventilation efficiency showed relatively superior results. In the leeward regions corresponding to R14, 24, 34, 44, and 54, there was almost no influx of fresh air, leading to the stagnant air zones. This can be attributed to the drag force effect of the crops on the airflow and the porous nature of the crop canopy implemented in the model.

The CFD model of the 2 ha Venlo greenhouse under tomato cultivation, with an external wind speed of approximately 3 m/s, exhibited internal and external airflow patterns resulting in air velocities ranging from 0 to 0.7 m/s and an average ventilation efficiency of 0.51 AER (min⁻¹) with a standard deviation of 0.25 considering the crop canopy. The overall ventilation efficiency of the greenhouse was calculated to be 0.51 AER (min⁻¹) which is very reasonable compared with previous studies on a plant-less 8-span (0.1 ha) Venlo greenhouse with only roof windows open, ventilation efficiencies of 0.28, and a 0.58 AER (min⁻¹) during the external wind speeds of 2.5 m·s⁻¹ and 5.5 m·s⁻¹, respectively [7].

When considering the summer season, the required ventilation rate can be determined based on the solar radiation in the reclaimed land and the set temperature difference between the inside and outside. Referring to a previous study [7], when the temperature difference between the inside and outside is set to 5 degrees, the required ventilation rates at solar radiation of 200, 400, 600, and 800 are suggested to be 0.3, 0.6, 0.9, and 1.2 AER (min⁻¹), respectively. The ventilation efficiency of 0.51 AER (min⁻¹), calculated through the simulation, appeared to be suitable only when the solar radiation is 300 W or less. These results describe a highly heterogeneous internal microclimate and an inefficient natural ventilation system for tomato cultivation. Under these conditions, an upgrade of the ventilation system is required to meet the appropriate growth conditions for tomato crops, such as a circulation fan, an internal section, and the vent design.

4. Conclusions

The authors of this study analyzed the airflow pattern and regional ventilation efficiency for designing a large-scale Venlo greenhouse in Saemangeum-reclaimed land in Korea. Understanding the air exchange mechanisms between the inside and outside of greenhouses is crucial for evaluating how the internal microclimate and ventilation efficiency interacted with the crop canopy. The design and validation of the CFD models considering the structural configuration, crop canopy, and external weather condition are critical steps before they can be used for a detailed analysis and making commercial decisions. In this study, it has been demonstrated that the model verification steps incorporating the improved grid independence test (GIT) and wall Y⁺ methods, an efficient procedure for design optimization, accuracy verification, and grid reduction in large-scale CFD models, can be developed and adopted. Even without extensive field data, these approaches yield accurate results that are statistically evaluated while significantly reducing the number of grid cells. This makes it feasible to model very large greenhouse complexes. Specifically, it has been shown that a grid resolution of 0.8 m within the greenhouse and plant domain is sufficient to obtain closely accurate results compared to using a finer grid size of 0.2 m. The accuracy verification results yielded low RMSE values (as low as 3.9%) and high R² values (0.968), resulting in a 38% reduction in the number of grid cells. This demonstrates the effectiveness of the proposed approach in achieving accurate results while reducing computational requirements. Through the CFD modeling analysis, it has been observed that the external wind directly influences the airflow direction inside the greenhouse, leading to negative pressure. This induces a stronger wind-wise air current above the crop canopy and a slow reverse flow inside the crop canopy, which aligns with findings from previous studies. When considering the air inflow parallel to the greenhouse, the average internal ventilation efficiency was shown to be 0.51 AER (min⁻¹). In the presence of 800 W·m⁻² of solar radiation, the results indicated that the required ventilation efficiency did not reach 1.2 AER (min⁻¹).

Future research is needed to find the optimal microclimate conditions for the cultivation of tomato crops according to the ventilation system. Utilizing the results of this study, it is essential to analyze the ventilation efficiency based on external weather conditions, including various wind directions and speeds. Furthermore, a microclimate analysis that encompasses sensible and latent heat emitted from crops and the ground is required. Based on this, we will derive improvement plans using ventilation structures, circulation fans, and windbreak facilities. The ultimate aim was to extrapolate the findings from the single unit to a greenhouse complex to explore the airflow pattern, ventilation rates, and fluid dynamics considerations caused by adjacent greenhouses when the design is expanded to the complex.