Intro to Modeling Evaporative Cooling

Nancy Bannach December 8, 2014
Share this on Facebook Share this on Twitter Share this on Google+ Share this on LinkedIn

When you think of evaporation, you probably think of the cup on your desk that spreads the aroma of coffee or tea. But evaporation is also a process with many industrial and scientific applications, ranging from meteorology to food processing. This blog entry is the beginning of a new blog series on modeling evaporative cooling. Here, we introduce the basic concepts using your coffee cup as an example.

Some Basic Concepts of Modeling Evaporation in COMSOL Multiphysics

Evaporation is a process that occurs if some liquid vaporizes into a gaseous phase that is not saturated with the liquid. We exemplify this process and its characteristic properties by using water as the liquid and air as the gaseous phase.

Let’s first define the saturation pressure, p_{sat}, at which the thermal equilibrium with the liquid or solid state is reached. It is strongly temperature dependent and there are many approximations out there, which are all very similar but not exactly the same.

The COMSOL Multiphysics simulation software uses the following approximation from Principles of environmental physics by J. L. Monteith and M. H. Unsworth (1990):

(1)

p_{sat}(T)=610.7 Pa \cdot 10^{7.5 \frac{T-273.15K}{T-35.85K}}

For ideal gases, it is easy to determine the saturation concentration at which the relative humidity is 100% with:

(2)

c_{sat}=\frac{p_{sat}(T)}{RT}

where R is the ideal gas constant.

The thermodynamic properties of moist air depend on the fraction of water vapor. A mixture formula is used to describe the properties with the proportional amount of dry air and water vapor. Assuming air behaves as an ideal gas, the density reads:

(3)

\rho_m=\frac{p}{RT}\left(M_a X_a+M_v X_v\right)

More details and references about the implementation of the moist air properties as used by COMSOL Multiphysics can be found in the Heat Transfer User’s Guide (located within the Heat Transfer Module).

Modeling Evaporative Cooling: A Coffee Cup Example

Before setting up the COMSOL Multiphysics model, let’s consider the effects causing the cooling of the coffee as it evaporates.

We assume a slight air draft around the cup (or beaker, since there is no handle in this case) that accelerates the cooling by transporting heat and removing water vapor from the surface. At the coffee-air interface, vapor escapes from the liquid into the air, causing additional cooling by evaporation.

Modeling evaporative cooling in a coffee cup.
Sketch of the participating effects in a coffee cup.

How to Implement the Evaporative Cooling Effect in a Model

The first step is to make use of the symmetry, which reduces the model size and thereby the computational time. For the slight air draft, we use the Turbulent Flow interface with a constant air velocity. A reasonable approximation here is to assume that the flow field won’t change with a change in temperature and moisture content. Hence, we calculate a stationary velocity field in an initial study.

What else do we need to model the evaporative cooling effect?

Thanks to the predefined moist air fluid type, it is easy to implement the evaporative cooling effect in a COMSOL Multiphysics model. To determine the temperature field, we add the Heat Transfer in Fluids interface to the model, whereupon the Multiphysics node appears.

With the Multiphysics node, you can build your non-isothermal flow model sequentially. That’s precisely what we’ll do here. We’ll start with the Turbulent Flow interface and add the multiphysics coupling one after the other. The Non-Isothermal Flow node defines the two-way coupling between the flow and the heat interfaces. Note that in this case, we do not need the strongly coupled approach since the flow field is assumed to be independent of the temperature or moisture content. Using the properties from the flow interface, the Non-Isothermal Flow node also accounts for the turbulence effects in the heat transfer interface.

Node settings defining non-isothermal flow properties.
Multiphysics node for Non-Isothermal Flow. The node settings define the non-isothermal flow properties: a common density for the heat transfer and flow interface, a turbulence model for heat, flow heating, and the names of the interfaces.

Let’s have a closer look at the steps for modeling evaporation. First, set up the coupling between heat transfer and vapor transport in order to accurately model the evaporative cooling effect and utilize the postprocessing variables that come along with the Moist Air feature of the Heat Transfer Module, such as relative humidity or moisture content.

