[FMDenaro]

Originally Posted by Jonas Holdeman

There is, of course, a simpler method for solving Simbelmyne's problem of a cylinder in a driven cavity. Yesterday I cobbled together some code from a LDC with a mesh I had generated for a MHD problem. It solves the stationary problem, and lacks a number of things but it apparently works. I use divergence-free finite elements derived as the curl of a simple cubic Hermite stream function (this might be called a stream function-velocity method). I use the usual BCs on the square boundary and u=v=0 on the cylinder surface. I should have used the condition that the stream function be constant on the cylinder surface but did not have the time and ambition to code that constraint. But solving the pressure-free NSE gives a surface stream function constant to within 1%.

It seems, however, that the omission of the surface constraint is not sufficient for use of quartic and quantic stream function elements. Maybe I will check this out sometime. The three figures of mesh, stream contours, and vector velocity plot are attached. I show the steady result because I don't know Simbelmyne's initial conditions.

Hello, I think that the constant value of the stream-function on the obstacle is mandatory, otherwise the mass is no longer conserved. This problem can be relevant for unsteady simulation...

[Simbelmynë]

Hey,

So I have made a small comparison of the two methods mentioned in the thread, prior to the suggestion by Jonas Holdeman. The simulation is on a 101x101 mesh, Re=100 based on the lid length, transient, explicit, with case1 representing the method where I try to minimize the error in a sequential manner for the boundary update (1st order approximation to the normal derivative). case2 use the integral condition (1st order approximation to the normal derivative). I use a quad mesh in Fluent as a comparison.

Assuming that Fluent gives the best result here we see that the first case is quite close and while being qualitatively correct the second case is quite far off.

Probably I am implementing the integral condition in an incorrect way (I calculate the stream function boundary condition prior to the Poisson iteration, but I have seen absolutely no difference in calculating it inside the Poisson loop, for a number of test runs on smaller meshes. In case1 I also calculate the stream function boundary condition prior to the Poisson loop).

[FMDenaro]

Originally Posted by Simbelmynë

Hey,

So I have made a small comparison of the two methods mentioned in the thread, prior to the suggestion by Jonas Holdeman. The simulation is on a 101x101 mesh, Re=100 based on the lid length, transient, explicit, with case1 representing the method where I try to minimize the error in a sequential manner for the boundary update (1st order approximation to the normal derivative). case2 use the integral condition (1st order approximation to the normal derivative). I use a quad mesh in Fluent as a comparison.

Assuming that Fluent gives the best result here we see that the first case is quite close and while being qualitatively correct the second case is quite far off.

Probably I am implementing the integral condition in an incorrect way (I calculate the stream function boundary condition prior to the Poisson iteration, but I have seen absolutely no difference in calculating it inside the Poisson loop, for a number of test runs on smaller meshes. In case1 I also calculate the stream function boundary condition prior to the Poisson loop).

I suspect something is wrong in the implementation, you should see a great difference when the zero circuitation constraint is computed along with the elliptic solver for Psi. Consider that the Psi value on the obstacle is a Dirichlet bc.s that dynamically adjusts its value depending on the iteration you have in the interior points of cavity. So, please:

1) provide the plot of the stream-function field
2) what about the fixed value of the Psi on the cavity walls?
3) what about the resulting value of Psi on the obstacle?
4) what about the residual of the circuitation after the Poisson solver converged?

Finally, I do not agree that Fluent can be assume as reference for the solution, did you check a grid refinement (100, 200, 400, etc..) and used a second order scheme?

[Simbelmynë]

Originally Posted by FMDenaro

I suspect something is wrong in the implementation, you should see a great difference when the zero circuitation constraint is computed along with the elliptic solver for Psi. Consider that the Psi value on the obstacle is a Dirichlet bc.s that dynamically adjusts its value depending on the iteration you have in the interior points of cavity. So, please:

