[vabishek]

Hello FOAMers,

I recently created a solver for steady-state, incompressible flow, multi-region conjugate heat transfer in a single rotating frame. The equations that are solved in the fluid and solid regions are listed in one of the attached images (equations.png). The variable V is the velocity in the rotating frame and Vrot in the solid equation is the convective velocity due to relative motion in the rotating frame of reference.

A cross-sectional view of the geometry can be found in the attached image. Symmetry is applied to the top and bottom walls and cyclic BCs are applied to the side walls. A volumetric heat source is applied. The rotating frame is fixed to the fluid region. The Solid-2 rotates with the fluid region. Whereas, "solid-1" remains stationary i.e. it rotates in the opposite direction relative to fluid region in the rotating frame. Hence, the convective term in the heat equation for this particular region is non-zero.

I ran a couple of test cases to make sure if everything is in order but I noticed some strange temperature profiles in the circumferential direction in "solid-1" region near the solid-1 - heat zone interface (The dotted line in geometry.png shows the line where the temp. profile was plotted). In the case where the solid-1 region was rotating with the fluid region I saw a perfectly symmetric temperature profile in the circumferential direction, as expected (rotaing_Tprofile.png). However, when the solid-1 region was set to stationary or set to rotate in the opposite direction relative to friction plate in the rotating frame, I am seeing some strange results. Of course, the temp. profile in this case will not be symmetric due to the shear action of the solid-1 region at the solid-1 - fluid interface. However, I would expect the slope at the two ends of cyclic boundaries to be similar. I don't see that (nonrotaing_Tprofile.png). At this point, I am not sure if I missing something or if indeed the slopes won't be same at both ends.

I believe this inconsistency arises from the convective term in the heat equation in solid region. I am not sure if the cyclic boundaries is the root cause or if my implementation of the heat equation with convective term is incorrect. The following code reflects the implementation of heat equation in the solid region:

tmp<fvScalarMatrix> TEqn
(
  fvm::ddt(T)
  + fvm::div(((linearInterpolate(SRFS.U())) & solidRegions[i].Sf()),T) 
  - fvm::laplacian(DT, T)
  == HeatSource
);
TEqn().relax();
TEqn().solve();

Any obvious mistake in the above implementation?

Any suggestions/pointers would really help. Please let me know if you need more information.

Thanks,

Abishek

[wyldckat]

Hi Abishek,

As I mentioned in the answer I sent you to the PM you sent me, without a test case and access to the source code, I'm not able to diagnose problem.

My guess is that the problem is not in the equations, but instead the problem might be in the boundary conditions that relate the static and the rotating regions.
The other possibility is that the flux disturbances caused by mesh motion or by the MRF/SRF are not being taken into account for the heat distribution equation.
Beyond this, my guess is that you should take a closer look into the tutorial "$FOAM_TUTORIALS/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger" in the latest versions of OpenFOAM.

Best regards,
Bruno

[wyldckat]

Hi Abishek,

Even though you've sent me the case via PM, I'll answer most of the questions here in public, without disclosing the case itself.

Originally Posted by vabishek

As I had mentioned earlier, the solver seems to work fine when all the solid bodies rotate at same speed and in the same direction. I have issues with cyclic boundaries when solid is stationary wrt to the fluid region. Both these folders have a Allrun script that you can run. I also attached a sample image for each of these cases in the respective tarballs. It shows the temperature profile in the circumferential direction (X-axis) near the solid1-heatsource interface.

In the case of rotating solid, the temperature profile is perfectly symmetric and the slopes of temp near the wall are similar/continuous.

OK, the images are the same as in the initial post.

Originally Posted by vabishek

However, when the convective term in solid equation is active i.e. stationary_solid case, the slopes of temperature becomes dissimilar/non-continuous (see stationary_solid/solid1_temp_profile_mean_radius.png). I plotted this profile along x-axis at mean radius of geometry at z = 0.0009. So, the starting point is (-0.00505630066618323, 0.0795000046491623, 0.0009) and the end point is (0.00505630066618323, 0.0795000046491623, 0.0009), to be precise. .

