Simcenter STAR-CCM+ De-icing and Defogging Simulations Using Fluid Film vs. Thin Film Best Practices

2022-12-15T21:56:59.000-0500
Simcenter STAR-CCM+ Simcenter STAR-CCM+ Virtual Reality Teamcenter Share Simcenter Cloud HPC Simcenter STAR-CCM+ Viewer Simcenter STAR-CCM+ Application Specific Solutions

Summary

This article shows the set up of Fluid Film and Thin film for defogging and de-icing simulations.


Details

For a long time, de-icing and defogging simulations in STAR-CCM+ were performed using the Thin Film model, which is based on empirical correlations for the evaporation of water. The Fluid Film model is physically more realistic, more flexible and more stable. This article will describe the setup of a headlamp defogging simulation both for the Thin Film and Fluid Film model and compare both approaches with each other.
Geometry Considerations
For this example, a generic car head lamp geometry is used. There are two light bulbs and several LED lights for the blinker and the daylight. The following picture shows a rendered image of the solid geometry:

Head Lamp Geometry Rendering
There is a small ventilation system that serves for the air exchange and the lights serve as a heat source. The transparent front cover is exposed to the environment and allows a transfer of heat. After a dark, cold night, the front cover has cooled down and a thin water film (consisting of millions of small droplets) has formed on the inside of the cover. Over the day, the environmental temperature will act as additional heat source. The components of the geometry are depicted in the following picture:

Head Lamp Parts

In special application cases you can also involve additional solar loads on the head lamp cover or an iced front cover. These topics are addressed only casually in this article.
Mesh Considerations
For this complex geometry, you may use the Polyhedral Mesher for the discretization of the volume. For the resolution of the boundary layers on the inside of the front cover, some prismatic cells will be generated by the Prism Layer Mesher. Mesh also the solid parts in order to account for the specific heat of the solid materials. As most of the solid parts can be considered as thin, you can select the Thin Mesher.

The meshing parameters for the current geometry are as follows:
Mesh Model Settings
Thin MesherCustomize Thickness Threshold: Enabled
Default Controls
ParameterValue
Base Size3 cm
Minimum Surface Size15 %
Prism Layer Stretching1.0
Thin Mesher Maximum Thickness5 mm
Prism Layer Total Thickness1.5 %
Custom Controls
Create Surface Controls in order to adjust the cell size distribution to local feature lengths and surface curvatures of specific surfaces. Add Part Controls in order to disable the Thin Mesher in fluid parts and to disable the Prism Layer Mesher in solid parts.

The mesh cells on the front cover surface have a edge length of about 2 mm. This size is too coarse in order to fully resolve the small droplets of the condensate which have a size of about 2 micro-meter. This will have an impact on the physical modeling, because surface tension effects need to be considered in a specific way. The resulting mesh will have about 327 000 cells. The following picture shows the cell surfaces on two cutting plane sections, colored by region (prism layer cells in red):
Head Lamp Mesh

Physical Properties
There are several physics continua to account for the different materials within the head lamp assembly. Whereas the solid physics continua and the halogen physics continuum differ from each other in the material and radiation properties, they remain the same both for Fluid Film and Thin Film approaches. For the setup of these two approaches, you just need to modify the Air physics continuum.