Heat transfer settings within the air domain.
Settings for heat transfer inside the air domain: (1) Coupling of the flow field (which is done via the Multiphysics node) for convective transport of the moist air. (2) Coupling to the Transport of Diluted Species interface, which gives us the correct input quantity of water vapor (3) to determine the moist air properties according to Equation 2.

The last thing to do is to set up proper boundary conditions. Here, we only discuss the boundary conditions connected to the evaporation. The rest is straightforward and can be read in the documentation of the model.

At the water surface, heat is released by evaporation. The heat of vaporization for water is approximately H_{vap}=2454\frac{kJ}{kg}\cdot M_w (it’s actually temperature dependent, but using a constant value here is a reasonable approximation) with the molar mass of water M_w=18.015\frac{g}{mol}. The amount of heat released depends on how much vapor escapes from the water surface into the air. This relates the heat source to the Transport of Diluted Species interface via the total flux in normal direction to the surface, which can be understood as the net flux of water vapor into air.

At the water surface, the relative humidity will always be 100%. Hence, the saturation concentration is reached. It defines the concentration of water vapor at the water surface according to Equation 2, with the saturation pressure determined by the Heat Transfer in Fluids interface. All in all, it is a strongly coupled phenomena that is implemented in no time.

Next, we take a look at the results of a time-dependent study over 20 minutes. The initial coffee temperature was 80°C and air at 20°C enters the modeling domain with a relative humidity of 20% causing the cooling. Below, you can see the resulting temperature and relative humidity distribution, both after 20 minutes.

A model depicting temperature distribution in the coffee cup.
Temperature distribution after 20 minutes.

An image showing relative humidity.
Relative humidity after 20 minutes.

The temperature is highest in the shadow zone of the beaker/coffee cup. Therefore, the relative humidity becomes very low.

Does evaporation have a strong influence on the cooling? We can find out by comparing the average temperature of the coffee — including evaporation — to that of the same model neglecting evaporation.

To do so, we’ll set up a third study solving for the Heat Transfer in Fluids interface only and disable the Boundary Heat Source node. The resulting plot clearly shows that cooling due to evaporation has a significant impact on the overall cooling:

Plot comparing average temperature over time.
Comparison of the average coffee temperature over time.

Conclusion and Next Steps

This blog post has shown the basic aspects you need to consider when modeling evaporative cooling. Keep these concepts in mind as we continue the series about evaporation modeling. Next up, we’ll take it a step further and answer the question of what happens in a porous material and how this can be implemented.

For now, feel free to download the Evaporative Cooling model shown here along with detailed instructions from our Model Gallery to try it out for yourself.



