CS/Syrthes Coupling

[AdamLarat]

Hi there!

I am trying to reproduce numerically some experimental results : some metal bars are heated by Joule effect and cooled by surrounding air buoyancy.
Thus, I prescribe a heating field within the bars on the Syrthes side and couple the heat transfer simulation with the Code_Saturne solver.
Unfortunately, my numerical setup cools the bars 15 to 30% more than what is expected by the experimental result.

1) Could someone explain me how the Syrthes/Saturne coupling is done ? From what I have read so far in the Saturne Manual, the coupling uses some kind of duality method that I don't grab entirely.

2) At first, if I would have to code a simplified version of the coupling myself, I would impose a heat conduction transfer at the boundary (BC) between the air and the metal as :
F_{BC} = h.(T_{air} - T_{metal}),
with the sign of the transfer coefficient 'h' depending on which side I am looking at.
Is there any such 'h' coefficient in the Syrthes/Saturne coupling I could modify to control the heat conduction efficiency?

Thanks !

[Yvan Fournier]

Hello,

The method used in this coupling (and the coupling with the newer code_saturne external thermal solver based on the CDO numerics) is as follows :

At each time step :
- The solid solver (Syrthes or code_saturne solid domain instance) sends a solid wall temperature.
- The fluid solver (code_saturne) sends a boundary cell fluid temperature and an exchange coefficient.

So the thermal code has a flux based on h.(Tw - Tf)
And the fluid a Dirichlet value based on Tw.

The "h" exchange coefficient is a numerical coefficient, not necessarily a "true" exchange coefficient in the engineering sense, as depending on the mesh refinement, the fluid temperature (which is the mean temperature in a boundary cell) may be closer to the wall temperature or the bulk fluid temperature, so for a fixed/converged flux (at mesh convergence), (Tw-Tf) would diminish and h increase in a corresponding manner to obtain the correct flux.

Compared to simply passing a flux, this involves a relaxation, which is not ideal for unsteady computations, but helps with robustness (for a truly conservative result, using code_saturne's internal coupling instead is better).

The computed exchange coefficient is based on the wall law modeling, so you may naturally get different results depending on the chosen wall law. Also, a "user defined" exchange coefficient may be provided, using the boundary condition code (icodcl) 15 if I am correct, but I do not have any example at hand for how this can be set up.

Best regards,

Yvan

[AdamLarat]

Thanks for your answer, Yvan. It helped a lot!
It appears the overcooling problem comes more from time-stepping than from the way the coupling is done. Nevertheless, I am still interested in understanding the heat transfer coupling, so unfortunately this will double the number of questions I want to address. Sorry for that!

1] Let's start with the thermal coupling : I am interested in the steady state, I am fine with a non-conservative relaxation method, as long as I understand the physical meaning of the equilibrium equations verified at steady state.
1.a) Is there a reference to the relaxation method you are using for the thermal coupling or is it home-made ?
1.b) I quite don't understand the equations verified by Tw and Tf at steady state. If \rho is the volumic heating in the solid part, and f = h.(Tw-Tf) is the heat flux at the solid surface, then at steady state one should get
\int_{SOLID} \rho = \int_{SOLID SURFACE} f (Eq. 1),
the heating produced by Joule effect inside the solid is compensated by the cooling capacity of the surrounding fluid. Since "h" varies with the simulation variables and particularly with Tf, as you claim, it is not clear for me the impact of (Eq. 1) on the temperatures. Or, to put it another way, is "h" designed in such a way that (Eq. 1) is guaranteed at steady state and why?