So essentially the temperature plot is near one of the solids and the axis used for the plot is aligned with the rotation direction itself.
The first problem I noticed is that the plot axis in the images do not match up to the dimensions you stated, namely "-0.0050563 to 0.0050563" doesn't match up with the plots that have "0.0 to 0.01xx"... unless you transposed the X axis by "+0.0050563".

From what I can figure out, the problem doesn't seem to be in the solver itself, the temperature plot seems to be reflecting the flow being rotated as a result of the Coriolis effect. Since you still defined in the file "constant/S1SRFProperties" that the fluid domain is rotating, even though you defined the other two solids are stationary, but the domain itself is still rotating. In addition, there is nothing stating if the region "HeatSource" is static or not, i.e. "constant/S4SRFProperties" doesn't exist... so I guess it's using the definition that is in "constant/SRFProperties"? Which it is also rotating?

In other words, even though the solids are not rotating, the heat source is rotating, so is the fluid domain, therefore the temperature distribution seems legitimate for the stationary solids, because the flow is being pushed against the heat source on one side and pulled away on the other side.

By the way, the solver wasn't in the shared folder. But it's for the best, because I hopefully diagnosed this correctly after looking at the mesh, how the parts are connected and how the boundaries and SRF settings were defined.

In addition, you might want to double-check if the diagnosis about the flow is correct. Use an already existing SRF solver to test the flow for the fluid part, i.e. to confirm if the flow profile with static solids is correct in your solver.

Best regards,
Bruno

[vabishek]

Hi Bruno,

Thank you for looking into the issue.

So essentially the temperature plot is near one of the solids and the axis used for the plot is aligned with the rotation direction itself.
The first problem I noticed is that the plot axis in the images do not match up to the dimensions you stated, namely "-0.0050563 to 0.0050563" doesn't match up with the plots that have "0.0 to 0.01xx"... unless you transposed the X axis by "+0.0050563".

Yes, that is correct. The x-axis is transposed by a value of "+0.0050563".

From what I can figure out, the problem doesn't seem to be in the solver itself, the temperature plot seems to be reflecting the flow being rotated as a result of the Coriolis effect. Since you still defined in the file "constant/S1SRFProperties" that the fluid domain is rotating, even though you defined the other two solids are stationary, but the domain itself is still rotating. In addition, there is nothing stating if the region "HeatSource" is static or not, i.e. "constant/S4SRFProperties" doesn't exist... so I guess it's using the definition that is in "constant/SRFProperties"? Which it is also rotating?

In the case of stationary solids, the regions "HeatSource and Solid2" rotate with the fluid domain. I account for this in the BC of the fluid region. The SRFVelocity at the walls of Fluid_to_HeatSource and Fluid_to_Solid2 is set to relative "yes". Since, they rotate with the fluid domain the rpm in relative frame is set to zero. However, Solid1 is stationary with respect to the fluid region. Hence, that region rotates in the opposite direction of the fluid domain but at the same speed. There is no separate constant/SRFProperties for the HeatSource region as I am using the same properties from region "Solid2" for HeatSource. This can be seen in the Allrun script where all region properties are set.

In other words, even though the solids are not rotating, the heat source is rotating, so is the fluid domain, therefore the temperature distribution seems legitimate for the stationary solids, because the flow is being pushed against the heat source on one side and pulled away on the other side.