For the solid continua, select following models: Three-Dimensional, Implicit Unsteady, Solid, Constant Density and Segregated Solid Energy. For some of the physics continua you should activate the Radiation model (Surface-toSurface Radiation (S2S), Gray Thermal Radiation (GTR)). The following table provides an overview about the radiation model and material settings for each physics continuum:
Physics ContinuumRegionTypeMaterial Radiation
AirAir_Main, Air_ShroudFluidsee belowS2S, GTR
FilamentDB_Filament, FB_FilamentSolidTungstenS2S, GTR
GlassCoverSolidFused QuartzS2S, GTR
HalogenDB_Halogen, FB_HalogenFluidAir, Ideal GasS2S, GTR
PolycarbonatBlinker_LED, Daylight_LED (1-3), DB_Envelope, FB_EnvelopeSolidPC-
PolyethersulfonBlinker, Daylight, DB_Housing, FB_Housing, DB_Reflector, FB_Reflector, LeftBank, MiddleBank, RightBankSolidPES-
PolypropylenHousingSolidPP-
Using either the Fluid Film model or the the Thin Film model, the physics model selection for the Air physics continuum have to be changed according to following table:
Models for the Air Physics Continuum
(models in brackets will be selected automatically)
Fluid FilmThin Film
  • Three Dimensional
  • Implicit Unsteady
  • Multi-Component Gas
  • Non-Reacting
  • Segregated Flow (Gradients, Segregated Species)
  • Ideal Gas
  • Turbulent (Reynolds-Averaged Navier-Stokes)
  • K-Epsilon Turbulence (Realizable K-Epsilon Two-Layer, Exact Wall Distance, Two-Layer All y+ Wall Treatment)
  • Segregated Fluid Temperature
  • Fluid Film
  • Gravity
  • Radiation
  • Surface-to-Surface Radiation (View Factors Calculator)
  • Gray Thermal Radiation
  • Multi-Phase Interaction
  • Optionally: Solar Loads
  • Three Dimensional
  • Implicit Unsteady
  • Gas
  • Segregated Flow (Gradients)
  • Ideal Gas
  • Turbulent (Reynolds-Averaged Navier-Stokes)
  • K-Epsilon Turbulence (Realizable K-Epsilon Two-Layer, Exact Wall Distance, Two-Layer All y+ Wall Treatment)
  • Segregated Fluid Temperature
  • Gravity
  • Radiation
  • Surface-to-Surface Radiation (View Factors Calculator)
  • Gray Thermal Radiation
  • Thin Film
  • Defogging
  • Optionally: De-Icing, Solar Loads
Multi-Component Gas Model: Air, H2O
Material properties:
  • H2O > Heat of Formation: 0.0 J/kg
Defogging model properties:
  • Film Density: 997.561 kg/m^3
  • Film Latent Heat: 2441630 J/kg
Models for the Fluid Film phase:
  • Shell Three-Dimensional (Flow, Laminar, Segregated Fluid Film Temperature)
  • Multi-Component Liquid (Segregated Species)
  • Constant Density
  • Non-reacting
  • Optionally for De-Icing: Melting-Soldification
Multi-Component Liquid model of the Fluid Film Phase: H2O only
Material properties:
  • Heat of Formation: -2441630.0 J/kg
  • Saturation Pressure: Antoine Equation
  • Dynamic Viscosity: 100000 Pa s (explanation)
Multiphase Interaction models:
  • Film-Physics Continuum Interaction
Set the Fluid Film phase in the Film-Physics Continuum Interaction model to your Fluid Film phase and set the Connectivity in the Component Mapping model to both H2O. Then add another model:
  • Evaporation/Condensation
Evaporation/Condensation model properties:
  • Implicit coupling: Disabled
  • Under-Relaxation Factor: 0.2
  • N_seeds: 10E9 /m^2

Heat of Formation Settings (Important !)

For the Fluid Film approach, the settings for the Heat of Formation property has been changed. For the evaporated water component of the air it is set to zero and for the liquid water component of the film it is set to the negative value of the latent heat of evaporization of water at 25 C. This helps to stabilize the solution. By default, the Heat of Formation property represents the amount of enthalpy that is evolved during the formation of 1 kg of the material from it's elements in their standard states. The phase enthalpies are calculated by following equations:

enthalpy equation

In the equations above, hv and hl stand for the enthalpies of the vapor and the liquid phase, Tsat is the saturation temperature at which the phase change takes place, h0v and h0l represent the heat of formation, cpv and cpl is the specific heat and T0v and T0l are the standard state temperatures. The latent heat hlat is then calculated as follows:

latent heat equation

If the heat of formation for the vapor phase is set to zero at the standard state temperature, the heat of formation for the liquid phase becomes equal to the negative of the latent heat to fulfill this equation.

(see also Simulating Physics > Materials > Setting Material Properties Methods > Using the Enthalpy Difference Method for Heat of Vaporization)

Saturation Pressure


The saturation pressure psat of water is calculated using the Antoine equation as shown below. The coefficients A = 11.949, B = 3978.205 and C = -39.801 are fitted for the logarithm base e. They are valid for a temperature range between 0 C and 374.15 C. The units of the temperature T is in Kelvin and the atmospheric pressure patm is 101325 Pa.
antoine equation
(see also User Guide > Modeling Physics > Modeling Materials > Setting Material Properties Methods > Using the Antoine Equation)