Comments

  1. Morteza Alipanah October 31, 2016   4:38 pm

    Hi,

    The way of how to simulate the evaporating water is explained here. can you please let me know how can we simulate the condensing water vapor mixed in the moist air into a liquid.

    Thanks

  2. Nancy Bannach November 4, 2016   6:13 am

    Dear Morteza,
    the way of modeling condensation is exactly the same. If condensation or evaporation occurs depends on the temperature. This determines the saturation pressure/concentration and hence the heat source/sink term. Just set the relative humidity at the air inlet to 1 and set a higher ambient temperature and lower water temperature.

    In this example we neglected the decreasing of the water level in the beaker due to the relatively small amount of water evaporating into the air. If you want to include such effects (also rising due to condensation) you can add a deformed mesh physics interface to track the water level according to the amount of water entering. This adds more complexity and numerical effort to your model. You can have a look at the following example of freeze drying:
    https://www.comsol.de/model/freeze-drying-3924

    If you need further assistance feel free to contact our support team https://www.comsol.de/support.

    Have a nice day,
    Nancy

  3. Omar Al-Rawi February 6, 2017   9:21 am

    Hi,

    Many thanks for this intensive explaination.

    You have used the Lagrange multiplier expression, c_lm, to represent the flux from the water surface. Does this expression include the convection flow effect in air in addition to the diffusive effect and can I check that? Also, why this expression gives results lower than if I define the flux as:
    Flux = -D*concentration gradient

    Thanks

  4. Nancy Bannach February 7, 2017   5:27 am

    Dear Omar,
    this is a kind of numerical question. I’m not a mathematician, so I explain it with my own words. At the surface you have a fixed concentration (Dirichlet boundary condition). In the equation COMSOL solves (weak form), the Dirichlet boundary condition vanishes, but c=c0 holds of course.
    Now, you only know the concentration at the boundary, but not the flux. So using the concentration gardient is more an approximate flux since you calculate this gradient with the solution for c inside the domain.
    With the Lagrange multiplier you add additional equations for the flux on this boundary, which answer the question which flux must be applied to get the concentration c=c0. And then c_lm gives the exact flux. Often this difference is small, but especially in multiphysics problems you should use this approach. Some interfaces have this active by default. This is the “Compute boundary fluxes” check box in the Interface settings “Discretization” section.

    The variable c_lm is the total flux. But since the velocity is zero at the ware surface (no-slip) wall boundary condition the flux is by diffusion. Of course, the convection in the air affects this flux, because heat and concentration is carried away and this affects the saturation concentration.

    I hope this answers your question.

    Have a nice day,
    Nancy

  5. Omar Al-Rawi February 24, 2017   5:57 am

    Hi Nancy,

    Many thanks for this explaination.

    Indeed, I have another question. Based on the Reynolds number you’ve mentioned (1500), which I don’t know how is calculated, the flow should be laminar. So, why did you impose the model with turbulent flow?
    My wondering is becuase that I’m modelling evaporation of two neighbouring droplets under forced convective environment inside a duct with size (W20, H10 and D20mm). I’ve imposed the air flow is modeled with laminar flow, and the air flow velocities are (0.02, 0.2 and 2m/s). Just with the last air flow velocity, the model doesn’t work although the flow is still with laminar regime. Therefore, I am asking that becuase that I have to change the flow regime from laminar to turbulent with low Re number.

    Regards,
    Omar

  6. Nancy Bannach February 24, 2017   7:24 am

    Hi Omar,

    unfortunately there is not one single critical Reynoldsnumber where flow turns from laminar into turbulent. The critical Reynoldsnumber strongly depends on the geometry. If you would remove the cup from the flow domain, then it could be laminar. But if you place an obstacle in a flow stream, you almost always introduce turbulence.
    The calculation of the Reynoldnumber uses a characteristic length. For a pipe, this is the radius or diameter and it is just a convention what you use. Hence, the value for the critical Reynolds number changes.

    Maybe another blog post is interesting for you (coincidentally I was the autor again ;-) ):
    https://www.comsol.com/blogs/characterizing-flow-choosing-right-interface/

    In general fluid flow is a very challenging topic and still a huge topic in numerics. For your specific problem I recommend to use our Support Center or Discussion Forum, which is dedicated to talk about individual tasks.

    Have a nice day,
    Nancy

  7. Israr Ahmad April 3, 2017   3:52 am

    Hi Ma’am
    I have been working on electronic cooling by electrowetting on dielectric (EWOD) technique. where droplet is manipulated through sequential actuation voltage over the electrode. However, the droplet is transported over the hot spot where constant heat flux has been provided. The droplet is supposed to carry the heat via conduction, convection and phase change and moved over another electrode. Till now, i have been able to solve single phase cooling process with laminar two phase flow phase field and heat transfer in fluid module but facing difficulty to model evaporative cooling.
    In this example, the draft of air itself over the surface transported the heat and vapor into it. But in my case, there is still air in between two plates. Can you please give me idea for modelling evaporative cooling in my model.

    With Regards

  8. Bridget Cunningham April 3, 2017   9:40 am

    Hello Israr,

    Thank you for your comment.
    For questions related to your modeling, please contact our support team.

    Online support center: https://www.comsol.com/support
    Email: support@comsol.com

    Best,
    Bridget

  9. Victor Casquero Pretell April 28, 2017   10:17 pm

    Hi Nancy,

    The K value for the transport of heat in the evaporation of water, was assumed value 100 or it can be estimated mathematically.

  10. Nancy Bannach May 2, 2017   4:03 am

    Dear Victor,

    I assume your question about K relates to the 5.3 version of the model using the Moisture Transport in Air interface. The approach discussed here is a little different.

    You can find an explanation about K in the related pdf document:
    “The evaporation rate must be chosen so that the solution is not affected if further
    increased. This corresponds to assuming that vapor is in equilibrium with the liquid.”

    If this does not answer your question, contact our support team for detailed discussions (support@comsol.com).

    Best regards,
    Nancy

  11. Omar Al-Rawi June 12, 2017   8:12 am

    Hi again Nancy,

    My question may not relate to this blog but hopefully you can help me.

    If I want to send your model, which involves three studies, to HPC, which command I can use into the shell file to run all these studies?

    Actually, I have used this command to run one study so will it be different in the case of your model?
    “comsol batch -np $NSLOTS -tmpdir $TMPDIR -inputfile input_model.mph -outputfile output_model.mph”

    Regards,
    Omar

  12. Bridget Cunningham June 12, 2017   9:16 am

    Hello Omar,

    Thanks for your comment.

    In regards to this question, please contact our Support team.

    Online Support Center: https://www.comsol.com/support
    Email: support@comsol.com

  13. Omar Al-Rawi August 29, 2017   5:46 am

    Dear All,

    I have other questions about this model, particularly the heat transfer interface.

    1- In the hea tansfer in fluid 1 subnode, the “Equivalent conductivity for convection” feature is activated, and I think is because the convection inside water is neglected. But the question is that base on what the temperature difference is defined as 3[K], and what does this parameter mean?

    2- In the heat transfer in fluid 1, 2 and solid, why did you select the temperature in the “Model input” related to the “ht” interface instead “nitf1″ multiphysics node?

    Regards,
    Omar

  14. elias Yazdanshenas September 19, 2017   11:31 am

    Dear all,
    I am looking for how to change the fluid phase in the condenser by losing heat. Please let me know if you have any idea about that.
    Regards,
    Elias

  15. Nancy Bannach September 20, 2017   5:56 am

    Dear Omar,
    I’m not sure to which version of the model you refer. If I open the model from the gallery the Temperature difference for Equivalent conductivity is set to “Automatic”. Otherwise a user-defined temperature difference can be a guess or maybe estimated from similar problems. The temperature difference determines the Rayleigh number which again is used to estimate the convective heat flux by Nusselt number correlations (check the manual for more details).

    About your second question (I refer to the version 5.3 model from the gallery):
    The temperature from nitf is not available for solids and fluid1, because there we don’t solve for Nonisothermal flow. Otherwise, in most cases the temperatures from nitf and ht are the same and you can rely on the automatic seetings.
    There could be a difference depending on how you build your model and studies. If you are unsure, you can plot nitf.T-ht.T. This should be zero.

    Dear Elias,
    I’m not sure if I understand your question. In general, if evaporation/condensation occurs depends on the temperatures and is calculated automatically. If you look for other phase changes, then the “Phase Change Material” feature could be interesting.
    But I recommend to contact our support team for a detailed discussion

    Online Support Center: https://www.comsol.com/support
    Email: support@comsol.com

    Best regards,
    Nancy

  16. Steven Spielman January 3, 2018   8:14 pm

    Does the “evaporation rate” parameter have any physical significance? Why does it have the name it does, rather than, say, “Mangione coefficient” or “turbo encabulator”. I’m serious. The documentation I’ve seen would make exactly as much sense, had the name been replaced with either of these options.

  17. Nancy Bannach January 4, 2018   3:09 am

    Dear Steven,

    “Evaporation rate” is a common name that is used to describe drying processes and many other processes that deal with evaporation, e.g. food processes.
    The evaporation rate is a property that depends on the material properties and the process itself.
    However, in this case it is just chosen to be large enough such that the approach is valid. This doesn’t mean that it has no physical meaning but the exact value doesn’t play a role for the physics described here as long as it is large enough.

    See also:
    https://www.comsol.de/model/evaporation-in-porous-media-with-small-evaporation-rate-21931

    Have a nice day,
    Nancy

Loading Comments...

Categories


Tags