Modeling In-Vehicle VOCs Distribution from Cabin Interior Surfaces under Solar Radiation

: In-vehicle air pollution has become a public health priority worldwide, especially for volatile organic compounds (VOCs) emitted from the vehicle interiors. Although existing literature shows VOCs emission is temperature-dependent, the impact of solar radiation on VOCs distribution in enclosed cabin space is not well understood. Here we made an early e ﬀ ort to investigate the VOCs levels in vehicle microenvironments using numerical modeling. We evaluated the model performance using a number of turbulence and radiation model combinations to predict heat transfer coupled with natural convection, heat conduction and radiation with a laboratory airship. The Shear–Stress Transport (SST) k- ω model, Surface-to-surface (S2S) model and solar load model were employed to investigate the thermal environment of a closed automobile cabin under solar radiation in the summer. A VOCs emission model was employed to simulate the spatial distribution of VOCs. Our ﬁnding shows that solar radiation plays a critical role in determining the temperature distribution in the cabin, which can increase by 30 ◦ C for directly exposed cabin surfaces and 10 ◦ C for shaded ones, respectively. Ignoring the thermal radiation reduced the accuracy of temperature and airﬂow prediction. Due to the strong temperature dependence, the hotter interiors such as the dashboard and rear board released more VOCs per unit time and area. A VOC plume rose from the interior sources as a result of the thermal buoyancy ﬂow. A total of 19 mg of VOCs was released from the interiors within two simulated hours from 10:00 am to noon. The ﬁndings, such as modeled spatial distributions of VOCs, provide a key reference to automakers, who are paying increasing attention to cabin environment and the health of drivers and passengers.


Introduction
Over the past decades, China has been experiencing the world's fastest growth in vehicle population. As a result, commuters inevitably spend a substantial amount of time in vehicle cabins due to increased traffic congestions and vehicle population especially in major cities [1]. Epidemiological studies show that long-time exposure to air pollutants is associated with increased risks of morbidity and mortality [2][3][4], especially high volatile organic compounds (VOCs) concentrations emitted from cabin interiors [5][6][7] that could lead to respiratory irritation and cancer [8]. Developing advanced methods to identify in-cabin emission sources [9] has become a public health priority for consumers, car producers and government. Some attempts have been made to provide indoor air quality (IAQ) guidelines for passenger cars. For example, the Ministry of Environmental Protection of the People's Republic of China has promulgated national standard the HJ/T 400 "Determination of Volatile Organic Compounds and Carbonyl Compounds in Cabin of Vehicles" and GB/T 27,630 "Guideline for air

The Turbulence Models
Turbulence modeling is a critical process for the numerical investigation of thermal environments. There have been some available studies to evaluate the performance of various turbulence models. Zhai et al. [42] compared eight turbulence models for predicting airflow and turbulence in enclosed environments and found that the Re-Normalization Group (RNG) k-ε model performed best among the Reynolds-averaged Navier-Stokes (RANS) models. Hussain et al. [43] used six RANS turbulence models to simulate the thermal environment in an atrium, and found the Shear-Stress Transport (SST) k-ω model provided comparatively better results. Li et al. [44] presented a numerical evaluation of the eddy viscosity turbulence models in terms of CFD modeling of convection-radiation coupled heat transfer in the indoor environment, and demonstrated a great performance of k-ω group models. In this study, five RANS turbulence models (including the standard k-ε model, the RNG k-ε model, the realizable k-ε model, the standard k-ω model, the SST k-ω model) and the Large Eddy Simulation (LES) model were selected to evaluate the prediction of the airflow and temperature distributions in the vehicle cabin.

The Radiation Models
Thermal radiation makes the temperature distribution more uniform in an enclosed space by transferring thermal energy from a hot surface to a cold one. Generally, thermal radiation accounts for 30-70% of the total heat transfer rate [45]. In this study, radiation heat transfer was taken into consideration for modeling the cabin thermal environment.