Other Material Properties


You should ensure that the material properties for the air and the water components are adjusted to the temperature level on which the simulation is going to run on. The default settings represent the material properties for a temperature of 26.85 C. Especially for benchmarking the simulation results against experimental data, these settings may become more important and you may want to pay much more attention to have the correct properties for the correct materials.

Evaporation/Condensation Model Settings


Implicit coupling is deactivated, because in this kind of application the evaporation takes places far below the boiling temperature. Linearize Film and Gas Energy can make the simulation more stable. The number of seed points and the minimum seed radius is needed to compute the condensation at dry walls. A good estimation for the number of seed points on head lamp front cover is one billion per square meter.

(see also: User Guide > Modeling Physics > Modeling Multiphase Flow > Modeling Fluid Film > Modeling Fluid Film Evaporation and Condensation > Evaporation-Condensation Model Reference)

Modeling the Contact Force between Droplets and Cover


Two approaches for modeling of the contact force between the water droplets and the front cover surface will be discussed. First, there is a Surface Tension model available, when you selected the Fluid Film approach and it can be used in order to account for the capillary pressure (orthogonal to bottom surface) and line contact force (parallel to bottom surface) effects between the water droplets and the surrounding continuum within the fluid film momentum equation. The impact of these effects on the droplets in the real world is important, as it forces the small droplets to stay attached to the wall and thus mitigating the most portion of motion within the fluid. This is in contrast to the Thin Film model, in which the liquid is always fixed at the wall and cannot move.

However, due to the relative coarse mesh settings, single droplets of the condensate are not resolved. An averaged, constant film thickness distribution is used instead and with that simplified phase interface, the surface tension effect would not be reproduced correctly. Because there are several sources to accelerate the fluid in the film, such as gravity or shear stress on the phase interface, another way has to be found to mitigate the motion of the water film. This can be achieved by increasing the dynamic viscosity value of the water film material properties to a very high value. This approach enables virtually a resisting force in the water film opposing the outer forces and it does not affect the mass, momentum or energy balance.

Use the Surface Tension model in cases, where you expect larger drops (e.g. rain drops), or where the mesh is sufficiently refined.

(see also User Guide > Modeling Physics > Modeling Multiphase Flow > Modeling Fluid Film > Modeling Surface Tension > Surface Tension Model Reference)

De-Icing


Using the Thin Film model, you can optionally add the De-Icing model, respectively the Melting/Soldification model when using the Fluid Film approach. If the boundary or initial conditions are set properly, this will allow the water film to freeze and solidify or to melt over the time. You have the possibility to specify the ice material properties and the latent heat of the melting/soldification process. As we assume the film layer to be thin, we do not need the "Remove Soldified Film Option" like in airfoil icing applications. The benefit of the Fluid Film model is that de-icing and defogging are allowed at one boundary at the same time, whereas with the Thin Film approach, you can either simulate de-icing or defogging, but not both at the same time at the same boundary. Neither Thin Film allows a phase change from ice to water or vice versa.

(see also: User Guide > Modeling Physics > Modeling Multiphase Flow > Modeling Fluid Film > Modeling Fluid Film Melting and Solidification > Melting-Solidification Model Reference)
Domain and Boundary Conditions
For the Fluid Film model, the initial film thickness is an initial value which you can set in the physics continuum. When you use the Thin Film model, then the initial thickness is considered as boundary setting (see "Initial Conditions").

For the velocity inlet, keep the default velocity value of 1 m/s. This setting can be adjusted in case you want to account for convection and/or diffusion of the evaporated water. The static temperature of the inlet and outlet are left at the default values of 26.85 C. When you use the Fluid Film model, set the mass fractions of Air to 1.0 for both the inlet and outlet.

In the Filament and LED regions, you can enable a "Total Heat Source" in the "Energy Source Options". This will allow you to account for the heat emitted by the lamps.

For the Cover > Outer Cover boundary, set the "Thermal Specification" condition to "Environment". This allows a heat transfer from the environment over the front cover into the computational domain.
Initial Conditions
The initial temperature of the fluid film, the front cover and the air is 3 C. All other initial temperatures remain at the default settings (300 K). The initial film thickness is 5 m. In the Fluid Film model, the initial film thickness is defined in the physics continuum Air> Models > Fluid Film Phase > Initial Conditions.