1) provide the plot of the stream-function field
2) what about the fixed value of the Psi on the cavity walls?
3) what about the resulting value of Psi on the obstacle?
4) what about the residual of the circuitation after the Poisson solver converged?

Finally, I do not agree that Fluent can be assume as reference for the solution, did you check a grid refinement (100, 200, 400, etc..) and used a second order scheme?

I will try to answer the questions. Case 2 is the integral condition approach (which I might have implemented incorrect). Case 1 is my other approach. Both cases have been tested on the same mesh for a simulation time of 10 s (well past the time to reach steady-state).

Q1-Q3: I have attached images with stream function and vorticity for the two cases.

Q4: I'll post it soon. I imagine that you are interested in the results when I use the integral condition inside v.s. outside the Poisson loop.

Q5 (fluent question): I ran the case with a second order scheme and I attach a (simple) mesh sensitivity comparison.

Thank you for your input.

[FMDenaro]

well, I see... here some couple of questions:

1) Case 1, have you computed the mass rate flowing throught the obsacle due to the numerical errors? I suppose you get not exactly the same value of Psi on the obstacle.

2) Case 2, you implemented the zero circuitation condition using a first order derivative and a mid-rule for the line integral, right? that corresponds to approximate the cylinder with a somehow shlightly greater radious of the cylinder. You should try a second order forward discretization for the normal derivative and a trapezoidal rule for the integral.

Then, generally from the picture (which are both physically plausible) I cannot say which is the best solution, what about the grid resolution you used for Fluent, have you ensured that the cell Reynolds number is lesser than 1 in all cells? do you set a second order central scheme or second order upwind??

[Simbelmynë]

1. I don't think I have any mass flow through the cylinder for Case 1 since the stream function is constant along the cylinder wall (all my nodes have the same stream function on the wall). The method I use sets a stream function on the cylinder wall based on one of the near-boundary nodes. In that cell the (1st order) normal derivative is zero and hence the v_tan condition is also satisfied there. Since it is not satisfied in all other nodes on the cylinder I try to find the node which causes the least maximum error for the normal derivative.

2. I used your lecture notes, eqn. 68 for the condition, where dX=constant and dY=constant so that

sum_of_all_near_wall_nodes_psi/number_of_nodes_on_the_wall=psi_on_the_wall

I'll try a second order method and see if that changes anything.

When I implement this inside the Poisson loop I just calculate the sum and set the value of Psi on the wall before I do each sweep with the SOR method that I use. I have a comparison of the integral condition (case 2) implemented inside the Poisson loop and outside the Poisson loop, for a 31x31 mesh in a couple of minutes.

Edit: OK so now I have attached a figure with the difference of having the integral condition inside (case2_a) and outside (case2_b) of the Poisson loop. This is at steady-state, I have not extracted any information during the transient stage.

[FMDenaro]

Ok, results are equal, do you have difference in the number of iterations for the convergence of the Poisson solver?

Furhtermore, from the pattern of coulour I see you have Psi_wall=0 fixed on the lid edges and you get (more or less) on the obstacle

Case 1 Psi_obs=-0.045
Case 2 Psi_obs=-0.035

right?
Differences appear quite meaningful... I am thinking about a test to check and assess the codes. What about if you shift the cylinder towards one of the walls? If you have the surface almost tangent to a wall, a correct solution should provide Psi_obs -> Psi_wall. This way you should have an idea of the best approach

[Simbelmynë]

Thank you for good suggestions, I'll try and see what happens!

I also tried to run case1 inside/outside the Poisson loop. While it gives good results for 101x101 nodes it is very far off for 31x31 nodes (does not capture major flow circulation on top of cylinder) when used outside the Poisson loop.

However, when used inside the Poisson loop it provides good results. I thought that it would take much longer time to run the case, but in fact it is much much faster when the Case1 BC is used inside the Poisson loop. So the SOR solver seems to like that kind of BC, even if I don't see how it is different compared to the integral condition with respect to convergence (both are Dirichlet / Global Neumann type conditions).