As I mentioned earlier, the temperature profile is plotted in the region "solid1" (closer to the solid1 - HeatSource & Fluid interface) which is stationary. Since, I am applying constant heat source to the right and left of the fluid region, shouldn't I expect similar slopes of the temperature curves as it moves along the direction of rotation. Say, if I move from right to left on the temperature plot, I would expect the temperature to go up in the region near the heat source, temperature would then drop near the fluid region and there will be increase in temperature as it gets near the heat source region, again. I would expect similar behaviour both in the case of rotating and stationary solid with the exception of slopes and symmetric curves. In the case of rotating solids, since all the regions including the fluid domain rotate at the same speed, you would expect a perfectly symmetric curve (that is exactly what we see in the first plot). However, when the solid1 is stationary with respect to fluid region, the symmetric profile is no longer expected due to the fact that the solid is now stationary. However, we would expect a similar slopes of temperatures near the heat source but I we don't see that. Remember that the boundaries on the sides have cyclic BCs, so if I repeat the temperature plot more than once by mirroring them I won't see any continuity of the temperature profiles. Doesn't this mean that there is a discontinuity of temperature profile at the cyclic boundary? That the first derivative is not being computed correctly? Or am I missing something here?

If I understand you correctly, we don't see similar slopes at the boundaries because the flow is being pushed against the heat source on one side and pulled away on the other side?

By the way, the solver wasn't in the shared folder. But it's for the best, because I hopefully diagnosed this correctly after looking at the mesh, how the parts are connected and how the boundaries and SRF settings were defined.

I just uploaded the solver and sent you the link via PM.

In addition, you might want to double-check if the diagnosis about the flow is correct. Use an already existing SRF solver to test the flow for the fluid part, i.e. to confirm if the flow profile with static solids is correct in your solver.

I have already done this simple exercise and everything in terms of fluid flow seems to be represented accurately.


Thanks,

Abishek

[wyldckat]

Hi Abishek,

Attached are 2 images that shows the plot at time "10", because I didn't want to run the complete simulation:

1. "Screenshot from 2015-11-15 20:55:55.png" - this image shows the plot you defined for the "solid1" region, and it also shows 2 marked (selected) points, showing on the right the respective values in a spreadsheet view that ParaView provides. These two points are at the cyclic boundaries. The values shown on the spreadsheet show that the temperature is not identical.
2. "Screenshot from 2015-11-15 20:59:59.png" - this is another image that shows the same exact plot, but this time I also loaded the left and right patches (the cyclic ones) for this region "solid1". As you can see, the temperature values are now identical.

The second image is the accurate one, because the values on the patches have the accurate values on said patches. The first image shows what happens when only the internal mesh is used, which doesn't have the temperature definitions on said patches, it only has the values for the centre of the cells... although in fact, ParaView uses the interpolated values that are on the vertexes of the cells, due to how ParaView interpolates from point data (at the vertexes).


The cyclic boundary condition only ensures that the values on each patch are identical. It does not ensure that the gradients near the patches are identical as well.

I hope this makes things clearer?

Best regards,
Bruno

[vabishek]

Originally Posted by wyldckat

Hi Abishek,

Attached are 2 images that shows the plot at time "10", because I didn't want to run the complete simulation:

1. "Screenshot from 2015-11-15 20:55:55.png" - this image shows the plot you defined for the "solid1" region, and it also shows 2 marked (selected) points, showing on the right the respective values in a spreadsheet view that ParaView provides. These two points are at the cyclic boundaries. The values shown on the spreadsheet show that the temperature is not identical.
2. "Screenshot from 2015-11-15 20:59:59.png" - this is another image that shows the same exact plot, but this time I also loaded the left and right patches (the cyclic ones) for this region "solid1". As you can see, the temperature values are now identical.

The second image is the accurate one, because the values on the patches have the accurate values on said patches. The first image shows what happens when only the internal mesh is used, which doesn't have the temperature definitions on said patches, it only has the values for the centre of the cells... although in fact, ParaView uses the interpolated values that are on the vertexes of the cells, due to how ParaView interpolates from point data (at the vertexes).


The cyclic boundary condition only ensures that the values on each patch are identical. It does not ensure that the gradients near the patches are identical as well.

I hope this makes things clearer?

Best regards,
Bruno

Hello Bruno,

