CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Running, Solving & CFD

obtain velocity from pressure profile

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree6Likes
  • 1 Post By geth03
  • 1 Post By geth03
  • 1 Post By geth03
  • 1 Post By geth03
  • 1 Post By geth03
  • 1 Post By Tobermory

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   July 26, 2021, 18:17
Default obtain velocity from pressure profile
  #1
New Member
 
Naive CFD
Join Date: Jun 2021
Posts: 24
Rep Power: 3
Deep111090 is on a distinguished road
Hi All,
I asked this question in another forum also, but I haven't reached any conclusion yet. I apologize if you have to read this again.
I am solving a 2D lid-driven flow using icoFOAM algorithm, one of the tutorial cases.


I decided to experiment with it, the First step of which involved the calculation of the velocity and pressure up to convergence.


In the second step I modified the code of the icoFoam solver by commenting out the pEqn.solve() line (i.e. not solving the pressure) and used the pressure obtained in the first step to calculating the velocity profiles. In this step, the velocity was initialed from zeros with boundary conditions similar to the first step.


My assertion was that the velocity profile obtained in the second step should be exactly similar to that of the first step since they correspond to the same pressure field.


But the results are very much Reynolds number dependent. At low Reynolds number, I do have similar velocity profiles as of both the cases (though not exact match ~differ by order of 1e-4). But at Reynolds's number greater than 200 the velocity profiles are not remotely similar, i.e. the velocity fields obtained for the second step are incorrect and exhibit discontinuities in their profiles.


Can you please guide me in the right direction as why this is happening.

Thanks in advance.
PS: if my question is not clear please feel free to ask for more details.
Deep111090 is offline   Reply With Quote

Old   July 27, 2021, 09:44
Default
  #2
Senior Member
 
Join Date: Dec 2019
Location: Cologne, Germany
Posts: 290
Rep Power: 6
geth03 is on a distinguished road
if you want to know what your steady-state velocity profile looks like when pressure is pre-determined, you do not need to engage with pressure-velocity coupled solvers.

if you know your pressure field, just solve

fvVectorMatrix UEqn
(
fvm::div(phi, U)
- fvm::laplacian(nu, U)
);

solve(UEqn == -fvc::grad(p));

as many times as you need to reach steady-state.
Deep111090 likes this.
geth03 is offline   Reply With Quote

Old   July 27, 2021, 09:52
Default
  #3
New Member
 
Naive CFD
Join Date: Jun 2021
Posts: 24
Rep Power: 3
Deep111090 is on a distinguished road
Thanks Geth for the response. I tried that a few times. The problem here is that the "phi" is not conserving the mass locally. Right now the question before me is more fundamental that is it even possible to get a solution by this method.

Essentially this problem summarises to two eqautions (momentum in x and y) and two unknowns(u_x,u_y).

But the given the equations are non linear is it mathematically possible to reach the solution (I may be wrong while stating this line)
Deep111090 is offline   Reply With Quote

Old   July 27, 2021, 10:04
Default
  #4
Senior Member
 
Join Date: Dec 2019
Location: Cologne, Germany
Posts: 290
Rep Power: 6
geth03 is on a distinguished road
mathematically, it is true that the NS-equation is non-linear.

but openfoam does not consider that, it linearizes the convection term in the NS-equation.
fvm::div(phi, U) : look at that expression. here you can see that velocity transports itself.

this way you can actually use an iterative solver, and not need to engage with non-linear systems.
Deep111090 likes this.
geth03 is offline   Reply With Quote

Old   July 27, 2021, 10:12
Default
  #5
New Member
 
Naive CFD
Join Date: Jun 2021
Posts: 24
Rep Power: 3
Deep111090 is on a distinguished road
yes..now i have question that how shall i calculate phi ,..i am doing it by fvc::interpolate(U)&mesh.Sf().
But by default the openFoam solves it by phi = phiHby-pEqn.flux();

Question 1: is there any difference (as both approaches give slighty different answers for phi).
Question 2: Will this phi be locally mass conserving. as it is not if i don't solve pressure eqaution and just feed pre-determined pressure?
Deep111090 is offline   Reply With Quote

Old   July 27, 2021, 10:26
Default
  #6
Senior Member
 
Join Date: Dec 2019
Location: Cologne, Germany
Posts: 290
Rep Power: 6
geth03 is on a distinguished road
before entering the piso loop, you create phi with fvc::interpolate(U)&mesh.Sf().

in the pressure-velocity coupling part, phiHbyA is the volume flux derived from the UEqn, which does not take into account the pressure part, so it is overguessed.

fvScalarMatrix pEqn
(
fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
);

the pEqn is derived when inserting the momentum eqn in the mass conservation eqn, basically you are looking for a pressure field that conserves your mass flow. also called pressure poisson equation.

phi = phiHbyA - pEqn.flux();

here you correct phiHbyA, bc it was overguessed previously, to solve for a volume flux that considers a pressure field that conserves mass flow.

i hope you understand by now what each line is doing.

so if you do not solve for pEqn, you will always substract the same amount from your previous volume flux, which will change its value, of course.