The flow field experiences some small disturbances from time to time (only minor wiggles of the velocity vectors), that I do not see when the integral condition is used (case2). It does not seem to cause stability problems though.

Simulation time in seconds (yeah it is Matlab with zero optimizations ):
Case2_a: 709
Case2_b: 773
Case1_a: 165 (!!!)
Case1_b: 2034

[FMDenaro]

I am not surprised of such difference, using SOR you are updating the new computed solution X(i)k+1 and substitute it as soon as it is available in the other components, so when you insert inside it converges faster

[Simbelmynë]

So after running a (much much longer) simulation using a second order approximation for the normal derivative in the integral condition it is clear that the order has a substantial impact on the results.

The converged solution (the Poisson solver residuals flattened after approx 8 seconds of simulation time) had a residual on the wall of 1e-15. I calculated the residual as the sum of the derivatives on the wall.

Compared to the 100x100 fluent simulation a 65x65 case gave quite nice results.

So right now I am leaning towards the integral condition (2nd order). However, since it is 10-20 times slower than the original scheme (minimization) it is a bit difficult to judge which I should use.

[FMDenaro]

Originally Posted by Simbelmynë

So after running a (much much longer) simulation using a second order approximation for the normal derivative in the integral condition it is clear that the order has a substantial impact on the results.

The converged solution (the Poisson solver residuals flattened after approx 8 seconds of simulation time) had a residual on the wall of 1e-15. I calculated the residual as the sum of the derivatives on the wall.

Compared to the 100x100 fluent simulation a 65x65 case gave quite nice results.

So right now I am leaning towards the integral condition (2nd order). However, since it is 10-20 times slower than the original scheme (minimization) it is a bit difficult to judge which I should use.

I think that using in the SOR the correct relaxation factor you will get a fast convergence ;-)

[FMDenaro]

Originally Posted by Simbelmynë

So after running a (much much longer) simulation using a second order approximation for the normal derivative in the integral condition it is clear that the order has a substantial impact on the results.

The converged solution (the Poisson solver residuals flattened after approx 8 seconds of simulation time) had a residual on the wall of 1e-15. I calculated the residual as the sum of the derivatives on the wall.

Compared to the 100x100 fluent simulation a 65x65 case gave quite nice results.

So right now I am leaning towards the integral condition (2nd order). However, since it is 10-20 times slower than the original scheme (minimization) it is a bit difficult to judge which I should use.

forget to say that the integral condition discretized at second order requires not only the second order derivative but also the second order discretization of the integral, for example the trapezoidal rule. And you shoud consider this discretization for the residual

[Simbelmynë]

Originally Posted by FMDenaro

forget to say that the integral condition discretized at second order requires not only the second order derivative but also the second order discretization of the integral, for example the trapezoidal rule. And you shoud consider this discretization for the residual

Yes I will try that, now it is midpoint, but I still see a substantial improvement.

Originally Posted by FMDenaro

I think that using in the SOR the correct relaxation factor you will get a fast convergence ;-)

Will try to vary the relaxation factor as well

[FMDenaro]

start trying a factor between 1.6-1.9, the integral condition changes in a relevant manner the structure of the matrix

[Simbelmynë]

Originally Posted by FMDenaro

start trying a factor between 1.6-1.9, the integral condition changes in a relevant manner the structure of the matrix

I notice some speed up from the default 1.7 to 1.8, but it is nowhere near the simulation speed of the other method (which, by the way, seems rather insensitive to the choice of relaxation parameter).

----

I have also tried to update the vorticity implicitly for the case1 and it results in a near-perfect match with the Fluent results, even on a much coarser mesh. Running the same test now on case2 (integral BC).

[FMDenaro]

Originally Posted by Simbelmynë

I notice some speed up from the default 1.7 to 1.8, but it is nowhere near the simulation speed of the other method (which, by the way, seems rather insensitive to the choice of relaxation parameter).

----