The Turbulence Models
Turbulence modeling is a critical process for the numerical investigation of thermal environments. There have been some available studies to evaluate the performance of various turbulence models. Zhai et al. [42] compared eight turbulence models for predicting airflow and turbulence in enclosed environments and found that the Re-Normalization Group (RNG) k-ε model performed best among the Reynolds-averaged Navier-Stokes (RANS) models. Hussain et al. [43] used six RANS turbulence models to simulate the thermal environment in an atrium, and found the Shear-Stress Transport (SST) k-ω model provided comparatively better results. Li et al. [44] presented a numerical evaluation of the eddy viscosity turbulence models in terms of CFD modeling of convection-radiation coupled heat transfer in the indoor environment, and demonstrated a great performance of k-ω group models. In this study, five RANS turbulence models (including the standard k-ε model, the RNG k-ε model, the realizable k-ε model, the standard k-ω model, the SST k-ω model) and the Large Eddy Simulation (LES) model were selected to evaluate the prediction of the airflow and temperature distributions in the vehicle cabin.

The Radiation Models
Thermal radiation makes the temperature distribution more uniform in an enclosed space by transferring thermal energy from a hot surface to a cold one. Generally, thermal radiation accounts for 30-70% of the total heat transfer rate [45]. In this study, radiation heat transfer was taken into consideration for modeling the cabin thermal environment.
The P-1 radiation model is the simplest case of the P-N model, which is based on the expansion of the radiation intensity into an orthogonal series of spherical harmonics functions. The directional Sustainability 2020, 12, 5526 4 of 19 dependence in radiative transfer equation (RTE) is integrated out, resulting in a diffusion equation for incident radiation. Equation (1) is obtained for the radiation flux q r , only considering scattering and absorption when modeling gray radiation.
An advection-diffusion equation is solved to determine the local radiation intensity G in the P-1 model.
where q r is the radiation heat flux, a is the absorption coefficient, σ s is the scattering coefficient, G is the incident radiation, C is the linear-anisotropic phase function coefficient, σ is the Stefan-Boltzmann constant. The Surface to Surface (S2S) radiation model is applicable for modeling radiation in situations where there are no participating media. All surfaces involved in radiation are assumed to be gray and diffuse, ignoring absorption, emission and scattering and preserving only "surface-to-surface" radiation. The energy flux leaving a given surface is composed of directly emitted and reflected energy. The reflected energy flux is dependent on the incident energy flux from the surroundings, which then can be expressed in terms of the energy flux leaving all other surfaces. The energy leaving from the kth adiabatic surface can be expressed as Equation (3).
where q out,k is the energy flux leaving the surface, ε k is the emissivity, ρ k is the reflectivity of surface k, q in,k is the energy flux incident on the surface from the surroundings. The Discrete Ordinates (DO) radiation model is regarded as the most comprehensive radiation model, which solves the RTE for a discrete number of finite solid angles, as shown in Equation (4). Each associated with a vector direction → s is fixed in the global Cartesian system (x, y, z). Accuracy can be increased by using a better discretization, while it may be CPU-intensive with many ordinates.
where → r is the position vector, → s is the direction vector, → s is the scattering direction vector, α is the absorption coefficient, n is the refractive index, I is the radiation intensity, which depends on the position → r and direction → s , T is the local temperature, Φ is the phase function, Ω is the solid angle. For the Discrete Transfer Radiation Model (DTRM), the main assumption is that radiation leaving a surface element within a specified range of solid angles can be approximated by a single ray. The energy source in the fluid due to radiation is computed by summing the change in intensity dI along the path of each ray ds that is traced through the fluid control volume.

The Solar Load Model (SLM)
Thermal energy due to incident solar rays is a very common but important phenomenon. In the SLM, the solar beam direction and irradiation were calculated according to the provided geographical place and the given time based on the solar load model's ray tracing algorithm. A two-band spectral model was used for direct solar illumination and accounted for separate material properties in the visible and infrared bands. A single-band hemispherical-averaged spectral model was used for diffuse Sustainability 2020, 12, 5526 5 of 19 radiation. The ground was considered to be dry bare land with a reflectivity of 0.2 [40]. Solar scattering was set to the default value of 1.

Contaminants Emission Model
For the VOCs emission of building materials, the quasi-steady-state ER can be described by the following equation [46]: where E is the emission rate factor, t is the emission time, δ is the material thickness, C 0 is the initial emittable concentration inside the material, D m is the diffusion coefficient. The correlation between C 0 and T can be described by the following equation [47]. The correlation between D m and T can be described by the following equation [48]. Equation (9) is obtained by substituting Equations (7) and (8) into Equation (6) and then taking the logarithm on both sides.
where T is the temperature in K. C 1 , C 2 , D 1 and D 2 are all constants determined only by the physical and chemical properties of pollutants To quantify the pollutants released by materials in a certain period of time, Equation (10) is included.
where M is the total quality of the released contaminants, t s is the start time, t e is the end time, E(t) is the emission rate, A is the surface area.

