Efficient Simulation of Convective Ovens in Automotive Paintshops

The surface treatment is the process in an automotive factory that consumes most energy, water, and chemicals, and produces most waste and pollution. The ventilation in the paint booths and the curing ovens are using most of the energy. To accurately simulate the heat transfer in automotive curing ovens is essential to make it possible to optimize the processes with respect to quality, cost, and environmental impact. Significantly less prototypes then needs to be physically tested, a more uniform curing improves the product quality, and thermal deformation and stresses can be reduced. In Ref. [1], we studied heat transfer in a section of a convective curing oven, and optimized the position of objects in the oven with respect to energy consumption and curing time. In this paper, we extend and continue that work to full scale curing ovens. This complements earlier novel tools for simulation of the spray and sealing processes [2–5]. The modeling and simulation of the convective ovens typically used in the automotive paintshops to cure the different paint layers is challenging and includes multiple scales, turbulent air flows, thin boundary layers, large temperature gradients, and long curing times. A brute force conjugate heat transfer simulation of an oven resolving all time and length scales would be enormously time and resource consuming. Therefore, mathematical modeling is needed to obtain realistic simulations times.

We present a novel approach developed in a research project together with the Swedish automotive industry, which makes it possible to accurately simulate a curing oven in almost real-time. The goal is to successfully predict the time-dependent object temperature to decrease the number of physical tests that need to be carried out, especially during the production preparation phase. It also allows the oven operator to investigate possible alternative settings of the oven, e.g., flow rates and temperatures. In the approach, the individual nozzles in an oven are simulated to estimate the local nozzle Nusselt number. The Nusselt number is a dimensionless number describing the strength of heat transfer. For a complete oven, the Nusselt numbers of each nozzle are combined to model the effect of the air flow on the solid object, and thereby model the heating. Furthermore, we utilize and extend the novel geometric routines in IBOFlow that efficiently and robustly compute the intersection between a triangular volume mesh and a hexahedral Cartesian mesh [6]. This allows us to accurately describe the solid geometry on a coarse background grid, and enables the efficient solution of the heating of objects inside the oven. The novel algorithm and separation of scales approach allow us to simulate on a standard workstation. This is in contrast with previous work on simulation of oven curing [7], where a Lattice Boltzmann solver in the fluid is coupled with a finite difference solver in the solid, which requires a large cluster to run. The simulation results are demonstrated to be in close agreement with measurements from automotive production, and they can also be utilized for multicriteria optimization [1]. The structure of this paper is as follows: In Sec. 2, a detailed description of the numerical model is given, including the equations, the separation of scales approach, and the conductive heat transfer model. In Sec. 3, the curing of a truck cab in an automotive oven is described in detail. Comparisons between simulation and measurements are given and analyzed. Finally, conclusions are given in Sec. 4.

The proposed numerical method is motivated by the fact that a complete time and scale resolved simulation using the Reynolds' averaged Navier–Stokes' equations together with conjugated heat transfer is very computationally demanding. This is especially true since our goal is to present a method where an entire curing process, up to one hour, can be simulated overnight on a standard workstation. As will be described in more detail below, the necessary time-step to resolve the impinging jet from a typical nozzle in an automotive curing oven is as small as $10−4 s$⁠, which would give more than 10⁷ time steps in a full oven. To resolve the boundary layers at the moving target, cell sizes less than $0.1 mm$ are necessary, whereas a full oven can be 100 m long, and have hundreds of nozzles. In this section, we will describe how we solve this problem by separating the scales while preserving a physics-based approach, localize the resolved simulations, and couple the localized simulations to the full oven scale. The idea is that the smallest time and length scales are only resolved for the localized simulations that are decoupled from the macroscale simulation. Detailed simulations for a range of physical scenarios are precomputed and stored in a database.

The numerical method has been implemented in the in-house multiphysics solver iboflow [8], extending earlier available software modules employed in the Virtual Paintshop. The fluid dynamics engine in IBOFlow is a colocated, segregated, incompressible Navier–Stokes solver on an octree-based Cartesian mesh, which uses the SIMPLE-C method for pressure–velocity coupling. All geometries are handled with help of an immersed boundary method, for further details, see Refs. [1,4,9], and [10].

The Reynolds' number, with air, from a nozzle in a typical automotive curing oven is in the range $25,000−100,000$⁠, well into the turbulent regime. Nozzles can be either circular or rectangular, and are used to increase the heat transfer at the surface. In Ref. [11], a comparison between the applicability of different turbulence models to estimate the Nusselt number for impingement heat transfer is performed. The recommendation is to use either Menter's $k−ω$ SST or Durbin's $v2f$ method. We use the SST turbulence model [12–14], which has lower computational cost and still captures the location of the secondary peak well. The secondary peak, which corresponds to the maximum heat transfer coefficient, can be seen in Fig. 1 and is a typical characteristic of the Nusselt number below a round jet [11,15,16]. An immersed boundary formulation of the SST model is given in Ref. [14], using the variable transformation $g2=1β*ω$⁠, which we use in our implementation. The motivation to use g instead of ω is that $g→0$ linearly in the viscous layer, and the immersed boundary method linearizes the solution variables toward the boundary. That is, in practice we use a k–g model, not a $k−ω$ model. The heat flux at the solid fluid interface is computed from the friction temperature and velocity. The approach is similar to the one in Refs. [15] and [16].

