Add a source term to simpleFoam
Hello All,
I want to add a force vectorField as body force to the equations in simpleFoam. These forces do not change in time, it will simply be a vectorField with most cells set at zero and some will have a non-zero force vector. Is it just as simple as putting: >>> eqnResidual = solve ( UEqn() == -fvc::grad(p) (new part: + source vector) ).initialResidual(); <<< Probably there is more involved. I searched this forum and surprisingly, similar questions were posted a couple of times but none of them were really answered. Hopefully this time some expert passes by? Brgds, Mark |
it is as simple as that ;)
you just have to make sure that the units match. the following lines are just because I may be tired, so if you find those advice to direct, please forgive me ;) lest say you want to read this field from "constant" folder. in createFields.H you just add: Code:
volVectorField my_field_name Code:
eqnResidual = solve you can do it in the code Code:
eqnResidual = solve hope it helps you a bit. Pawel |
Hi Mark and Pawel
I do not agree that it is as simple as that, since if "momentumPrediction" is set to "no", then the momentum equation is not solved, and hence the body force is not taken into consideration. Thus my experience is that in the pEqn.H one need to modify the momentum corrector as follows: U -= (fvc::grad(p) - my_field / rho )/AU; However, you might need to consider how this contribution is added. I have used a constant in space, varying in time body force, thus the divergence of the momentum equation gives a 0 (zero) contribution to the pressure equation. As your bodyforce in varying in space, the term will take all the way to the pressure equation, I am not certain that it is sufficient to do as above. Best regards, Niels |
Hello Pawel and Niels,
First of all thanks for your contributions. Niels, you explain about momentumcorrector but we are talking about simpleFoam, which always uses a momentum corrector, right? In that case, does it still make sense? Next: now I modified the simpleFoam solver as described by both of you. I tried it on the pitzDaily example (a backward facing step) and add some force on the upperwall, directing downwards. The liquid is not moving downward at all. While writing this I realise this might be caused by the fact that the fluid is incompressible. Does exerting a force on a wall makes sense? I better take some time to test it with a bodyForce on the internal fluid (what is what I actually want to do anyway). Finally dimensions: Pawel, you are right. I need to supply force per unit volume and divided by rho, which gives: F/V/rho = (kg m s-2) * (m-3) * (kg m-3)^-1 = m/s2. Once again thanks for your help, Mark |
Hello Niles and Mark,
Niles, I assumed that "momentumPredictor" is set to "true". Now I have a question- if we assume that we solve the momentum equation, and that we add, for example, a gravity mass force, we get the UEqn like: Code:
tmp<fvVectorMatrix> UEqn Code:
eqnResidual = solve Code:
volScalarField AU = UEqn().A(); Mark, first- check the boundary type- it should be "fixedValue" (I found myself forgeting about it several times) As for the second thing- I belive that by adding this additional field you simply introduse a velocity source term to the UEqn. For me it seems that by putting some value on the wall would be just like applying a velocity inlet on that wall. (But it would be better if someone veryfied my statement). Regards, Pawel |
Hi Pawel
In the pressure equation, U is estimated as U = UEqn().H()/AU; hence based on all the terms in the UEqn, thus not on neither grad(p) nor my_fields. This means that the pressure correction to U does not include any of these, and therefore you need to have it in the second term as well. One of my colleagues has been calculating laminar and turbulent oscillating wave boundary layers using OpenFOAM and he achieve perfect comparison with experimental/theoretical data using the method I have described above. Best regards, Niels |
Thenks Niles for pointing that out.
Fortunately I added my source terms into velocity equation so they were included while pressure calculations. Now I will pay more attention to that matter. Regards, Pawel |
Ah, that explains why you didn't need to add it in the bottom.
Good luck to both of you:) and enjoy the up-coming weekend, Niels |
Hi,
I am also trying to bring in a source term in my simpleFoam solver, following the code snipped below. Quote:
But unfortunately the compiler complains that rho is not defined. So what do I have to do to get thr rho field? Thanks and best regards, Nick |
Hello Nick!
The way I see it, there are two ways of including rho in your case: 1) when you create your field, divide it by rho at the beginning, set the dimensions as [myField/rho] and solve just: Code:
UEqn() == -fvc::grad(p) + myFieldDivRho Code:
volScalarField rho Hope it helped a bit ;) best regards, Pawel |
All times are GMT -4. The time now is 17:20. |