BUG: Coupled simulation of single pin heat sink is unstable

[ChN]

Dear code-saturne community,

I would like to validate the result of code-saturne 3.0.4 coupled with syrthes 4.1.1.
Unfortunately, this seems to be a very tough job, because any simulation I ran just went bananas.
Right now, I don't know how to get a steady state simulation of a natural convection for a single pin heat sink.
The base plate has 50°C and the walls of the fluid are supposed to be fixed at 0°C.
The simulation starts off with sound values, but after some time, it gets unstable ;( -> SYRTHES' mnx-file.
I have tried different types of algorithms: laminar (due to the fact that natural convection would produce small velocities in air, IMHO, this should be the right choice), k-epsilon, k-omega, ... none of those could do the job...

I hope, someone can tell me, what I am doing wrong here...

I have the files stored in this dropbox folder:
https://www.dropbox.com/sh/75gy3kzo70dh ... ln44a?dl=0

Best regards,
Dr. Christoph Neumaier

[Brian Angel]

Hello,

I will take a look and come back to you.

Regards,

Brian Angel.

[Brian Angel]

Hello,

I've run your set up with separate calculations for Code_Saturne and SYRTHES to see if I can get the codes running separately.

SYRTHES:
I set the BaseFin condition to 50°C and the CouplingFinfaces to 0°C, both for heat exchange type BCs with h = 1W/m2K. SYRTHES runs for 10000 iterations with a very small amount of heat transfer into the solid but there are no stability problems.

Code_Saurne:
I set the coupling BC to 50°C, the external walls to 0°C and intialised the field at 0°C. I also set the clipping of the temperature scalar to a minimal value of 0°C and a maximal value of 50°C (I see that your set up uses the default values of +/-1e12). The simulation is running and has done some 300 iterations with numerous clips at each iteration for the minimum temperature which is negative in some 5 or so cells. I suspect that this might be where the problem is coming from.

Can you set your clipping of the temperature scalar to the above min and max values, rerun your calculation and let me know what happens?

Just to note that I'm currently running v3.2.1 of CS which is different to you.

Regards,

Brian Angel.

[Brian Angel]

Hello,

Following on from my previous post, I left the Code_Saturne calculation running and the flow variables have obtained numbers of the order of e+11 after some 600 time steps. I noticed during this calculation that the maximum Courant number increases from one time step to the next, even when the flow variables have reasonable values.

So, I have just started a new calculation and reduced the time step from dt = 1.0s to dt = 0.1s. The solution has run for 2000 time steps and the max Courant number is stable at around 5.1. I have compared the Courant number at 200 time steps for dt = 1.0s and 2000 time steps for dt = 0.1s. For the former, the max Courant number = 8467 whilst for the latter the max Courant number = 5.1. The values of the flow variables also show large differences between the two solutions at 200 and 2000 time steps, with the latter being more stable.

So, can you reduce the time step to 0.1s as well as set the temperature clipping limits and let me know what happens.

Regards,

Brian Angel.

[ChN]

Dear Brian,

thank you for your first efforts! I will try what you have found out and I will return ASAP.

BTW: What also seems weird to me is that in the coupled simulation the velocity at the boundaries of the fluid is pointing right at the surface in direction of the gravity vector (at the wall I would have expected to have only vectors parallel to the wall). Its magnitude is somewhat larger than that inside at the pin, see attached picture ;(... But let's advance step by step!

Best regards,
Christoph

[ChN]

Dear Brian,

First of all: Thanks for the help in the case of stability. This seems to be solved now with your time-stepping analysis in the separate domains. -> Is there a general rule, how high the CourantNb should be allowed to rise? Until recently, I thought it must stay below 1. But according to your answer, it is allowed much higher than 1 (but stable).
I found out that I could also use the adaptive time-stepping in SYRTHES which would keep this simulation part from injecting too high rises in temperature over time into the code-saturne part (preventing instabilities)...
Does anyone know, how these two types time-stepping are synchronized?

BUT: Nevertheless: I am still getting these weird velocity vectors at the boundary walls... The velocity should be PARALLEL, and not perpendicular...
I already tried the following: using an inlet and force the velocity at the bound to be 0m/s in the perpendicular direction and in some other arbitrary angle. Neither of the two helped.
I found out that I could also use the SIMPLE algorithm in the fluid-domain but this has not been the issue.
Then I changed the direction of the gravity (pointing in a direction with z and y component and a length of the usual 9.81m/s^2). This lead to the change in direction of the "plume", but even here: The velocity at the boundary is WRONG. I do not have walls which were supposed to be leaky...
Thus, if I use the streamline-tool in paraview, it shows that not only right at the boundary the velocity has a vector pointing in the wrong direction, but also the interpolation between the realm with the main convection cell is taking place and the boundary is wrong (there is no boundary layer in the result / no crust which I could snip off in the post-processing in paraview). -> boundary condition at the wall is buggy!

Best regards,
Christoph

[Brian Angel]

Hello Christoph,

The CFD codes that I have used to date, both open source and commercial, can remain stable when the Courant number is greater than 1 and still give good results. A general rule of thumb that I use is an average Courant number of the order of 5 but with a maximum CFL which can be greater than this. However, as per usual, there are always exceptions to this.

As far as I know, the time stepping between Code_Saturne and SYRTHES is not synchronised as the solution of these two codes is sequential, but someone with more experience might be able to give further details or correct me on this.

Concerning the near wall velocity vector problem, I haven't come across this one before. Is it not just a question of visualisation and the Paraview process for plotting vectors? The process I currently use in Paraview to plot vectors (in the volume or on a section) is as follows:

1) Convert cell data to point data from the imported results file
2) Extract the fluid domain from the point data
3) Plot the vectors using glyph