where $uτ$ and $Tτ$ are the friction velocity and friction temperature, respectively.

According to experimental correlations, see overview in Ref. [20], Nu scales with an exponent of the Reynolds' number. For a laminar jet the exponent is 0.5, but for the turbulent case most authors use a higher value from the range $0.5−0.8$⁠. We use the exponent 0.56, which was determined to work well in the present setting after testing.

Our approach is based on separation and localization. We localize the simulations to individual nozzles to allow us to separate the boundary layer scale from the oven scale. The scale separation contains three steps: motivation, local description, and local to global coupling.

To allow us to localize the heat transfer, we first have to investigate how two impinging jets interact for typical nozzle size, Reynolds' number, and nozzle to target distance. We therefore perform a test simulation with two nozzles with 10 cm diameter separated by 40 cm located 40 cm from a steel plate. The setup is shown in figure. These are typical sizes and distances for a curing oven. The simulations are performed at Reynolds' number 65,000. Two simulations are performed: first the Nusselt number profile for one nozzle is computed, second the Nusselt number profile of both nozzles is computed. We compare the results of simulating both nozzles at the same time with using the one nozzle simulation twice, but shifting the results. The comparison is shown in Fig. 2, where it is clearly seen that for a typical curing oven configuration the heat transfer impact of the individual nozzles is independent within the accuracy of the Nusselt number simulation [11,20]. The peak at r = 0 for the two nozzle setup can be ignored since it represents a small area, i.e., a small part of the total heat transfer. This allows us to restrict our attention to one nozzle at the time and simulate them independently. When two or more Nusselt number profiles impact the same region we set the local heat transfer coefficient to their maximum.

There are many experimental and numerical studies on the simulation of heat transfer from an impinging jet, see, e.g., Refs. [11,15,16], and [20] and references therein. In Ref. [20], a thorough review of experimental correlations is given.

To compute Nusselt number profiles, we use a setup similar to Ref. [16], but instead of using a fixed flux at the lower domain boundary we set isothermal conditions. In contrast to the approaches in Refs. [15] and [16] that are similar to ours, we do not use rotational symmetry. This allows us to study rectangular and circular nozzles in the same framework, and also to perform nozzle interaction studies as above. Large rectangular nozzles are common in curing ovens, especially in the hold zone.

As noted by Refs. [15] and [16], the specific solid material has small influence, especially at the high Reynolds' numbers studied here; in curing ovens, the Reynolds' number is usually in the range $25,000−100,000$⁠. The influence of the solid material, as investigated in Ref. [16], is smaller than the expected accuracy of a Nusselt number simulation with SST; in Ref. [11] the accuracy at the stagnation point is estimated to be around 20%.

In our simulation setup, we position a 5 cm-thick steel plate with a $2.4×2.4 m2$ cross section area at the bottom of the domain. We set isothermal boundary conditions at the lower boundary of the plate. At the top of the domain, the nozzle is put in the middle, for the rest of the upper domain boundary, we use adiabatic slip wall, to minimize the influence and the mesh requirement. The remaining domain boundaries are pressure outlet. The simulations are performed with y⁺ approximately 1 at the plate, and CFL = 1 corresponding to a time-step in the order of 0.1 ms. A sketch of the setup is shown in Fig. 3.

To ensure that the grid resolution described above is sufficient, we perform a grid study of the Nusselt number profile under a nozzle with 10 cm diameter at Reynolds' number 65,000 and nozzle to plate distance 40 cm. The result is shown in Fig. 4. The ratio between the grids in the wall normal direction at the wall is 2, corresponding to y⁺ bounded by 2, 1, and 0.5, respectively.

The results for 7 cm and 10 cm diameter nozzles at Reynolds' numbers 45,000 and 65,000, respectively, are shown in Fig. 5. These results are at 30, 40, and 50 cm distance. For each diameter, we generate Nusselt number profiles for a range of distances and store all the results for varying diameters and distances in a database. The chosen Reynolds' numbers are in the middle of the range typically used and are directly applicable for the case described in Sec. 3.

To describe an oven, we set the location and orientation of the nozzles. To set a specific operating state of an oven we assign the mass flow rate to each nozzle, and the temperature to each zone in the oven. A zone consists of many nozzles, often grouped together in a periodically repeating pattern. The zones are typically distinct, and represent ramp zones, where the mass flow rates are high and the object is heated rapidly, and hold zones, where the mass flow rates are lower.

