[sturgeon]

I'm trying to add a force to a modified simpleFoam. I managed to get the code to execute without any floating point errors but the solution is odd and non-physical so I think I may have made a mistake - in particular, I notice odd effects across regions where the cell size changes, making me think I've mistakenly multiplied by volume/area somewhere, or missed this out.

My implementation is in the UEqn is:

volScalarField u("u", (vector(1,0,0) & U));
volScalarField v("v", (vector(0,1,0) & U));
volScalarField w("w", (vector(0,0,1) & U));
volScalarField Fx("Fx",(f*v-l*w*J));
volScalarField Fy("Fy",-f*u);
volScalarField Fz("Fz",(pow(J,2)-1)*(l*u/J+beta*(T-TRef)*mag(g))+l*u/J+J*zeta*(p-pow(w,2)));

volVectorField F("F",Fx*vector(1,0,0)+Fy*vector(0,1,0)+Fz*vector(0,0,1));

surfaceScalarField Ff = fvc::interpolate(F)&mesh.Sf();

        solve(UEqn == fvc::reconstruct
   (
                (
                     fvc::interpolate(rhok-1)*(g & mesh.Sf())
					+ Ff
					 - fvc::snGrad(p)*mesh.magSf()
                )
				) 
				);

        fvOptions.correct(U);
    }

I am fairly certain my addition of the force F is wrong, but every other variation refuses to compile. I feel like if it's set up properly, I should just be able to have "fvc::interpolate(F)&mesh.Sf();" instead of Ff within UEqn (see this thread: Minor issue compiling solver).

My understanding of "fvc::interpolate(rhok-1)*(g & mesh.Sf())" is that rhok is extrapolated to the cell surfaces, then multiplied by g crossed with the mesh surface, to give the rho*g buoyancy term (ignore the -1 in my code)

g is a uniformDimensionedVectorField, whereas fvc::interpolate(F) returns a volumeVectorField, so I don't think "fvc::interpolate(F)&mesh.Sf();" is doing the same thing as "(g & mesh.Sf())". But I can't get any other variation to compile. I thought rewriting F as a vectorField rather than a volumeVectorField would fix it (so the code is just surfaceScalarField Ff = F&mesh.Sf()) but I can't get it to compile like that, which I find quite confusing since I thought g in this context was a vectorField, unless there's a subtle difference between vectorField and uniformDimensionedVectorField that I am missing?

Anyway, I feel like I am missing something obvious, can anyone help me out?

Thank you

EDIT:

Looking up uniformDimensionedVectorField, the only reference I can see is related to the g term, and if I add "Info << g << nl" in my code, the console shows it's just a dimensioned vector, despite the field name?