Thank you for your explanation. Things are much clearer now. If I understand correctly, I need to write my own piece of code that would take care of gradients near the patches as part of the cyclic boundary condition? Any pointers? I know it is not that straightforward but I am trying to get my head around to see if something like that is feasible.

Thanks,

Abishek

[wyldckat]

Originally Posted by vabishek

If I understand correctly, I need to write my own piece of code that would take care of gradients near the patches as part of the cyclic boundary condition? Any pointers? I know it is not that straightforward but I am trying to get my head around to see if something like that is feasible.

Quick answer: I can't imagine another boundary condition that would do what you're trying to enforce, at least not from what I can remember seeing in OpenFOAM.

What I forgot to mention is that from what I could figure out, the results seemed physically correct. Have a look at the temperature profile in the fluid region and you'll also see that it's more intensely transferring heat more to one side than the other.
Because if most of the solid regions are stationary and only one solid is rotating, then that means that there is practically a moving wall going over one side of the fluid region, resulting in a swirl effect being induced inside the fluid region.

[vabishek]

Originally Posted by wyldckat

What I forgot to mention is that from what I could figure out, the results seemed physically correct. Have a look at the temperature profile in the fluid region and you'll also see that it's more intensely transferring heat more to one side than the other.
Because if most of the solid regions are stationary and only one solid is rotating, then that means that there is practically a moving wall going over one side of the fluid region, resulting in a swirl effect being induced inside the fluid region.

Yes, that seems to be case. I am carrying out further tests to make sure I fully understand the physics behind this convective phenomena, and that the results make total sense.

Thanks for your help on this, Bruno. I appreciate it.

Regards,

Abishek

[sadsid]

Originally Posted by vabishek

Hello FOAMers,

I recently created a solver for steady-state, incompressible flow, multi-region conjugate heat transfer in a single rotating frame. The equations that are solved in the fluid and solid regions are listed in one of the attached images (equations.png). The variable V is the velocity in the rotating frame and Vrot in the solid equation is the convective velocity due to relative motion in the rotating frame of reference.

A cross-sectional view of the geometry can be found in the attached image. Symmetry is applied to the top and bottom walls and cyclic BCs are applied to the side walls. A volumetric heat source is applied. The rotating frame is fixed to the fluid region. The Solid-2 rotates with the fluid region. Whereas, "solid-1" remains stationary i.e. it rotates in the opposite direction relative to fluid region in the rotating frame. Hence, the convective term in the heat equation for this particular region is non-zero.

I ran a couple of test cases to make sure if everything is in order but I noticed some strange temperature profiles in the circumferential direction in "solid-1" region near the solid-1 - heat zone interface (The dotted line in geometry.png shows the line where the temp. profile was plotted). In the case where the solid-1 region was rotating with the fluid region I saw a perfectly symmetric temperature profile in the circumferential direction, as expected (rotaing_Tprofile.png). However, when the solid-1 region was set to stationary or set to rotate in the opposite direction relative to friction plate in the rotating frame, I am seeing some strange results. Of course, the temp. profile in this case will not be symmetric due to the shear action of the solid-1 region at the solid-1 - fluid interface. However, I would expect the slope at the two ends of cyclic boundaries to be similar. I don't see that (nonrotaing_Tprofile.png). At this point, I am not sure if I missing something or if indeed the slopes won't be same at both ends.

I believe this inconsistency arises from the convective term in the heat equation in solid region. I am not sure if the cyclic boundaries is the root cause or if my implementation of the heat equation with convective term is incorrect. The following code reflects the implementation of heat equation in the solid region:

tmp<fvScalarMatrix> TEqn
(
  fvm::ddt(T)
  + fvm::div(((linearInterpolate(SRFS.U())) & solidRegions[i].Sf()),T) 
  - fvm::laplacian(DT, T)
  == HeatSource
);
TEqn().relax();
TEqn().solve();

Any obvious mistake in the above implementation?

Any suggestions/pointers would really help. Please let me know if you need more information.

Thanks,

Abishek

How you calculated Vrot in your solver?