To determine the heat flux on the object at any time during the curing, we use ray-tracing. For each node of the surface triangulation of the object, we shoot a ray toward each nozzle to determine if it is visible from that nozzle. If it is visible, we compute the distance and angle, this allows us to apply the Nusselt numbers computed in Sec. 2.3.2. For each node, based on the previous analysis of nozzle interaction, we compute the maximum heat transfer coefficient over all visible nozzles. The reference temperature in Eq. (11) is allowed to vary with position in the oven, using either the nominal set temperature or a reference measurement. Note that this temperature can be significantly different in different parts of the oven, e.g., the heating and cooling zones. The corresponding heat flux is applied to the solid.

The body of a car or truck cab consists, to a large extent, of $0.5−2 mm$-thick sheet metal. In order to resolve such a geometry on a Cartesian grid, we use the volume fraction method developed in Ref. [6], which allows us to describe the volume and area fractions locally. The method only needs a surface mesh of the object. The approach is based on computing a least squares plane approximation of the surface in each cell, and computing the volume and area fractions from the resulting polyhedra. In Fig. 6 an example of an object with two sheet metal parts is shown together with the Cartesian mesh and the volume fractions representing the object.

To improve the accuracy of the volume and area fraction computation, it has been significantly improved and extended from the base approach described in Ref. [6]. Rather than computing the volume fractions directly on the grid, we refine it recursively during the computation. The refined cells themselves are never stored, instead their volume and area fractions are added to the cell above them in the recursion. The recursion is stopped as soon as the two opposite surfaces of a sheet metal part are in different cells or a maximum recursion depth is reached. For refined cells that do not intersect the object, we use ray-tracing to determine whether they are inside or outside the object. The convergence of the volume for the cab studied in Sec. 3 is shown in Fig. 7. For the results, we use recursive level 6.

For multimaterial parts the materials are averaged in the cells and on the faces. This allows us to accurately describe the volume averaged density and specific heat capacity of each cell, and the area-averaged thermal conductivity on each face of each cell. This is sufficient to accurately describe the heat transfer inside the solid objects. The area and volume fractions correspond to the variable α_s in Eq. (7). The time scale necessary to resolve the conductive heat transfer with given fluxes is on the scale of seconds, drastically decreasing the number of time steps necessary to simulate a curing process.

Those parts of the surface mesh that are within the bound are marked as interior, all other parts are marked as exterior.

Finally, we add the conductive fluxes (13) to the interior triangles.

To validate our approach, we simulate the curing of a Scania R20H cab in a convective curing oven at the paint shop in Oskarshamn, Sweden. A drawing of the oven is shown in Fig. 8 and its discretized representation in Fig. 9. It has 306 circular and 24 rectangular nozzles. The measurements are performed on a dry cab with eight probes attached to it. The location of the probes on the cab is shown in Fig. 10. Probe 1 is an air probe used to track the air temperature inside the oven. For the validation, we compare simulation results with the measurements of probes $2−8$⁠, which are positioned on the cab itself. The probes are positioned to give an accurate description of the heating of the cab, including areas such as beams with thicker material. To ensure proper curing, the resulting oven curves should match the specification of the paint manufacturer. In particular, the minimum time above paint specific critical temperatures must be ensured. This time can be postprocessed from the results of the simulation.

The oven is 47 m long and the total curing time, including cooling, is 40 min. The maximum temperature in the oven is 190 °C. The ramping zone contains 110 10 cm circular nozzles and 24 $0.456×0.915 m2$ rectangular nozzles. The cooling zone contains 196 7 cm nozzles. In the ramp zone, the nozzle Reynolds' number for the circular nozzles is 88,000 and in the cooling zone the nozzle Reynolds' number is 120,000. A drawing of the oven is shown in Fig. 8 and the discretized version in Fig. 9. In the discretized version, the elevators are modeled as horizontal zones. The eight zones in the oven with different flow and temperature conditions are shown in Fig. 9.

To simulate the curing in the discretized oven, we use nozzle Nusselt number profiles generated as in Sec. 2. The conductive heat transfer is simulated on an isotropic Cartesian grid with 1,629,538 cells with size 6.25 mm. The volume was computed using six recursive levels as motivated by the grid study in Fig. 7. All triangle-cell intersection areas are computed and marked for flux type given by Eqs. (11) or (13), split at 5 cm, and Eq. (15). The simulation was performed with the time-step 1 s. The large time-step is possible due to the multiscale approach. The detailed simulations are one-way coupled to the macroscale; therefore, the time-scales decouple. The results of the simulation compared with the measurements are shown in Figs. 11–17. As seen in the figures, the simulations closely capture the temperature profiles from the measurements. The pointwise mean deviation between measured and simulated temperatures for the seven probes is shown in Table 1. The least accurate probe (3) is positioned on a thin part in the front of the cab, where position is important, and the projected fluxes give a larger error compared to a resolved CFD simulation. For the measurements, the pointwise values can vary up to 10°C between measurements, i.e., more than $5%$⁠. On average the results, except for probe 3, are within the measurement accuracy.