In order to set the initial layer thickness when using the Thin Film model, set Regions > Air_Main > Boundaries > FluidFilmFront [Interface 1] > Physics Conditions > Thin Film Option to "Fog Layer" and in Physics Values > Fog Layer Thickness enter the value 5 m.
Analysis Controls
During the simulation run, you can monitor the balances of mass and energy. For this purpose, you have to consider the inflows and outflows of the domain as well as the different sources for phase mass and energy, such as:
  • Evaporated vapor from the water film
  • Latent heat of evaporization of the water film
  • Heat emitted by the lamps
  • Heat transferred over the cover
  • Heat transferred over the fluid/solid interfaces

Mass Balance

While the mass flux of the evaporated water in the Fluid Film model is stored in the "Film Evaporation Rate of H2O" field function, no such function is available for the Thin Film model. To track the evaporation mass flux in this case, use following setup:
  1. Create a user field function "Film Mass" with the definition: ${WindshieldFilmDensity} * ${FogLayerThickness} * mag($${Area})
  2. Create a sum report named "Current Film Mass" of field function "Film Mass" in FluidFilmFront [Interface 1]
  3. Create a new Field Sum Monitor named "Stored Film Mass" of the field function "Report: Current Film Mass" in FluidFilmFront [Interface 1], set the trigger to "Time Step" and activate "Enable Sliding Sample Window. Set the Slinding Sample Window Size to 2.
  4. Create a new maximum report named "Stored Film Mass" of field function "Sum of Report: Current Film Mass" in FluidFilmFront [Interface 1]
  5. Create a new expression report with the definition: (${StoredFilmMassReport} - 2 * ${CurrentFilmMassReport}) / ${TimeStep}. The current film mass has to be subtracted twice from the stored film mass, because the stored value contains the sum of the film mass of the last time step and the current time step. Dividing the difference by the time step gives the film evaporation rate.
With reports providing the values for the inlet and outlet mass flows and the film evaporation rate report, the mass balance can be reported by:

${MassFlowInletReport} + ${MassFlowOutletReport} - ${WaterFilmEvaporationRateReport}

You can create a monitor and plot from the mass balance report, as depicted in the following picture. For the Thin Film model, you can see some more "fluctuations" in the mass balance, resulting from the additional vapor fraction that is added to the domain during the evaporation. The Fluid Film mass balance is smoother, because the liquid film mass is already incorporated into the set of conservation equations.

mass balance

Energy Balance

The film evaporation rate report can also be used to track the latent heat of evaporization for the energy balance. Using the Fluid Film model, the latent heat of evaporization is incorporated in the boundary heat flux over the fluid film interface. The power of the lamps are applied to the corresponding solid regions, therefore the boundary heat flux over the fluid/solid interfaces will account for them. An expression for the power balance could be defined as follows:

${PowerInReport} + ${PowerOutReport} + ${PowerfromSolidsReport} + ${PowerFilmandCoverReport}

Using the Thin Film model, the latent heat of evaporization is not incorporated into the boundary heat flux and an additional expression has to be used to account for it. The defintions of the single reports are given in the table below:

${PowerInReport} + ${PowerOutReport} + ${PowerfromSolidsReport} + ${PowerCoverReport} + ${PowerFilmReport}
NameTypeField FunctionPart
Power InSurface IntegralBoundary Advection Heat FluxInlet
Power OutSurface IntegralBoundary Advection Heat FluxOutlet
Power from SolidsSurface IntegralBoundary Heat FluxAll Fluid/Solid Interfaces (without cover interface)
Power Film and CoverSurface IntegralBoundary Heat FluxFluid Film Interface
Power FilmExpression${FilmEvaporationRateReport} * ${FilmLatentHeatReport}
Power CoverSurface IntegralBoundary Heat FluxCover interface
The monitor plots of the power balance are depicted in the following picture. The energy balance of the Thin Film simulation is much more oscillating than in the Fluid Film simulation. This indicates a more stable convergence behavior for the Fluid Film approach.
energy balance

Post-Processing

Film Thickness Distribution

Most important engineering value for this type of analysis is the film thickness distribution. The field functions "Fluid Film Thickness" (Fluid Film model) or "Fog Layer Thickness" (Thin Film model) account for that.