I have also tried to update the vorticity implicitly for the case1 and it results in a near-perfect match with the Fluent results, even on a much coarser mesh. Running the same test now on case2 (integral BC).

ok, there exists a theoretical estimation for the optimal value but it requires the spectral radius of the transition matrix.
I suggest to do a test, do only the first time step and check convergence for value going from 1 to 2, step .1. Once found the minimum refine the step a try around this value. Usually 3 digits of the value allow to use an optimal relaxation factor during all the simulation

[Simbelmynë]

Originally Posted by Jonas Holdeman

There is, of course, a simpler method for solving Simbelmyne's problem of a cylinder in a driven cavity. Yesterday I cobbled together some code from a LDC with a mesh I had generated for a MHD problem. It solves the stationary problem, and lacks a number of things but it apparently works. I use divergence-free finite elements derived as the curl of a simple cubic Hermite stream function (this might be called a stream function-velocity method). I use the usual BCs on the square boundary and u=v=0 on the cylinder surface. I should have used the condition that the stream function be constant on the cylinder surface but did not have the time and ambition to code that constraint. But solving the pressure-free NSE gives a surface stream function constant to within 1%.

It seems, however, that the omission of the surface constraint is not sufficient for use of quartic and quantic stream function elements. Maybe I will check this out sometime. The three figures of mesh, stream contours, and vector velocity plot are attached. I show the steady result because I don't know Simbelmyne's initial conditions.

Hey Jonas,

Is it possible for you to plot the x-velocity along the vertical axis through the geometric center of the cavity, for Re=100 (based on the cavity length)? It would be interesting to see how your simulation compares to the other results in this thread.

[Jonas Holdeman]

Originally Posted by Simbelmynë

Hey Jonas,

Is it possible for you to plot the x-velocity along the vertical axis through the geometric center of the cavity, for Re=100 (based on the cavity length)? It would be interesting to see how your simulation compares to the other results in this thread.

Attached is my result for 41x41 mesh, Re=100. It is visually identical to some of your results.

[Simbelmynë]

Originally Posted by Jonas Holdeman

Attached is my result for 41x41 mesh, Re=100. It is visually identical to some of your results.

Thank you Jonas. I think your simulation is more accurate and I assume that it converges to that particular solution for larger grids. It is also nice to have a completely different method producing similar results. Do you have any reference to the stream function velocity method that you use? Also, have you tested your BC approach against flow over cylinder or similar, where we have some experimental data to compare with?

[Jonas Holdeman]

Originally Posted by Simbelmynë

Thank you Jonas. I think your simulation is more accurate and I assume that it converges to that particular solution for larger grids. It is also nice to have a completely different method producing similar results. Do you have any reference to the stream function velocity method that you use? Also, have you tested your BC approach against flow over cylinder or similar, where we have some experimental data to compare with?

References are:
A Hermite finite element method for incompressible fluid flow.
Int. J. Numer. Meth. Fluids, 64:376--408, 2010.

Computation of incompressible thermal flows using Hermite elements.
Comput. Meth. Appl. Mech. Engrg., 199:3297--–3304, Dec 2010.

A velocity-stream function method for three-dimensional incompressible fluid flow.
Comput. Meth. Appl. Mech. Engrg., 209-212:66--73, Jan 2012.

I have done steady flow over a cylinder, but not unsteady flow. I do not recommend using just the u=v=0 BC on the cylinder. It works for my cubic Hermite elements with stream function and velocity DOFs, but that was just to get a quick result from existing (Matlab) code. But, the method is very robust with regard to BCs. In the flow over a step, you can truncate the mesh through a recirculation bubble with no ill effects (see first reference). I am working on a correct implementation of constant psi constraint on cylinder.

A neat thing about the method is that the solution gives not just the velocity, but also the stream function for plotting, without (or perhaps in place of) the pressure. Other quartic and quantic elements include second derivatives of psi DOFs, from which you can find the vorticity at the nodes directly.