Hence when you plot vectors, it would seem that one is only plotting the fluid domain data (no boundary data) at the nodes of the mesh which might explain why you see vectors going into the wall boundary.

Another line of thought is how refined is your mesh in the near wall regions?

Regards,

Brian Angel.

[ChN]

Hello Brian,

I made some tests on this topic....
Although I added the viscous layer hypothesis in my mesh (FluidMeshWithViscousLayer.png), I still get the velocity-error at the bottom (fluid_Velocity_Error.png).
This may also be seen in the Streamlines (fluid_Velocity_ErrorStreamlines.png)

BTW: In paraview, it is only possible to get the velocity vectors iff you decompose the data into the fluid block and then you may recalculate the values from the cell-center to the nodes (what is much nicer, because then data is interpolated).
In order to get rid of the block selection, one may chose in code-saturne to leave out the boundaries and the intermediate layer to syrthes. -> only one block remains which is unique -> paraview knows what to do.

[ChN]

Here are my project files so you may rebuild the error...
I think, since the last time, I changed to the "SIMPLE" solver instead of the "SIMPLEC" solver, but just because it seemed to be producing more reliant data. Can anyone tell me the difference between these two?


Moreover: If I make the timestepping to small in SYRTHES, syrthes4ensight can nomore handle the amount of datapoints... (10000 points or so.) -> This seems to be another bug. Dear SYRTHES guys: Please fix this by adding >>list<< structures instead of >>array<< structures... (IMHO, I assume this could be the cause)

Sorry to uncover one bug after another, but I thought I had a simple example ...

I hope you can help me go on with this problem.

In order to optimize my structure later on:
Is there a possibility to stop the simulation from user defined functions (SYRHTES or code-saturne), when I am in steady state (for example: wait for 10secs in physical time to start the system and then stop the simulation if in the last 5 simulation steps temperature changes less than 0.2K and velocity vectors have a distribution of less than 2%)?

Thank you in advance,
Christoph

[Brian Angel]

Hello,

I will download your files and take a look next week and come back to you.

Regards,

Brian Angel.