2] Time-stepping : in this project, I have gone through multiple difficulties, one of them being the stability of the simulation. I already posted on this forum about it (viewtopic.php?p=18018#p18018). During the stabilization of my simulation, I fixed both Saturne and Syrthes time step Δt at a same fixed value, which appears to be way larger than the CFL time-step prescribed by Saturne. Now that my simulation is stable, I have been able to reduce this fixed time step, to simply notice that the overcooling I was observing is just resulting from the numerical error in time : the steady state temperature increases with decreasing Δt, until converging to an estimated value of Tw(t=\inf), which contains only "spatial error". I am happy with this!
But now, I don't quite understand how to properly drive the time stepping of my simulation and make it as efficient as possible.
- If I let Saturne drive the time-stepping by enabling its pseudo-steady algorithm (viewtopic.php?t=3224), "<time_passing>1</time_passing>" in the setup.xml file if I understand quite well, the time-step Δt_{SOLID} prescribed on the Syrthes side does not seem to matter. The solution is the same whatever the value of Δt_{SOLID} is.
- On the contrary, if I fix the fluid time step Δt_{FLUID}, the global time-stepping of the coupled simulation is driven by both
Δt_{FLUID} and Δt_{SOLID}.
2.a) Could you explain me quickly the time marching process of the coupling? Do both code march forward at their own time step and synchronize regularly? It seems to me that the coupling occurs at each time step, but that the time-steps on each side may be different, which results in different "simulation times" in the solid and fluid sides. However, when the pseudo-steady algorithm is chosen on the Saturne side, then the time steps seem to synchronize. Where could I found more information on that specific and meticulous topic?

Once again, sorry for the density of my questions. I am trying to be as clear as possible.

Best regards,

Adam

[Yvan Fournier]

Hello,

Yes, for 1), you should have equilibrium at the steady state. The lower the conductivity h, the higher (Tw-Tf) will need to be to obtain a given flux, the steady solid temperature will be higher if h is lower, lower if h is higher.

For 2), for a steady solution, convergence will often be faster when decoupling the tie steps. So that should be is automatically chosen for the pseudo-time step. For a transient computation, coupling the time steps is needed, so with a fixed or adaptive (non-local) time step, the tie steps can be synchronized.

I am not sure whether this is set automatically, but calling

Code: Select all

cs_coupling_set_sync_flag(PLE_COUPLING_TS_MIN);

in cs_user_parameters.cpp will ensure the smallest time step of coupled codes will be used for example.
Possible enum values are in libple/src/ple_coupling.h.

Best regards,

Yvan

[AdamLarat]

For 1), you told in your first answer that the 'h' coefficient varies during the simulation as a relaxation flux to micmic conservativity. I perfectly understand that the solid temperature at steady state varies with a fixed value for 'h', but I am going to use the "automatic" coupling with the varying 'h'. In this context I am not sure what the final wall temperature physically represents.

For 2), I have noticed myself that the convergence to steady state may be accelerated when decoupling the time steps. But I have also noticed that the "varying time step" marching method on the Saturne side, chooses the most adapted time step and provides a numerical solution with a very low time error.
Unfortunately, when using the "varying time step" in Saturne, the time stepping in Syrthes is automatically set to the same value. Is there any way I could benefit of both sides : an automatic time stepping on the fluid side which is chosen according to the CFL constrain and a large solid time step, which provides a way faster convergence to the steady state?

[Yvan Fournier]

Hello,

1) We assume that for a pseudo-steady computation, h and temperatures will all stabilize, so the flux expressed as "h(Tw-Tf)" should stabilize. Ideally, in a linear wll law region in a mesh convergence study, the flux should reach a fixed value (assuming wall law computation is also ideal), but (Tw-Tf) will be proportional to the mesh refinement dx, so h will be proportional to 1/dx.

The final wall temperature should be the actual wall temperature (solid side) reached at equilibrium, assuming the turbulent conductivity is computed with a good enough accuracy. And that depends on the quality of the wall model.

2) For the varying time step in code_saturne, the synchronization of time steps may be a bug (it was never synced by default in older versions, and is probably synced too aggressively now).

If you add

Code: Select all

cs_coupling_set_sync_flag(PLE_COUPLING_TS_INDEPENDENT);

in the cs_user_parameters function, it should do the trick.

Best regards,

Yvan