Model Validation and Discussion
The temperature and flow field in the airship is a result of the interplay among multiple factors including solar radiation, earth reflection, infrared radiation, external forced convection and internal natural convection, which is similar to the thermal environment of a parked car. Therefore, the experimental results of Li et al. [49] were selected for model validation. In their study, the transient thermal behaviors of an airship under different solar radiation were revealed. An airship model with a spherical tank type was built in a closed laboratory. The body was covered with 0.1 mm polyimide film, and shaped by thin metal sheets. Solar irradiation was supplied by a TRM-PD solar simulator and measured by a XLP12-1S-H2 heat flow meter. Eighteen T-type thermocouples were arranged in 18 different locations to obtain the hull and inner gas temperatures. Eight points are located on plane 1 including from point 1 to point 6 and point 8. For plane 2, there are also eight points, from point 9 to point 14 and point 16. Point 18 and point 19 are distributed on plane 3. Plane 1 and plane 2 are symmetrical about plane 3 with 150 mm axial distance. In Figure 2, a full-size computational domain was modeled with the same physical dimensions as the experimental airship. The computational domain was discretized using the structured mesh with five refined inflation layers applied close to the solid surfaces. The grid independence was achieved at 0.2 million mesh elements, as shown in Table 1. To evaluate the CFD approach to model convection-radiation coupled heat transfer in the Sustainability 2020, 12, 5526 6 of 19 airship, the hull and inner gas temperatures were compared with experimental data in terms of the accuracy and computational cost. Six turbulence models and four radiation models were included. In each configuration, one turbulence model was combined with or without a radiation model. Finally, thirty CFD cases were obtained totally.
Sustainability 2020, 12, x FOR PEER REVIEW 6 of 20 was modeled with the same physical dimensions as the experimental airship. The computational domain was discretized using the structured mesh with five refined inflation layers applied close to the solid surfaces. The grid independence was achieved at 0.2 million mesh elements, as shown in Table 1. To evaluate the CFD approach to model convection-radiation coupled heat transfer in the airship, the hull and inner gas temperatures were compared with experimental data in terms of the accuracy and computational cost. Six turbulence models and four radiation models were included.
In each configuration, one turbulence model was combined with or without a radiation model. Finally, thirty CFD cases were obtained totally. Figure 2. The airship computational domain according to [49].
The hull material is polyimide film with thermal conductivity λ = 0.32w/( • ). The total solar absorptivity of the external surface is 0.45 with a 0.81 absorptivity in the infrared spectrum. The coupled solver with pseudo-transient relaxation was applied for the solution of the momentum, energy, and turbulence equations. The total calculation adopted 600 s according to the experimental sampling time. The results of 30 designed test trips with different combinations of radiation and turbulence models are obtained. Except for the air temperature measured at point 17 and 18, the rest are measured at the airship surface. Since the solar simulator is installed above the airship and keeps both axes parallel, the temperature distributions on plane 1 and 2 are basically the same. The temperature of six locations including point 1, point 3, point 6, point 17 and point 18 are shown in Figure 2. The airship computational domain according to [49]. The hull material is polyimide film with thermal conductivity λ = 0.32w/(m·K). The total solar absorptivity of the external surface is 0.45 with a 0.81 absorptivity in the infrared spectrum. The coupled solver with pseudo-transient relaxation was applied for the solution of the momentum, energy, and turbulence equations. The total calculation adopted 600 s according to the experimental sampling time.
The results of 30 designed test trips with different combinations of radiation and turbulence models are obtained. Except for the air temperature measured at point 17 and 18, the rest are measured at the airship surface. Since the solar simulator is installed above the airship and keeps both axes parallel, the temperature distributions on plane 1 and 2 are basically the same. The temperature of six locations including point 1, point 3, point 6, point 17 and point 18 are shown in Figure 3. The measured data at each point is selected as the benchmark. The temperature of the measuring point closer to the light source is higher. Among them, the temperature at point 6 is the highest, reaching more than 90 • C, from where the temperature drops along with the body surface. Due to the shelter of the airship, point 1 is not directly exposed to the sunlight, resulting in a similar temperature as that of the surroundings with 16.7 • C. Sustainability 2020, 12, x FOR PEER REVIEW 7 of 20 Figure 3. The measured data at each point is selected as the benchmark. The temperature of the measuring point closer to the light source is higher. Among them, the temperature at point 6 is the highest, reaching more than 90 ℃, from where the temperature drops along with the body surface. Due to the shelter of the airship, point 1 is not directly exposed to the sunlight, resulting in a similar temperature as that of the surroundings with 16.7 ℃. To evaluate the performance of model combinations, a method called the total temperature error is proposed to assess the accuracy of numerical simulations. The total temperature error is calculated according to Equation 11. In addition, the computational cost is obtained according to the CPU calculation time when using a computer with an AMD Ryzen 7 1700X, 3.40 GHz CPU, 32 GB of RAM and 8 compute nodes.
where is the simulated temperature of the point, is the measured temperature of the point. In the experiment, 18 points were measured totally.
In Figure 4, comparing the temperature errors with and without radiation, the predicted air and solid surface temperature profiles agree better with the experimental results when the effects of thermal radiation are accounted for in the numerical investigation. On the contrary, the total error when ignoring the surface radiation may be twice that when including the surface radiation. When the turbulence model is fixed and combined with different radiation models, it is found that the S2S and DO models perform best, while the DTRM model has the lowest accuracy. Among all turbulence models, the standard k-ω and SST k-ω models have very close accuracies to predict the temperature To evaluate the performance of model combinations, a method called the total temperature error is proposed to assess the accuracy of numerical simulations. The total temperature error is calculated according to Equation (11). In addition, the computational cost is obtained according to the CPU calculation time when using a computer with an AMD Ryzen 7 1700X, 3.40 GHz CPU, 32 GB of RAM and 8 compute nodes.
where T isimulated is the simulated temperature of the i th point, T imeasured is the measured temperature of the i th point. In the experiment, 18 points were measured totally. In Figure 4, comparing the temperature errors with and without radiation, the predicted air and solid surface temperature profiles agree better with the experimental results when the effects of thermal radiation are accounted for in the numerical investigation. On the contrary, the total error when ignoring the surface radiation may be twice that when including the surface radiation. When the turbulence model is fixed and combined with different radiation models, it is found that the S2S and DO models perform best, while the DTRM model has the lowest accuracy. Among all turbulence models, the standard k-ω and SST k-ω models have very close accuracies to predict the temperature profiles. The SST k-ω model has a clear advantage in predicting accuracy compared to the LES model; it works the best for the high Rayleigh number buoyancy-driven flow [50]. Figure 4 reveals an M-shape trend of the calculation time according to the order of the combined model. The radiation model has a more significant impact on computation time than the turbulence model. The S2S and P1 radiation models require much less computation time than the DO model, although some small disparities exist but still are comparable. The DTRM is not compatible with parallel processing, so it will consume Sustainability 2020, 12, 5526 8 of 19 more time to model radiative heat transfer. The above model validation and comparison suggest that the SST k-ω and S2S model are the best choices to predict the in-vehicle thermal environment under solar radiation. Some existing studies [51] also show that the SST k-ω model performs slightly better than the RNG k-ε model when simulating convection-radiation coupled heat transfer.
profiles. The SST k-ω model has a clear advantage in predicting accuracy compared to the LES model; it works the best for the high Rayleigh number buoyancy-driven flow [50]. Figure 4 reveals an Mshape trend of the calculation time according to the order of the combined model. The radiation model has a more significant impact on computation time than the turbulence model. The S2S and P1 radiation models require much less computation time than the DO model, although some small disparities exist but still are comparable. The DTRM is not compatible with parallel processing, so it will consume more time to model radiative heat transfer. The above model validation and comparison suggest that the SST k-ω and S2S model are the best choices to predict the in-vehicle thermal environment under solar radiation. Some existing studies [51] also show that the SST k-ω model performs slightly better than the RNG k-ε model when simulating convection-radiation coupled heat transfer. To validate the theoretical correlation between the ER and surface temperature, some experimental results in the previous literature were used. Table 2 displays the measured emission rates of five pollutants from the car mat under three varied temperatures and total volatile organic compounds (TVOC) from PBS-C at five different temperatures. The detailed experimental description can be found in reference [52] and [53].  To validate the theoretical correlation between the ER and surface temperature, some experimental results in the previous literature were used. Table 2 displays the measured emission rates of five pollutants from the car mat under three varied temperatures and total volatile organic compounds (TVOC) from PBS-C at five different temperatures. The detailed experimental description can be found in references [52,53]. The unit of ER for the volatile organic compounds (VOCs) and total volatile organic compounds (TVOC) is µg·m −2 ·h −1 and mg·m −2 ·h −1 , respectively.
The linear curve fittings are illustrated in Figure 5. All R 2 are greater than 0.95, indicating a satisfactory correlation between the VOCs emission rate and temperature proposed in Equation (9).
In the follow-up study, the little impact of existing in-air VOCs on the emission rate is ignored during the emission period [54]. Moreover, the materials maintain their original appearance and properties without any bake out treatment, which indicates abundant VOCs to be volatilized in a quasi-steady state. Besides, the effect of relative humidity on the emission factor is ignored because of the constant relative humidity in that environment [35].  The linear curve fittings are illustrated in Figure 5. All are greater than 0.95, indicating a satisfactory correlation between the VOCs emission rate and temperature proposed in Equation 9. In the follow-up study, the little impact of existing in-air VOCs on the emission rate is ignored during the emission period [54]. Moreover, the materials maintain their original appearance and properties without any bake out treatment, which indicates abundant VOCs to be volatilized in a quasi-steady state. Besides, the effect of relative humidity on the emission factor is ignored because of the constant relative humidity in that environment [35]. Figure 6 presents the vehicle model established referring to the original data of a hatchback, in which the engine and luggage compartment were ignored. In addition, some interior components including safety belts, electrical equipment and console details were neglected due to the fewer effects on airflow. The details of the dimension are shown in Table 3. The passenger compartment was designed with eight air conditioning inlets and two outlets. Four inlets are located on the center console. Among them, No.1 and No.2 are, respectively, located on the middle, facing the gap between the two bucket seats. No.3 and No.4 are located on the sides. No.5 and No.6 are, respectively, arranged in the feet space of the driver and passenger. No.7 and No.8 are located on the armrest box, facing the bench seat. The air outlets No.1 and No.2 are located on the rear board. No additional airflow inlet and outlet are included anymore under the well-sealed assumption. The vehicle computational domain was discretized using unstructured tetrahedron grids [55]. Inflation layers were used at the interfaces between the air volume and the solids. The case study adopted 3.8 million elements after converting the domain to the polyhedral meshes.  Figure 6 presents the vehicle model established referring to the original data of a hatchback, in which the engine and luggage compartment were ignored. In addition, some interior components including safety belts, electrical equipment and console details were neglected due to the fewer effects on airflow. The details of the dimension are shown in Table 3. The passenger compartment was designed with eight air conditioning inlets and two outlets. Four inlets are located on the center console. Among them, No.1 and No.2 are, respectively, located on the middle, facing the gap between the two bucket seats. No.3 and No.4 are located on the sides. No.5 and No.6 are, respectively, arranged in the feet space of the driver and passenger. No.7 and No.8 are located on the armrest box, facing the bench seat. The air outlets No.1 and No.2 are located on the rear board. No additional airflow inlet and outlet are included anymore under the well-sealed assumption. The vehicle computational domain was discretized using unstructured tetrahedron grids [55]. Inflation layers were used at the interfaces between the air volume and the solids. The case study adopted 3.8 million elements after converting the domain to the polyhedral meshes.   Figure 7 shows the thermal energy transfer for the cabin. External heat enters the cabin through three ways including heat conduction, heat convection and thermal radiation. When the vehicle is parked under the sunlight for a soaking period, some solar radiation enters the passenger compartment passing through the windows, some is reflected by the solid envelope, the rest is absorbed. Solar radiation leads to a considerable thermal load through heating the envelope and interiors. In addition, the cabin exchanges the heat with the external environment through the coupled convection and radiation. Due to the uniform temperature distribution, the airflow cycles are driven by buoyant force, creating the natural convection in the cabin. The scorching air is trapped inside the cabin due to the lack of openings, resulting in the greenhouse effect [56]. In this study, the simulation was based on the city of Hangzhou (118°21′-120°30′E, 29°11′-30°33′N), the capital of Zhejiang Province located along Southeast coast of China, characterized by long and hot summers. The ambient conditions were chosen on June 21 (summer solstice). The windshield orientation was to the south.  Figure 7 shows the thermal energy transfer for the cabin. External heat enters the cabin through three ways including heat conduction, heat convection and thermal radiation. When the vehicle is parked under the sunlight for a soaking period, some solar radiation enters the passenger compartment passing through the windows, some is reflected by the solid envelope, the rest is absorbed. Solar radiation leads to a considerable thermal load through heating the envelope and interiors. In addition, the cabin exchanges the heat with the external environment through the coupled convection and radiation. Due to the uniform temperature distribution, the airflow cycles are driven by buoyant force, creating the natural convection in the cabin. The scorching air is trapped inside the cabin due to the lack of openings, resulting in the greenhouse effect [56]. In this study, the simulation was based on the city of Hangzhou (118 • 21 -120 • 30 E, 29 • 11 -30 • 33 N), the capital of Zhejiang Province located along Southeast coast of China, characterized by long and hot summers. The ambient conditions were chosen on June 21 (summer solstice). The windshield orientation was to the south.

