SUGAR BEET CLAMP CFD MODELLING: HEAT OF RESPIRATION AND CONVECTION ONLY

OpenFOAM version: 9

Case: clamp_11

Solver: clampPimpleFoam v3

BACKGROUND

So far, the clamp modelling project has developed to the point where basically all the components of the model are in place. However, and as noted in the last post in the Results section, issues remain. The rate of heat transfer in the porous zone is off – the air very quickly cools to the temperature of the beets, but the beets do not seem to change temperature between 60 seconds and 60 minutes despite a large difference between the air and clamp temperatures. This is unexpected given that old experiments run by NBR showed that it took 6 hours for the core of a beet to cool from 12 degrees to -2 degrees when the air was -3. I’d expect that a 300K – 274K = 26 degree difference between initial temperatures of the beet mass and air would results in greater temperature changes over this period.

To see how the model can be improved, the entire system is being reviewed. Also, the experts at the LTH Department of Energy Sciences have been consulted.

OTHER FIXES

To start with, a few minor and not-so-minor changes have been identified.

• Initial conditions on U and T is better if it reflects the fixed inlet condition
• p was behaving badly in early iterations, so the nOrthagonalCorrectors was increased from 0 to 1. nCorrectors was also increased from 2 to 3. This seems to have fixed that problem.
• Porosity was not included in the source heat transfer term (between solid and fluid) in the Ts and Tf equations. This was corrected by dividing these components by the porosity (por) in Tf, and 1-porosity (invPor) in Ts. The new scalar invPor was added as a lazy/ amateur way of avoiding the issues that dividing directly by 1-por gives in the areas outside the clamp where por = 1; that is, it avoids the error that is thrown when one divides by zero. By making invPor a constant, and multiplying this constant by clampDummy = 0 in the open air space, the right result is gained, and the error avoided. The problem with this lazy programming is that it introduces the risk that a future change in porosity that is implemented in the setFieldsDict is not implemented in the transportProperties file.

HEAT OF RESPIRATION ONLY

Running the model with heat of respiration only was a way of testing if this part was working properly and not the cause of the above noted issue (eg, was Qr producing “cold” instead of heat?). The end result was that this wasn’t the source of the issue, but Qr needed refinement. Also, some other future improvements to the model would be desirable (see end of this post). The value was Ks was revised (from 0.53 to 0.554) after revisiting the work of Tabil (2003).

FIRST RUN

To get the model to run with heat of respiration only, Ks and Asf were set to 0 – there is not convection and no transfer between phases. With inlet U = 0 0 0, the model runs super quickly, which is nice. But it runs quickly to infinity. That is, with Qr as an exponential model, once the clamp starts to warm up, the rate of respiration grows, which produces more heat, and so the viscous cycle we often talk about is entered. In the end, the clamp becomes a supernova. Fun, but probably not something you want to be remembered for – causing the destruction of the earth, if only in a computer simulation. So there were two courses for action – revise the rate at which thermal energy is ‘created’, or revise the parameters that move heat out of the clamp. Both were done.

REVISING Qr

Qr is completely defined within the solver and as such, each update needs a recompilation (run wmake) of the solver.

Qr was originally sourced from Shaaban (2014) – a conference paper written as part of a MSc program, and was given as exp(25.292 – 6291(1/T)). It was this value that was used in the FIRST RUN above. When revisiting this work, a new version of Qr was found in the 2020 MSc thesis: exp (25.472 – 6291(1/T)) / 182. The “/ 182” makes quite a difference!

Before this new equation was found, it seemed relevant to change the functional form of the model. Exponential growth in respiration is both dangerous in the modelling (creates supernovas), and goes against what is known about plants; there is a maximum rate they can metabolise, and this is often reached around 20 to 25 degrees. So, a Logit expression was derived that matched the curve of Shaaban well between 0 and 15 degrees (our likely operating temperatures), and had a maximum equivalent to that of Shaaban at 20 degrees. RMSE was used to guide the setting of parameters. The final model was Qr = 45.650 / (1 + exp(-0.145*(Ts-283))).

Using the logit version of Qr, and with 275 K initial temperatures, the clamp region went up a little over 0.1 degree over 12 hours. Historical data shows an increase of 0.5 degree over 24 hours with a starting temp of 1 degree, and that’s with an ambient temp of around 0 C and probably with some transfer out of the system. Also, my calculations suggest around 600 W/m3 are needed to gain 1 C, but Qr was only going at 23 W. These both suggest a multiplication factor is needed (not a change to the new Qr equation). In the end, a multiplication factor of 4 worked quite well. This needs further testing, makes me wonder where the ‘/ 182’ came from, and shows that ‘my calculations’ are not very good.

Taking the ‘/182’ version of Qr, the Logit equivalent that minimses RMSE between 0 and 18 degrees with the max set to that of 25 degrees is = 0.431 / (1 + exp(-0.116*(Ts-289))). This wasn’t added to the solver.

HEAT OF RESPIRATION AND CONVECTION

When the convection coefficient – Ks – was returned to the original value of 0.554, the model looked to behave well.

However, it became clear that the model does have the boundary issues Morteza had flagged once or twice as something to look out for. In the above video, a gradient can be seen at the edge of the clamp, but really this shouldn’t exist when only respiration and conduction within the solid phase is what is included in the model. The issue seems from the fact that the Ts field covers the whole domain, but doesn’t currently have any connections to the rest of the model expect through conduction with the clamp region at its boundary – it basically stays the initial value. One way of fixing this issue is to break this connection between Ts in the clamp region and Ts not in the clamp region, but I’m not yet sure how to do that. An alternative would be to use this to advantage. Given it is unlikely that the temperature difference between Tf (air) and Ts (beet) at the edge of the clamp will ever be high, setting Ts in the region outside the clamp region to Tf would help avoid this issue, and would speed up heat transfer at the clamp edge. It is not certain that this is actually needed. It is certain that more knowledge is needed to code this. Also, it is certain that the pros at LTH are pros.

One of these pros suggested that the Ts – Tf used in solver for the heat transfer between phases should be made implicit (maybe semi implicit) instead of explicit like it currently is. That is, it shouldn’t rely on values from the previous time period as that would increase the risk of error following through with each iteration. This issue doesn’t seem to be root of all issues. When the sections of the solver that include Ts – Tf are included again (Asf > 0), then the results seem reasonable while U = 0. At the same time, if Asf = 30.04 (the modelled value for specific surface area), the simulation blows up – too much energy is being transferred in each iteration.

JUST A LITTLE VELOCITY

And it all falls apart. When the inlet velocity is set to anything above 0 m/s, the simulation stalls. In a 12 hour simulation, dTs/dt starts well for the first couple of hours but then stops progressing again (which it definitely shouldn’t do, and doesn’t do when U = 0).

Before working out how this can happen, the first step is to make Ts and Tf implicit…

IMPROVEMENTS FOR THE NEXT MODELLING PROJECT

Need to improve: soil, water, porosity along the ground? I did have some good ones, but they’ve slipped away for now…