to sum up, you solve for U with known pressure field, calculate then the volume flux without it and then correct the volume flux again with the known pressure field.
Deep111090 likes this.
geth03 is offline   Reply With Quote

Old   July 27, 2021, 12:23
Default
  #7
New Member
 
Naive CFD
Join Date: Jun 2021
Posts: 24
Rep Power: 3
Deep111090 is on a distinguished road
Hi I have tried the aforementioned method...my code looks like



fvVectorMatrix UEqn
(
fvm::div(phi, U)
- fvm::laplacian(nu, U)
);

solve(UEqn == -fvc::grad(p));


phi = fvc::interpolate(U)&mesh.Sf();


I got following results for Reynolds number 250. The phi is not conserving mass in this case. Hence the unrealistic velocity profiles. The results for using pre-determined pressure are attached as images
Attached Images
File Type: jpg Mag1.jpg (30.3 KB, 13 views)
File Type: jpg U_x.jpg (28.1 KB, 7 views)
File Type: jpg U_y.jpg (30.3 KB, 7 views)
Deep111090 is offline   Reply With Quote

Old   July 27, 2021, 12:27
Default
  #8
New Member
 
Naive CFD
Join Date: Jun 2021
Posts: 24
Rep Power: 3
Deep111090 is on a distinguished road
The actual results should look like as attached. I am not sure what is causing this phi to be not conserving mass locally. Even though we are solving from pressure which is predetermined and known to be correct.
Attached Images
File Type: jpg Mag2.jpg (30.2 KB, 5 views)
File Type: jpg U_x2.jpg (28.9 KB, 3 views)
File Type: jpg U_y2.jpg (30.8 KB, 3 views)
Deep111090 is offline   Reply With Quote

Old   July 28, 2021, 02:23
Default
  #9
Senior Member
 
Join Date: Dec 2019
Location: Cologne, Germany
Posts: 290
Rep Power: 6
geth03 is on a distinguished road
how many times did you execute the code?
only once?
did you try to solve it many more times? in a loop for example?
geth03 is offline   Reply With Quote

Old   July 28, 2021, 02:24
Default
  #10
New Member
 
Naive CFD
Join Date: Jun 2021
Posts: 24
Rep Power: 3
Deep111090 is on a distinguished road
Hi.. i looped it many times. the momentum euations had resuidals in range of 1e-6 and the continuty errs were high for local.
Deep111090 is offline   Reply With Quote

Old   July 28, 2021, 05:42
Default
  #11
Senior Member
 
Join Date: Dec 2019
Location: Cologne, Germany
Posts: 290
Rep Power: 6
geth03 is on a distinguished road
i think the main problem here is phi. the volume flux which is not known but is treated as known and then re-calculated again for a known pressure field.

i think as U_wall or U_top becomes larger and therefor Re becomes larger, the difference or discrepancy of this approach becomes more obvious.

the main problem of a non-linear set of equations is that for large systems it is much slower than the iterative approach, given an unknown pressure field. the way it is coded in OF both U and p are determined with the iterative approach or segregated, within a loop, to ensure mass conservation.

your problem is, that you already know your pressure field. so you don't need to the approach of fluid transporting itself, i guess.

i don't know if OF is capable of solving non-linear systems.
Deep111090 likes this.
geth03 is offline   Reply With Quote

Old   July 28, 2021, 13:35
Default
  #12
New Member
 
Naive CFD
Join Date: Jun 2021
Posts: 24
Rep Power: 3
Deep111090 is on a distinguished road
Thanks for the response.
Right now i am under the impression that it can not be done without solving for continuity..which seems to be tricky as I am unable to figure out the way to enforce continuity on "phi" without solving the Poisson's equation. Can you suggest some way to enforce continuity on each cell locally without solving Poisson's .
Deep111090 is offline   Reply With Quote

Old   July 30, 2021, 03:48
Default
  #13
Senior Member
 
Join Date: Dec 2019
Location: Cologne, Germany
Posts: 290
Rep Power: 6
geth03 is on a distinguished road
unfortunately not.

i think without the pressure-velocity coupling, your solution of the velocity for a given pressure field will be highly dependant of your initial values for velocity if you try to solve for velocity in the manner how it is implemented in openfoam.

if you have time, could you use different initial velocity fields to see how sensitive your end solution is for the initial velocity field?

you could
1) initialize your new solver with the end solution of the pressure and velocity field of icoFoam to see if things change
2) initiliaze with fields from a time you get with icoFoam, lets say endTime/2 or 0.75 x endTime.
Deep111090 likes this.
geth03 is offline   Reply With Quote

Old   July 30, 2021, 10:47
Default
  #14
New Member
 
Naive CFD
Join Date: Jun 2021
Posts: 24
Rep Power: 3
Deep111090 is on a distinguished road
Thanks Geth. You are right.