The Thermal Setup
The cabin enclosure is composed of body structure, door, floor and glass, exposed to outdoor climatic conditions with sun and wind directly; the mixed wall thermal boundary was uniformly

The Thermal Setup
The cabin enclosure is composed of body structure, door, floor and glass, exposed to outdoor climatic conditions with sun and wind directly; the mixed wall thermal boundary was uniformly applied to the shell. According to the local summer weather condition, the ambient temperature of 38 • C was set as free steam and external radiation temperature for the whole cabin body, except for the driver's foot space adjacent to the engine cooling water. The appropriate thermal resistance across the wall thickness was imposed according to the wall thickness and material properties. The details are presented in Table 4. Besides, the shell conduction approach was utilized to model conduction in the planar direction of steel body with good thermal conductivity. All solid surfaces were considered stationary walls with No-slip conditions. The convective heat transfer coefficient was calculated based on the empirical formula in Equation (12) [29]. The inlets and outlets of the HVAC system were treated as the wall when the vehicle kept ventilation off. h = 1.163 4 + 12v 0.5 (12) where v is the wind speed relative to the parked vehicle with 0.2 m/s, h = 10.89 W/m 2 K.

The Radiation Setup
The windshield, side window, and rear window were treated optically as semi-transparent walls; all other surfaces were considered opaque. All surfaces participated in radiation heat transfer. The emissivity of interior surfaces was assumed to be 0.95 [30], and 0.88 for the windows [57]. The optical properties of the cabin surfaces are listed in Table 5.