You can display the current spatial distribution of the water film using a scalar scene. Use a copy of the water color bar with decreased opacity at the left end to clear out the very thin film thickness values. Following videos show the transient results of the different modeling approaches.

https://videos.mentor-cdn.com/support/videos/1200/6b39a9c3-d9d8-4688-adef-c9e6ef3e250a-en-US-video.mp4
https://videos.mentor-cdn.com/support/videos/1200/e1ccb357-463c-4a48-98c9-6aa15f9f3845-en-US-video.mp4

Some point probes were used to track the transient evaporation on distinct locations. This can be useful for comparisons with experimental data. The following plot shows the transient film thickness on seven locations as well as an averaged value for the Fluid Film model. This plot shows that the evaporation rate depends on the spatial position and can change over time.

Fluid Film Thickness

Cover Dry Area Fraction

You can define that a cell of the cover surface is considered as dry, if the fluid film thickness is below a certain value. Let's assume this threshold is about 1 m. Below this value, following vector field function returns the area of a cell, above the threshold, it returns a zero-vector:

${FluidFilmThickness} < 1e-6 ? $${Area} : [0,0,0]

Using a Sum Report, you can evaluate the sum of dry area magnitude in the fluid film region. Create another report to sum up the total area of the whole fluid film region. Then, create an Expression Report to return the dry area fraction. Use following expression to avoid floating point exceptions:

${CoverTotalAreaReport} > 0 ? ${CoverDryAreaReport} / ${CoverTotalAreaReport} : 0.0

The report value can be added as annotation to scenes or monitored and plotted.
cover dry area

Relative Humidity of Air

The evaluation of the relative humidity of the air using the Thin Film model is simply done by using the inline field function "Vapor Humidity Value". Such a field function is not available for the Fluid Film model, but it can be established as user field function by following steps:
  1. Go to Tools > Parameters and create scalars for the parameters A, B and C of the Antoine equation
  2. Create a new scalar field function "Partial Pressure Water", set dimensions to "Pressure" and the definition:
    ${AbsolutePressure} * ${MoleFractionH2O}
  3. Create a new scalar field function "Saturation Pressure Water", set dimensions to "Pressure" and the definition:
    exp(${Antoine_ParameterA} - ${Antoine_ParameterB} / (${Temperature} + ${Antoine_ParameterC})) * 1e5
    The multiplication with 1e5 is for the conversion of units bar into Pascal
  4. Create a new scalar field function "Relative Humidity [%]" with definition:
    ${PSat_Water} == 0 ? 9999 : ${PPart_Water} / ${PSat_Water} * 100
    This formulation will avoid floating point exceptions, if the saturation pressure of water should be evaluated as zero somewhere
This procedure is attached as JAVA macro to this article.

For visualization, create a new scalar scene and then create a new resampled volume of the air region and add it to the scalar displayer. Select "Relative Humidity [%]" as field function, or "Vapor Humidity Value", respectively. Adjust the minimum and maximum values of the scalar field to 0 and 100 and edit a copy of the color bar by setting the opacity of the left top marker to zero. This will clear the volume from low humidity values:
color bar settings

Add a solution time annotation to the scene and also a geometry displayer that you can make transparent. You can also change the Volume Rendering Settings in the scalar displayer to "Local Lighting" mode and increase the rendering quality to 0.75. The result will look similar like this:
relative humidity resampled volume

Performance Considerations

Setting the Time Step

The evaporation of water is a slow process, controlled by diffusion. Also the heat transfer processes are slow. If the inlet velocity is small and convection is not dominant, then you can select a larger time step of 0.1 s for the Thin Film model. The Fluid Film model appears to be more stable and allows also a larger time step of 1 s without changing the solution significantly.

Steady Simulation of Air Flow, Unsteady Evaporation

The evaporation of water does not affect the flow solution very much. Therefore it can be suitable, first to solve the steady state flow solution and afterwards switching to a transient simulation in which the evaporation of the film is solved and the flow solution is frozen. In this case, the time step can be increased which will result in a reduction of the solver time, but more manual work is required to switch the time model. Though keep in mind that a correct temperature field is key to retrieve a correct behavious for the condensation / evaporation processes.


KB Article ID# KB000040407_EN_US

Contents

SummaryDetails

Associated Components

Design Manager Electronics Cooling In-Cylinder (STAR-ICE) Job Manager Simcenter STAR-CCM+