Message Boards Message Boards


Model the Coffee cooling problem with WSM?

Posted 3 years ago
10 Replies
20 Total Likes

Greetings while I am new to WSM I have a few years physics experience.

It seems to me that a the classic 'when do I add milk/cream? is a perfect opportunity for illustrating the modelling capabilities of WSM. However, I am yet to see this application..

Too simple perhaps?

I have plenty of examples of mathematical modelling. this conference article covers all the possibilities with 5 models, each more complex to allow for fat globules on surface insulation radiation, colour of surface of cup, material of cup, volume of cup, surface area etc..

I am hoping experienced WSM guru will not find this too trivial and apply their talents to develop an all tim classic WSM model..

regards Gary high school physics/IT teacher melbourne Australia

** the need for this? we have just licensed every high school teacher (20k)and their students (200k) on whatever OS they have for the entire wolfram suite Wolfram Alpha PRO Mathematica Wolfram SystemModeler Wolfram Cloud...

10 Replies

This doesn't answer your question but I couldn't resist posting that a classmate of mine at Reed College, Andrew Case, wrote a B.A. dissertation modelling the behaviour of milk drops falling into a cup of coffee, and if I recall correctly, even did some of the work in Mathematica

BTW why not Mathematica? See: The Coffee Cooling Problem. You can download the code freely there.

enter image description here

Posted 3 years ago

Greetings Vitally i was aware of the Demonstration, however the differential equation and approximate exponential over estimates higher temperatures and underestimates lower temperatures. if you 'stand back' far enough observation data and the calculated data points are close, but a not very good fit.

an 'ideal' cooling curve does not represent real behaviour of a container. e.g.. radiation, convection, convection, insulation of milk fat layer

Better than a quadratic or quartic approximation however, not an accurate predictor of behaviour as the model does not take into account colour, container or surface interface (i.e. insulated with milk fat).

The demonstration is good as far as it goes.

I was after something like the kettle WSM example which calibrates the theory with real data to achieve a far better fit (~ 2%)

regards Gary (39?C here today!)

WSM can do this easily and is perfect for that type of system. I would start with one of the heat transfer examples and modify it. You can even add some of the thermal components if you want to model circulating air (i.e. someone turns on a fan!). Look at the thermal examples to see how to do that. The basic case is very much like the examples.

If you have problems you can post your example for some more help.


Hi Gary!

I went through some of the models in the paper you posted. Overall, it seems like SystemModeler could be very good fit for this, allowing students to add and remove different effects and see the results. The models I show below are attached to this post. I'd love to hear more about the milk adding problem also, seems like an interesting scenario, but I am not sure I understand what the object of the students would be.

As for the paper you linked:

Experiment 1 can be described using standard components from the Modelica.Thermal.HeatTransfer package. The pot will be HeatCapacitor component, the ambient temperature will be modeled using FixedTemperature and the convection is modeled using a ThermalConductor, which follows Newtons law of cooling.

enter image description here

The G parameter in the ThermalConductor is equivalent to the k parameter they use. From what I could tell, the paper did not include any measurement of the heat capacitance or ambient temperature so I went with 3 dl of water and 20 degrees Celsius. However, both of these would probably need to be higher to fit their experimental data.

To create experiment 2 I first duplicated experiment 1 by selecting it and pressing Ctrl+D, you can also right click and select Duplicate. Experiment 2 requires a component like the ThermalConductor but one that has an exponent that causes nonlinear behaviour in the heat flow. No such component exist in the Modelica Standard Library, but we can easily create one. I created a new component to be used in experiment 2 by dragging the normal ThermalConductor into Experiment2.

enter image description here

And gave it a new name, "ArbitraryExponentConductor"

Now I had to modify it to use the exponent. After opening the new component I first added a new parameter by right clicking the parameter view and and selecting Insert > Parameter

enter image description here

I used the name x as in the paper and used type Real.

enter image description here

Now I had to modify the equations so I went into the Modelica Text View (Ctrl+3) and changed the line:

 Q_flow = G * dT;


 Q_flow = G * dT ^ x;

dT corresponds to the temperature difference (tc-ts) in the paper.

Going back into Experiment 2, I changed the normal ThermalConductor by right clicking it and selected Change Type. In the dialog, I gave the name of the new type (CofeeCooling.Experiment2.ArbitraryExponentConductor). You can also drag the component from the component browser directly into the field.

enter image description here

or of course, delete the component, drag the new one in and make new connections.

For experiment 3, you need to add a some more stuff. Start by douplicating experiment 1. Connect the ThermalConductor to a new HeatCapacitor instead of the FixtedTemperature. That heat capacitor will be the pot, while the original one will be the coffee. The first ThermalConductor then represents equation 1 in the paper, transfer of heat from coffee to the pot. Add another ThermalConductor and connect it between the pot and the FixedTemperature to represent equation 5. Also add two BodyRadiation components and connect them from each capacitor to the FixedTemperature. These will represent all the radition effects described. They are bidirectional so they represent two equations each (3,4 and 6,7). For evaporation, I created a custom component which is described by the equation

  port.Q_flow = k * port.T;

Where k is the product of the P, l and v parameters described in the paper. You could add individual parameters for each of them instead, as described in the text for experiment 2.

Connect the evaporation to the coffee capacitor.

enter image description here

The 4:th experiment is much like the the 3:d one. I modified the evaporation component to have the equation

  port.Q_flow = k * port.T ^ z;

instead. However, with the exponent they present, 3.7, the temperature will go down very fast for all but very small values of k. Without knowing what parameter they had for l, v and P, it is hard to see if that approach is sound. Do you have an idea if they published more information on the other parameters they fitted?

Posted 3 years ago

Greetings Patrik

i follow the development of your far.