In this paper, a novel framework for simulating the heat transfer inside convective curing ovens is presented. The framework combines localized resolved RANS simulation coupled with background grid conductive simulation using projection of Nusselt number profiles. It also includes radiosity based radiation fluxes. Utilizing a background grid decreases the necessary number of degrees-of-freedom for the conduction simulation by several orders of magnitude. The localization to a single nozzle at the time alleviates the need for a multigrid approach and a large computational cluster. The very efficient implementation gives a major improvement of computational speed compared to earlier approaches, and makes it possible to perform detailed simulations in close to real-time on a standard computer. This fact allows the engineers to replace physical prototypes with virtual ones to shorten the lead time in product development, to reduce energy consumption, and to avoid future unforeseen technological and environmental problems that can be extremely costly if they are discovered at the end of the production line, or even worse by the customer.

A validation is performed for a truck cab cured in an oven at the Scania plant in Oskarshamn, Sweden. Overall, the agreement between simulations and measurements is very good, almost within the measurement uncertainty. The conclusion from this and other performed case studies is therefore that the simulations can be used to predict the outcome of the process, optimize process parameters, and detect areas with insufficient curing. The framework is integrated in the ips software¹ as an oven simulation module, complementing the other virtual paintshop tools. The standard tests carried out by the automotive manufacturers are on dry objects, which is consistent with the recommendations from the paint manufacturers. Our initial goal has therefore been to validate such tests as demonstrated in this paper. The natural next step is to include transient tracking of the paint layer thickness and solvent concentration, to allow simulations of the curing itself. Future work also includes to analyze the effect of the oven curing on the adhesive joints from hemming processes.

This work has focused on convective ovens. To extend the approach to also include IR ovens, that are commonly used in repair shops, would be straight forward. The radiosity fluxes already included would be complemented by radiation sources. The impact of these radiation sources would be based on their view factors.

This work was supported in part by the Swedish Governmental Agency for Innovation Systems, VINNOVA, through the FFI Sustainable Production Technology program and the projects” Virtual PaintShop—Simulation of Oven Curing” and “Virtual Verification of the Hemming Process,” and in part by the Production Area of Advance at Chalmers University of Technology. The computations were partly enabled by resources provided by the Swedish National Infrastructure for Computing at Chalmers Centre for Computational Science and Engineering (C3SE) partially funded by the Swedish Research Council through Grant Agreement No. 2018-05973.

• VINNOVA (Grant No. 2016-03371; Funder ID: 10.13039/501100001858).

• Swedish Research Council (Grant No. 2018-05973; Funder ID: 10.13039/501100004359).

$cfp$ =

fluid specific heat capacity (J/kgK)

$csp$ =

solid specific heat capacity (J/kgK)

D =

nozzle diameter (m)

F =

view factor

g =

turbulent time scale (⁠$s$⁠)

h =

heat transfer coefficient (⁠$W/Km2$⁠)

J =

radiosity (⁠$W/Km2$⁠)

k =

turbulent kinetic energy (J/kg)

Nu =

Nusselt number

p =

pressure (Pa)

P_k =

production term (J/s)

Pr_t =

turbulent Prandtl number

q_i =

heat flux (⁠$W/m2$⁠)

q_r =

heat radiation flux (⁠$W/m2$⁠)

q_w =

wall heat flux (⁠$W/m2$⁠)

$qexts$ =

external heat source (⁠$W/m3$⁠)

$Pω$ =

production term (⁠$kg/s2$⁠)

t =

time (s)

T_f =

fluid temperature (K)

T_ref =

reference temperature (K)

T_s =

solid temperature (K)

T_w =

wall temperature (K)

$Tτ$ =

friction temperature (K)

$u$ =

fluid velocity (m/s)

$uτ$ =

friction velocity (⁠$1/s$⁠)

α^s =

solid volume fraction

ε =

thermal emissivity

κ_f =

fluid thermal conductivity (W/mK)

κ_s =

solid thermal conductivity (W/mK)

μ =

dynamic viscosity (Pa·s)

μ_t =

dynamic eddy viscosity (Pa·s)

ρ_f =

fluid density (⁠$kg/m3$⁠)

ρ_s =

solid density (⁠$kg/m3$⁠)

σ =

Stefan–Boltzmann's constant (⁠$W/m2K4$⁠)

ω =

specific dissipation rate (1/s)