conservative versus non conservative form heat advection - strange boundary effects
81 views
Skip to first unread message
mauro13
Mar 31, 2016, 10:55:36 AM3/31/16
to moose-users
Dear moose users,
I've been encountering some strange behaviour while solving a coupled pressure and heat transfer problem (darcy flow plus advection and diffusion, no variable fluid density).
In my problem, I have two blocks (porous units) which are cut by two intersecting faults (here represented by lower dimensional surface elements, i.e. the traces are depicted in the figures).
I'm advecting a non-isothermal from from the back to the front (from 40 to 80 degrees) via a constant pressure gradient imposed on the respective sides. The front thermal boundary is not fixed.
While everything works fine when I do solve the coupling via a non conservative form for the advective kernel, I do observe some strange boundary effects if relying on a fully conservative form (negative temperatures along the front boundaries propgating within the model).
I followed the basic outlines of some previous posts (https://groups.google.com/forum/#!topic/moose-users/yWi8bdfZWvU), and did calculate the additional residual and jacobian coming from the advective energy term.
This indeed gave a better result. Now the temperature range is OK, the results at 3 OPs agree (see the pdf's).
However, a strange boundary behaviour remains localized around the boundary still (i.e. *.png's).
Since it is also mesh size dependent (the bc is an integrated one), I can try to minimize by imposing local refinement though without eliminating it from the simulation.
I'm wondering whether some of you already encountered the same problem and if you arrived at a better (if at all) solution.
Thanks in advance, and sorry for the long post.
I've been encountering some strange behaviour while solving a coupled pressure and heat transfer problem (darcy flow plus advection and diffusion, no variable fluid density).
In my problem, I have two blocks (porous units) which are cut by two intersecting faults (here represented by lower dimensional surface elements, i.e. the traces are depicted in the figures).
I'm advecting a non-isothermal from from the back to the front (from 40 to 80 degrees) via a constant pressure gradient imposed on the respective sides. The front thermal boundary is not fixed.
While everything works fine when I do solve the coupling via a non conservative form for the advective kernel, I do observe some strange boundary effects if relying on a fully conservative form (negative temperatures along the front boundaries propgating within the model).
I followed the basic outlines of some previous posts (https://groups.google.com/forum/#!topic/moose-users/yWi8bdfZWvU), and did calculate the additional residual and jacobian coming from the advective energy term.
This indeed gave a better result. Now the temperature range is OK, the results at 3 OPs agree (see the pdf's).
However, a strange boundary behaviour remains localized around the boundary still (i.e. *.png's).
Since it is also mesh size dependent (the bc is an integrated one), I can try to minimize by imposing local refinement though without eliminating it from the simulation.
I'm wondering whether some of you already encountered the same problem and if you arrived at a better (if at all) solution.
Thanks in advance, and sorry for the long post.
Peterson, JW
Mar 31, 2016, 11:09:56 AM3/31/16
to moose-users
Interesting images!
I didn't understand all that you wrote -- we might need some LaTeX/equation PDFs to understand exactly the governing equation and boundary conditions that you are solving. For example, I'm not sure exactly what is meant by "non-conservative form for the advection kernel" in this context.
Also, the thread you linked is about boundary conditions for the incompressible Navier-Stokes equations where the pressure is a very different beast, we need to be a little careful about applying those comments to different sets of equations. Finally, if you can share your code and input files that would be very helpful as well, as there may be a bug or something inconsistent with the formulation that we might see from that.
John
mauro13
Mar 31, 2016, 11:39:47 AM3/31/16
to moose-users
Dear John,
Thanks for your prompt answer.
Indeed, I refer to a Navier-Stokes post, but I convert (in my head) to a laminar darcy flow problem.
In the conservative form, I do solve for:
while in the non conservative case I do solve for:
In this second case, I split the energy term (I do need to call two kernels), while I do treat the energy into one single kernel in the first case. I guess this might cause the problem, since the first situation might introduce an additional condition on grad_test, that I tried to include by calculating the additional residual from the boundary term, which I do not need in the second case.
More on the set up: the problem is a benchmark I'm trying to set up in order to test our moose implementation. Indeed, the results look fine when compare to other software (see png when I do compare the result against OpengeoSys).
The main problem, is that (my) non conservative form will fail for fluid density driven problems (i tested against the Elder benchmark) because of the assumption on the advective term.
I would be more than happy to share the code (though bit messy at this stage) and input files. In this regard let me know what is the best platform for sharing.
Thanks alot again for your help and support,
Mauro
Thanks for your prompt answer.
Indeed, I refer to a Navier-Stokes post, but I convert (in my head) to a laminar darcy flow problem.
In the conservative form, I do solve for:
while in the non conservative case I do solve for:
In this second case, I split the energy term (I do need to call two kernels), while I do treat the energy into one single kernel in the first case. I guess this might cause the problem, since the first situation might introduce an additional condition on grad_test, that I tried to include by calculating the additional residual from the boundary term, which I do not need in the second case.
More on the set up: the problem is a benchmark I'm trying to set up in order to test our moose implementation. Indeed, the results look fine when compare to other software (see png when I do compare the result against OpengeoSys).
The main problem, is that (my) non conservative form will fail for fluid density driven problems (i tested against the Elder benchmark) because of the assumption on the advective term.
I would be more than happy to share the code (though bit messy at this stage) and input files. In this regard let me know what is the best platform for sharing.
Thanks alot again for your help and support,
Mauro
Sebastian Schunert
Mar 31, 2016, 12:49:59 PM3/31/16
to moose...@googlegroups.com
When you apply the divergence to the first term you only apply it to T, but qf depends on space, too. Or am I off here?
--
You received this message because you are subscribed to the Google Groups "moose-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to moose-users...@googlegroups.com.
Visit this group at https://groups.google.com/group/moose-users.
To view this discussion on the web visit https://groups.google.com/d/msgid/moose-users/069df515-9667-41b6-998f-f47f50668e46%40googlegroups.com.
Peterson, JW
Mar 31, 2016, 1:03:03 PM3/31/16
to moose-users
On Thu, Mar 31, 2016 at 10:49 AM, 'Sebastian Schunert' via moose-users <moose...@googlegroups.com> wrote:
When you apply the divergence to the first term you only apply it to T, but qf depends on space, too. Or am I off here?
I agree, it appears that the "non-conservative" form has dropped a
T * div( \rho_f c_f \vec{q}_f) term
This is sometimes dropped because it is zero by the mass conservation equation, i.e. div( \rho_f c_f \vec{q}_f)=0, but you are not using that form of the mass conservation equation.
John
mauro13
Apr 1, 2016, 3:03:29 AM4/1/16
to moose-users
Hi to all,
and thanks for your replies.
Indeed the two forms are equivalent only in a divergent-free case (which is not my case).
To correct for the mass (pressure) term in the non conservative case, I do add off diagonal entries in the jacobian while solving for the heat advective part: T * \rho_f \c_f * ( k / visc) * grad_p_f (in this particular case densities are constant).
The missing of this term in the physics, I guess is the main reason while this version does not work well for high density dependent problems like the Elder benchmark.
mauro
and thanks for your replies.
Indeed the two forms are equivalent only in a divergent-free case (which is not my case).
To correct for the mass (pressure) term in the non conservative case, I do add off diagonal entries in the jacobian while solving for the heat advective part: T * \rho_f \c_f * ( k / visc) * grad_p_f (in this particular case densities are constant).
The missing of this term in the physics, I guess is the main reason while this version does not work well for high density dependent problems like the Elder benchmark.
mauro
mauro13
Apr 1, 2016, 7:12:04 AM4/1/16
to moose-users
I might have an additional question regarding best way (if there's any) to calculate the pressure gradient term.
I tried two different options (giving the same results): one by calculating in the material class and the other by calculating in the kernels (where needed). Was wondering which one should be preferred.
Thanks
I tried two different options (giving the same results): one by calculating in the material class and the other by calculating in the kernels (where needed). Was wondering which one should be preferred.
Thanks
Peterson, JW
Apr 1, 2016, 11:42:36 AM4/1/16
to moose-users
On Fri, Apr 1, 2016 at 1:03 AM, mauro13 <mauroc...@gmail.com> wrote:
Hi to all,
and thanks for your replies.
Indeed the two forms are equivalent only in a divergent-free case (which is not my case).
To correct for the mass (pressure) term in the non conservative case, I do add off diagonal entries in the jacobian while solving for the heat advective part: T * \rho_f \c_f * ( k / visc) * grad_p_f (in this particular case densities are constant).
In the Jacobian only? That doesn't seem right.
The missing of this term in the physics, I guess is the main reason while this version does not work well for high density dependent problems like the Elder benchmark.
Then you should not expect a consistent result between the conserved and non-conserved problems, even for smooth solutions. You are simply solving a different problem in the two cases.
But I think the main difference between the two cases most likely has to do with your boundary condition for the conserved form, and in particular its implementation in MOOSE. You had asked previously what the best way to share your code would be, and the answer is GitHub. If you started your app by forking Stork, just push your code to your Stork fork in a branch. Otherwise, finding some other way to push it to GitHub is the best way for us to be able to take a look.
John
Peterson, JW
Apr 1, 2016, 11:51:36 AM4/1/16
to moose-users
On Fri, Apr 1, 2016 at 5:12 AM, mauro13 <mauroc...@gmail.com> wrote:
I might have an additional question regarding best way (if there's any) to calculate the pressure gradient term.
I tried two different options (giving the same results): one by calculating in the material class and the other by calculating in the kernels (where needed). Was wondering which one should be preferred.
From what I understand of your equations, the pressure is a primary variable, so you get access to its gradient by calling coupledGradient() and storing a const reference. See, for example $MOOSE_DIR/test/src/kernels/FDAdvection.C
John
mauro13
Apr 1, 2016, 12:24:39 PM4/1/16
to moose-users
Dear John,
Thanks for your support.
Thanks for your support.
Well, I do not see where else to add (Residual?!)
From what I understand of your equations, the pressure is a primary variable, so you get access to its gradient by calling coupledGradient() and storing a const reference. See, for example $MOOSE_DIR/test/src/kernels/
FDAdvection.C
Indeed, I do calculate the pressure gradient in the darcy eq. by coupledGradient, but that I can do either in the material class (and then passing the full darcy term as a material property to the kernel), or as you suggest in the Kernel. I tried both, but I think you suggest the latter one is the best option.
About the github: I did not fork stork, but one of my colleague did. Unfortunately he's home now. I will wait for him to come back on Monday and use his branch to share the code.
Thanks again,
Mauro
About the github: I did not fork stork, but one of my colleague did. Unfortunately he's home now. I will wait for him to come back on Monday and use his branch to share the code.
Thanks again,
Mauro
Peterson, JW
Apr 1, 2016, 12:48:46 PM4/1/16
to moose-users
In the Jacobian only? That doesn't seem right
Well, I do not see where else to add (Residual?!)
Yes... the residual and Jacobian must always be consistent. If you add a term in one, you need to add a corresponding term to the other.
From what I understand of your equations, the pressure is a primary variable, so you get access to its gradient by calling coupledGradient() and storing a const reference. See, for example $MOOSE_DIR/test/src/kernels/FDAdvection.C
Indeed, I do calculate the pressure gradient in the darcy eq. by coupledGradient, but that I can do either in the material class (and then passing the full darcy term as a material property to the kernel), or as you suggest in the Kernel. I tried both, but I think you suggest the latter one is the best option.
Shouldn't make any difference which way you do it.
About the github: I did not fork stork, but one of my colleague did. Unfortunately he's home now. I will wait for him to come back on Monday and use his branch to share the code.
Even better, create your own account on GitHub and push to that. You can re-fork stork or just start a new repo, doesn't make much difference to me.
John
mauro13
Apr 4, 2016, 8:09:52 AM4/4/16
to moose-users
Dear John,
As discussed we moved our branch on github : https://github.com/ajacquey/fennec/tree/devel
The file you will need to run the simple benchmark I was discussing you will find in examples/2faults2units (there are a lot of unused files as well, sorry).
The implementation of the two forms for the energy component are in the class (H)THMaterial classes and (H)/THKernels.
Thanks a lot for your help and feedbacks.
Mauro
As discussed we moved our branch on github : https://github.com/ajacquey/fennec/tree/devel
The file you will need to run the simple benchmark I was discussing you will find in examples/2faults2units (there are a lot of unused files as well, sorry).
The implementation of the two forms for the energy component are in the class (H)THMaterial classes and (H)/THKernels.
Thanks a lot for your help and feedbacks.
Mauro
Reply all
Reply to author
Forward
0 new messages