Its all my hunch and i would appreciate to be corrected. The problem here is coupling. It is way the openfoam is and use rhie chow interpolating in spirit of it not explicitly. I may be wrong but from what i can see is that openfoam sort of spread the elements of the rhie chow in entrie solution algorithm which causes problem if we comment some sections out.
The reason why it should not work with phi = fvc::interpolate(U)& mesh is because there is no coupling between phi and pressure due to this interpolation.
The reason it is not working with phi = phiHbyA-pEqn.flux() was al ready stated by you.

Last edited by Deep111090; July 31, 2021 at 14:32.
Deep111090 is offline   Reply With Quote

Old   July 30, 2021, 13:27
Default
  #15
Senior Member
 
Join Date: Apr 2020
Location: UK
Posts: 450
Rep Power: 11
Tobermory will become famous soon enough
I came here from your thread on Rhie Chow, and can understand now why you were asking about that. To reiterate:

Quote:
how shall i calculate phi ,..i am doing it by fvc::interpolate(U)&mesh.Sf().
But by default the openFoam solves it by phi = phiHbyA - pEqn.flux();
These 2 approaches generally do not give the same answer; the simple interpolation approach is incorrect, and will not conserve mass. The OpenFOAM solvers updates their face fluxes, phi, based on the momentum interpolation (also called Rhie Chow) approach which is the second method you mentioned above.

Coming back to your original experiment (which I approve, by the way - the best way to learn is to tinker - so well done!), you have only commented out the p.solve() line in the corrector loop ... the code is still calculating the momentumless face flux and then trying to correct this with the pressure gradient - but you have broken the link between the coefficients in the pEqn and the pressure values since you're not solving the pressure equation. This means that there's no pressure velocity coupling, as was observed above, and so I can see why the solution might go off the rails.

For example, in the normal solver, if you perturb the velocities in a few cells, then these will generate a pressure field that will ultimately dampen them down again since there is nothing there to maintain the pressure perturbation. In your solver you can perturb the velocities and yet no pressure perturbation occurs, and so no restoring force occurs. The only stabilising effect is viscosity ... which explains the Re dependence.
Deep111090 likes this.
Tobermory is offline   Reply With Quote

Old   July 30, 2021, 13:56
Default
  #16
New Member
 
Naive CFD
Join Date: Jun 2021
Posts: 24
Rep Power: 3
Deep111090 is on a distinguished road
Thanks for your response Tobermory. It seems that now we can reach some conclusion from ideas given by yourself and Geth. I am not sure which coefficients have have links broken. As the rAU remains the same through out all the iterations. And judging from the discussion on pEqn.flux() in this thread ( Description of flux() method ). The pEq.flux() depends upon (1/Ap) which is rAU and other thing it depends on is p.
Deep111090 is offline   Reply With Quote

Old   July 31, 2021, 06:44
Default
  #17
Senior Member
 
Join Date: Apr 2020
Location: UK
Posts: 450
Rep Power: 11
Tobermory will become famous soon enough
Quote:
but you have broken the link between the coefficients in the pEqn and the pressure values since you're not solving the pressure equation
To understand this, start by looking at the corrector loop. The pEqn is set up as follows:
Code:
fvScalarMatrix pEqn
   (
      fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
   );
This builds the fvMatrix for the pressure equation, populating the coefficients using phi, H and A - ie these depend on the velocity and flux fields, as well as the geometric constants (if mesh is stationary). The pEqn().solve would then solve the pressure field to give a p field that is "in equilibrium" with the matrix coefficients. You have commented out the solve, so the p-field and the pEqn matrix coefficients are no longer in balance - that is the broken link that I am referring to.

Later on in the corrector loop, the face flux is corrected with
Code:
   phi = phiHbyA - pEqn.flux();
where pEqn.flux() pulls the data directly from the pEqn fvMatrix (to avoid continuity errors) ... the problem is that the p values are not correct - they are not in balance with the face fluxes, and so will not preserve continuity.

Have a think of the analogy that I mentioned in an earlier post: if you were to perturb the velocities in a few cells, then you would expect to see some impact on the pressure field, which itself would generate some kind of feedback correction to the velocity field, smoothing out the perturbation. In your setup, you have disabled that feedback, so the velocity is free to wander.
Tobermory is offline   Reply With Quote

Old   July 31, 2021, 11:33
Default
  #18
New Member
 
Naive CFD
Join Date: Jun 2021
Posts: 24
Rep Power: 3
Deep111090 is on a distinguished road
Thanks I now have a cleared idea abt this.

Last edited by Deep111090; July 31, 2021 at 21:42.
Deep111090 is offline   Reply With Quote

Reply

Tags
cfd, icofoam, openfoam

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Wind tunnel Boundary Conditions in Fluent metmet FLUENT 6 October 30, 2019 12:23
serial udf to parallel udf radioss Fluent UDF and Scheme Programming 10 January 19, 2019 08:56
Unable to obtain the actual velocity profile of a sinusoidal udf Gourabc FLUENT 0 March 4, 2017 04:49
Boundary Conditions : Total Pressure or Velocity Gearb0x OpenFOAM Running, Solving & CFD 2 February 28, 2011 21:18
obtain non-uniform velocity profile Atit Koonsrisuk CFX 6 July 9, 2007 06:04


All times are GMT -4. The time now is 02:48.