I get stuck on the next step.

where an amount of cold milk/cream is added to the coffee(assuming black) after a specified time usually t=0 and t=300seconds volume ~10ml Cp will be different to water . skim milk 3.97; 'regular' milk 3.77 compared to water 4.18 kJ/?C/kg

the total volume changes from 200ml to 210ml assuming complete mixing the 'new' temperature will decrease , however the mixed Cp is unknown so the standard calculation of ratios is difficult to use.

the contrasting treatments are to add milk at the beginning, or to wait 5 minutes (300 seconds) and add milk after allowing the coffee to cool and the milk to warm(!) ambient is 20?C

once a working model has been achieved, can the details be calibrated using 'real' data? similar to the WSM example of the cooling kettle?

Looking good, so far!

My approach for High school would be to research the history of Newton and differential equations. Identifying the scope and constraints on the 'ideal' solution.

Then introduce the concept of modelling and calibration. the main idea to highlight the benefits of theory and experience, where there is 'never the ideal', however there can be very close approximations.

in physics anything better than 5% is considered a 'good fit', 2% a 'close' fit and <1% always the goal.


Hi Gary,

I see two approaches you could take when adding the milk to the coffee. Either you could have everything collected into a single component that has a event in it corresponding to the addition of the milk, or you could have a separate component that specifies the addition of milk as a flow over time. I explored both those scenarios in the attached model. Would be very interesting to hear your input which one would be better from a teaching perspective.

Approach 1, with a discrete event in the coffee involves creating a copy of the HeatCapacitor component and adding some parameters that separates the heat capacity of the milk, the amount of milk added, when it is added, etc. As you say, the mixed Cp is unknown. A naïve initial approach could be to just add the two heat capacities together. If C is the total heat capacity of the coffee, with or without milk, you could add an equation that says:

  C = if time > AddTime then Ccoffee * Vcoffee * 1000 + Cmilk * Vmilk * 1000 else Ccoffee * Vcoffee * 1000;

The coefficients are just there to convert the different units.

The temperature is a bit more difficult, since it varies continuously over time it is state so and can't be as easily changed as the capacity (which changes values only at one discrete interval).

What you have to do with states is using the reinit(var,newValue) function to reinitialize variable var, to the new value newValue. If you mix fluids together, the new temperature is the new total enthalpy divided by the new heat capacity:

t = (m1 c1 t1 + m2 c2 t2 + ... + mn cn tn) / (m1 c1 + m2 c2 + ... + mn cn)

(from Engineering Toolbox)

In Modelica, we could reinitialize the temperature when the simulation time exceeds the time when the milk should be added, using the following:

  when time > AddTime then
    reinit(T, (Ccoffee * Vcoffee * 1000 * T + Cmilk * Vmilk * 1000 * MilkTemperature) / C);
  end when;

I replaced the coffee heat capacitor from experiment 1 in the previously attached models. You could however use the new coffee component in any of the other experiments. This results in a fairly compact model:

enter image description here

If milk is added after 300 seconds, it produces the following simulation:

enter image description here

Approach 2 is having a short flow of milk, instead of a instantaneous addition of milk. The benefit of this is you could create your own addition strategy. For example, you could add half of the milk at the beginning, and half after 300 seconds. Or any arbitrary strategy. For now, I focused on doing it as a 1 second pulse.

An input is added to the coffee, corresponding to the flow of milk:

enter image description here

The volume of the milk in the coffee is no longer a parameter but increases with the flow:

  der(Vmilk) = u;

And the capacity increaces with the milk volume:

  C = Ccoffee * Vcoffee * 1000 + Cmilk * Vmilk * 1000;

Adding milk will increase the enthalpy in the system, but the increased heat capacity will still cause a drop in temperature:

  T = H/C;
  der(H) = port.Q_flow + Cmilk * 1000 * u * MilkTemperature;

With H being the enthalpy.

The milk component is simply a pulse from Modelica.Blocks.Sources.Pulse that has some additional parameters.

enter image description here

Everything taken together, we now have an additional component in the coffee cooling model:

enter image description here

As it should, this approach gives a similar plot as the first one. The only difference is that the milk is added over a duration of 1 second. As the duration approaches zero, the two approaches would converge.

enter image description here

You could use this approach to fit parameters, using the methodology from the electric kettle example. I am a bit unsure if I would recommend trying fitting one of the higher order models that was suggested in the paper. With so many unknowns, you run the risk of significantly over-fitting you model. Though for a student lab, it could be a very good exercise to try and isolate each of the factors and explore them separately. I.e, cover the mug with an isolating material to avoid radiation and evaporation and take some data. Then put it in a open vacuum flask to study the evaporation without the conduction.

Posted 3 years ago

Patrik, your detailed reply is very helpful.

In high school the purpose of the exercise is exploration within limits. These examples provide alternatives which can usefully be explored and verified by direct experiment.

limits of theory and limits of experiment (accuracy and errors) are all important lessons which can be reinforced with open ended exploration.

Comparison of different coffee pouring strategies can easily be measured and now modelled. Different cup materials and configurations can now also be accounted for. surface area, covered/uncovered, material glass, ceramic, metal, colour, double walled, vacuum insulated..

I appreciate your comment about 'over fitting the model' however in physics this is a constant temptation because we can.. Knowing when to stop is an important skill. regards Gary Melbourne (ambient air temp = 32?C today..)

Hi Everyone, I know that I am not actually answering the question, which is already answered, I believe. Here is something that I do in my lectures regarding that same problem. I don't use WSM, but just standard Mathematica and a rather simple approach based on Newton's law for cooling. Please see: Coffee & milk problem with Arduino and Newton's law of cooling

Posted 3 years ago

Nicely done Marco.

Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
or Discard

Group Abstract Group Abstract