The Emission Model Setup
Due to the common HVAC operation during driving and brief natural ventilation when getting off, it is reasonable to expect a very low concentration of contaminants left in the cabin. Moreover, the existence of contaminants has no impact on the airflow and concentration dispersion. Under such reasonable assumptions, defining the user-defined source (UDS) equation is computationally less expensive compared to the multi-component Eulerian approach when modeling the gas transport [51]. Therefore, the unsteady variation, convection, diffusion, and generation in the domain were calculated using a UDS coupled with the user-defined functions (UDFs) based on the existing flow parameters.
where ρ is the density of air, S is a scalar representing the contaminant concentration, Γ is the molecular diffusivity of S.

Results and Discussion
Chemical mass balance results demonstrated that carpet and seats are the most important VOCs source inside a new vehicle [58], so the bucket, bench seats and floor were chosen as the conventional VOCs sources in the cabin. Besides, the contaminant emission behavior of the dashboard and rear board were additionally investigated due to their prominent representation of high solar exposure.
Solar flux and average temperature variations and distributions with respect to soaking process are shown in Figures 8 and 9, respectively. The transmitted solar flux shows a decrease tendency on window No.3 and No.4, whereas with a sustained growth for the other envelope throughout the soaking period. Most direct solar irradiation (over 350 W) enters the cabin through the windshield and rear glass, falling on the dashboard and rear board. Only a small portion of the sun's rays (below 40 W) passes through the side windows. At noon, the sun moves directly above the vehicle with about 180 • solar incidence angle, causing a similar solar load on both side body. In Figure 8b, the surface temperature rises steadily with the increase of solar intensity. The highest average temperature with more than 60 • C occurs at the dashboard and rear board, while the floor has the lowest temperature with about 44 • C. The average temperature of bucket seat No.1 is slightly higher than that of No.2. This is due to the reason that the former receives more solar load from both the windshield and left glass (Figure 9a). The temperature on the two bucket seats tends to be the same at noon because of the comparable sun exposure (Figure 9b). The heat is difficult to transfer around the interior surface by conduction due to low thermal conductivity, exacerbating the thermal imbalance and local overheating [59].    Solar flux and average temperature variations and distributions with respect to soaking process are shown in Figure 8 and Figure 9, respectively. The transmitted solar flux shows a decrease tendency on window No.3 and No.4, whereas with a sustained growth for the other envelope throughout the soaking period. Most direct solar irradiation (over 350 W) enters the cabin through the windshield and rear glass, falling on the dashboard and rear board. Only a small portion of the sun's rays (below 40 W) passes through the side windows. At noon, the sun moves directly above the vehicle with about 180° solar incidence angle, causing a similar solar load on both side body. In Figure  8b, the surface temperature rises steadily with the increase of solar intensity. The highest average temperature with more than 60 ℃ occurs at the dashboard and rear board, while the floor has the lowest temperature with about 44 ℃. The average temperature of bucket seat No.1 is slightly higher than that of No.2. This is due to the reason that the former receives more solar load from both the windshield and left glass (Figure 9a). The temperature on the two bucket seats tends to be the same at noon because of the comparable sun exposure (Figure 9b). The heat is difficult to transfer around the interior surface by conduction due to low thermal conductivity, exacerbating the thermal imbalance and local overheating [59].  Figure 10 displays the temperature distributions of airflows in the driver plane. The air temperature near the hot interior surface is substantially higher than that which is far away from it; this is because of where the thermal boundary layers exist. In addition, the upper air owns a higher temperature because of the hot air rising and more heat sources. The hotter air gradually develops towards the floor, forming obvious temperature stratification phenomenon. The temperature in the driver's head position reaches 59 • C. At 10:00 am, the airflow crosses the bucket to the rear compartment, while it turns into a flow recirculation in the located temperature layer at 11:00 am and noon. Our analysis shows that the temperature distribution is primarily affected by solar radiation and airflow itself. Due to the heat exchange and direct solar heating, the temperature rises as much as 30 • C for direct exposed cabin surfaces and 10 • C for shaded ones, respectively, which supports the claim that solar radiation is a necessity in the thermal environment simulation.  Figure 10 displays the temperature distributions of airflows in the driver plane. The air temperature near the hot interior surface is substantially higher than that which is far away from it; this is because of where the thermal boundary layers exist. In addition, the upper air owns a higher temperature because of the hot air rising and more heat sources. The hotter air gradually develops towards the floor, forming obvious temperature stratification phenomenon. The temperature in the driver's head position reaches 59℃. At 10:00 am, the airflow crosses the bucket to the rear  Figure 11 shows the VOCs distributions in the whole cabin at 10:00 am and noon, respectively. And the driver plane was selected in the compartment to analyze the VOCs distributions under solar radiation in Figure 12. It is clearly noticed that the VOCs concentration at the hotter contaminant source is remarkably larger due to the strong dependence between the VOCs emission and the surface temperature. The dashboard and rear board are exposed to the strongest sunlight uniformly, therefore releasing more TVOC at the per unit area. The VOCs emission from the higher temperature area at the carpet increases approximately five-fold compared to the unexposed region. From the view of time, the concentration is about 3-4 times higher at noon than 10:00 am. Besides, the paths of concentration distribution and dispersion are significantly different in the four cases. A remarkable pollutant plume above the hotter surfaces is observed in Figure 12, which confirms that the near-wall thermal buoyancy flows are captured by the adopted turbulence and radiation models. The higher VOCs concentration distributes below the driver's knees, which is difficult to diffuse above the dashboard and rear board by natural convection (Figure 11a,b). Most TVOC is concentrated on the driver's head and above. As the concentration increases, a small amount of TVOC moves towards the floor under the driving of concentration difference. On the roof, the pollutants released from the dashboard and rear board tend to form a bridge of high concentration contaminants (Figure 11c,d). All pollutants gather near the surface sources and diffuse throughout the surroundings without ventilation. It demonstrates that thermal buoyancy and natural convection play an important role in dissipating pollutants from the contaminant source to the adjacent air in the enclosed environment.  Figure 10 displays the temperature distributions of airflows in the driver plane. The air temperature near the hot interior surface is substantially higher than that which is far away from it; this is because of where the thermal boundary layers exist. In addition, the upper air owns a higher temperature because of the hot air rising and more heat sources. The hotter air gradually develops towards the floor, forming obvious temperature stratification phenomenon. The temperature in the driver's head position reaches 59℃. At 10:00 am, the airflow crosses the bucket to the rear compartment, while it turns into a flow recirculation in the located temperature layer at 11:00 am and noon. Our analysis shows that the temperature distribution is primarily affected by solar radiation and airflow itself. Due to the heat exchange and direct solar heating, the temperature rises as much as 30℃ for direct exposed cabin surfaces and 10℃ for shaded ones, respectively, which supports the claim that solar radiation is a necessity in the thermal environment simulation. Figure 11. VOCs concentration emitted from carpet at 10:00 am (a) and noon (b); TVOC released from dashboard, rear board and seats at 10:00 am (c) and noon (d). Figure 11 shows the VOCs distributions in the whole cabin at 10:00 am and noon, respectively. And the driver plane was selected in the compartment to analyze the VOCs distributions under solar radiation in Figure 12. It is clearly noticed that the VOCs concentration at the hotter contaminant Figure 11. VOCs concentration emitted from carpet at 10:00 am (a) and noon (b); TVOC released from dashboard, rear board and seats at 10:00 am (c) and noon (d). Table 6 lists the amount of pollutants released by materials from 10:00 am to noon, which was calculated based on the Equation (10). There is a total of 35.08 µg VOCs emitted from carpet and 19 mg TVOC released from other interiors, which far exceeds the national standards of many countries. The seats become the largest source of pollutants due to the larger surface area. Moreover, as illustrated in Figure 10, above the seats are the places where pollutants are most likely to accumulate. It will pose a great threat to health if the level of pollutant exposure cannot be mitigated effectively by natural ventilation or HVAC systems when drivers and passengers re-enter the car. on the driver's head and above. As the concentration increases, a small amount of TVOC moves towards the floor under the driving of concentration difference. On the roof, the pollutants released from the dashboard and rear board tend to form a bridge of high concentration contaminants ( Figure  11c and 11d). All pollutants gather near the surface sources and diffuse throughout the surroundings without ventilation. It demonstrates that thermal buoyancy and natural convection play an important role in dissipating pollutants from the contaminant source to the adjacent air in the enclosed environment. Figure 12. The pollutant distributions on the driver plane from 10:00 am to noon.  Table 6 lists the amount of pollutants released by materials from 10:00 am to noon, which was calculated based on the Equation 10. There is a total of 35.08 μg VOCs emitted from carpet and 19 mg TVOC released from other interiors, which far exceeds the national standards of many countries. The seats become the largest source of pollutants due to the larger surface area. Moreover, as illustrated in Figure 10, above the seats are the places where pollutants are most likely to accumulate. It will pose a great threat to health if the level of pollutant exposure cannot be mitigated effectively by natural ventilation or HVAC systems when drivers and passengers re-enter the car.