This would make sense because I can do something like "Ff = F[1]&mesh.Sf()" and it compiles fine, so I need a single vector, not a vector field. I am thinking perhaps g&mesh.Sf() is the cross product of the vector g to every vector within mesh.Sf(), and I need to do the same thing but element by element, since F is a field rather than a singular vector? But this makes me think my original method of "surfaceScalarField Ff = fvc::interpolate(F)&mesh.Sf();" was correct. (Personally I am hoping not since I can't find any other areas where I could have made a mistake) Does anyone have any insight? Cheers

[sturgeon]

Sorry for the bump but can anyone confirm if my method is correct or not? Thank you

[TomasDenk]

Hi Sturgeon,

I'm not able to pin point where in your code is the mistake, because there are several symbols that are not defined in the code snippet. What I would do first if I were in your situation is dimension analysis. The equation for U is momentum equation and its dimension is N*m. So it seems that all components of your vector F should be in N/m. Make sure that all terms in Fz (which is sort of complex) are od the same dimension too.

I would also warn you not to use practice like

beta*(T-Tref)*mag(g)

very often. As long as g is defined as vector, I would prefer using it in a general way. Imagine a user that decides to define

g = (0, -9.81, 0) m/s^2

One term in your equation respects the choice

fvc::interpolate(rhok-1)*(g & mesh.Sf())

the other does not

beta*(T-Tref)*mag(g)

This is very confusing.

A while ago I wanted to added time averaged effect of a stirrer into momentum equation. My vector had correct dimension and worked perfectly:

tmp<fvVectorMatrix> tUEqn
(
    fvm::div(phi, U)
  + MRF.DDt(rho, U)
  + turb.divDevRhoReff(U)
 ==
    fvOptions(rho, U)
  + stirrerMomentum
);
fvVectorMatrix& UEqn = tUEqn.ref();

Good luck with finding the error,
Tom

[sturgeon]

Originally Posted by TomasDenk

Hi Sturgeon,

I'm not able to pin point where in your code is the mistake, because there are several symbols that are not defined in the code snippet. What I would do first if I were in your situation is dimension analysis. The equation for U is momentum equation and its dimension is N*m. So it seems that all components of your vector F should be in N/m. Make sure that all terms in Fz (which is sort of complex) are od the same dimension too.

When I was writing the code sometimes it would compile fine, but then upon execution of the solver I'd get a dimension error for the forces. So I thought if it was dimensionally inaccurate then it wouldn't be able to run at all?

I would also warn you not to use practice like

beta*(T-Tref)*mag(g)

very often. As long as g is defined as vector, I would prefer using it in a general way. Imagine a user that decides to define

g = (0, -9.81, 0) m/s^2

One term in your equation respects the choice

fvc::interpolate(rhok-1)*(g & mesh.Sf())

the other does not

beta*(T-Tref)*mag(g)

This is very confusing.

Yeah thanks, I am aware of this but can't really find a way around it aside from restricting choice of z = vertical axis.

Thanks for your input

[TomasDenk]

Yeah thanks, I am aware of this but can't really find a way around it aside from restricting choice of z = vertical axis.

My idea was to include gravity term in all three components of F with respect to acceleration components:

Fx   beta*(T-TRef)*(vector(1,0,0)&g))
Fy   beta*(T-TRef)*(vector(0,1,0)&g))
Fz   beta*(T-TRef)*(vector(0,0,1)&g))

[sturgeon]

Originally Posted by TomasDenk

Hi Sturgeon,

I'm not able to pin point where in your code is the mistake, because there are several symbols that are not defined in the code snippet. What I would do first if I were in your situation is dimension analysis. The equation for U is momentum equation and its dimension is N*m. So it seems that all components of your vector F should be in N/m. Make sure that all terms in Fz (which is sort of complex) are od the same dimension too.

Thanks for your help so far. I'm looking at the dimensions as you suggested and I'm a bit confused.

You say my vector F should be N/m. If I output the dimensions of F to the console, I'm getting ms^-2. But if I look at the other terms in the momentum predictor (so I use reconstruct individually on fvc::interpolate(rhok-1)*(g & mesh.Sf()) and fvc::snGrad(p)*mesh.magSf() and output their dimensions) they are also ms^-2. Surely my F must therefore be consistent with this?

[TomasDenk]

Originally Posted by sturgeon

You say my vector F should be N/m. If I output the dimensions of F to the console, I'm getting ms^-2. But if I look at the other terms in the momentum predictor (so I use reconstruct individually on fvc::interpolate(rhok-1)*(g & mesh.Sf()) and fvc::snGrad(p)*mesh.magSf() and output their dimensions) they are also ms^-2. Surely my F must therefore be consistent with this?

If what you write here is true, we might have found the problem. Your equation contains these terms:

fvc::interpolate(rhok-1)*(g & mesh.Sf())
+ Ff
- fvc::snGrad(p)*mesh.magSf()

The dimensions of the first and the third term are ms^-2. If the dimension of F is also ms^-2, then considering

Ff = fvc::interpolate(F)&mesh.Sf()

the dimension of Ff is m^3s^-2. As far as I know mesh.Sf() contains face normals with size of face area and dimension m^2. Can you check if Ff has the dimension I suspect?

[sturgeon]

Originally Posted by TomasDenk

If what you write here is true, we might have found the problem. Your equation contains these terms:

fvc::interpolate(rhok-1)*(g & mesh.Sf())
+ Ff
- fvc::snGrad(p)*mesh.magSf()

The dimensions of the first and the third term are ms^-2. If the dimension of F is also ms^-2, then considering

Ff = fvc::interpolate(F)&mesh.Sf()

the dimension of Ff is m^3s^-2. As far as I know mesh.Sf() contains face normals with size of face area and dimension m^2. Can you check if Ff has the dimension I suspect?

Yeah, it does, as does the dimensions of

fvc::interpolate(rhok-1)*(g & mesh.Sf()

and

fvc::snGrad(p)*mesh.magSf()

Cheers