Conclusions
The in-vehicle VOCs exposure is a public health concern worldwide. This study explored the thermal environment of an in-vehicle cabin under solar radiation and provided an early effort to quantify the spatial distribution of VOCs released from interior surfaces. In this study, the in-cabin thermal environment was simulated by considering in-cabin natural convection-conduction coupled with solar radiation. The performance of different combinations of turbulence and radiation models was validated with measured temperature data in a reduced-scale airplane cabin. Combining the SST k-ω turbulence model with a surface-to-surface radiation model (S2S) performed best in terms of computational accuracy among 30 different combinations, including turbulence models such as the standard k-ε, RNG k-ε, realizable k-ε, standard k-ω, SST k-ω and LES model and radiation models such as the P-1, S2S, DO, DTRM model. Our findings suggest that solar radiation plays a critical role in determining the temperature distribution in the cabin, which can increase as much as 30 • C for direct exposed cabin surfaces and 10 • C for shaded ones. The maximum average temperature was observed over 60 • C on the dashboard and rear board that were exposed the direct sunlight. The lowest temperature was found on the cabin floor at 44 • C. Such a high cabin temperature profile considerably lowers the thermal comfort and promotes VOCs emissions due to a strong temperature dependence. The dispersion of VOCs strongly depended on the local emission rate and airflows. The dashboard and rear board were shown to have a larger emission rate than that of the seats and floor because of the higher surface temperature. The VOC plume from the seats rose upward towards the ceiling, whereas the VOCs from the floor stayed below the seats. From the 2-h simulation period from 10:00 am to noon, there was a total of 35.08 µg VOCs and 19.02 mg TVOC accumulated throughout the cabin. With increasing attention to the cabin environment and the health of drivers and passengers, the findings, such as modeled spatial distributions of VOCs, provide automakers with an important design guide that could improve the current ventilation design for the summer time.