[Donton]

Hello!!
I am trying to write a C++ code to simulate 2D flow past a fixed circular cylinder using the immersed boundary method with finite volume and staggered grid, fractional step method. I am trying to use the unsteady outflow boundary condition:
du/dt + U_1 * (du/dx) = 0
dv/dt + U_1 * (dv/dx) = 0

Now my question is how to calculate this U_1 so that the mass flow rate through the inlet boundary and the outlet boundary remain the same. Can anyone please suggest any useful paper or book or thesis regarding this outflow boundary condition implementation?

[selig5576]

The way you should do it (on a collocated grid) is as follows

Given we have



You can calculate the integral just by using sums. Typically you should also add a additive correction factor.

EDIT: I noticed you said staggered grid. In that case, just use the velocity instead of the face velocity.

[sbaffini]

It probably depends from the specific method you are using. In the cell centered FV context that U1 actually is the massflux on the outlet faces (there is a constant rho missing).

The idea is very simple, but I'll put it down in the most general way (so that it can be useful in general). You know the entering massflux as you have the velocity bcs at all the inlet and wall boundaries. Once you integrate the massflux on all those boundaries you have the total mass flux entering your domain, say, Min. You want that the massflux on all the other boundaries, Mout, exactly cancels the entering one, Min. So, Mout = -Min.

Let's say that you have N outlet boundaries and you want each of them to handle a fraction wi of the overall Mout. Each face j of the total Mi faces on each outlet boundary i has area A_ji and the total area for each outlet boundary i is A_i = SUM(A_ji,j=1,Mi).

You can then write down for each outlet boundary i:

-wi*Min = SUM(A_ji*rho_ji*U_ji,j=1,Mi)

However, you want to use a single velocity and a single density for all the faces on boundary i. Let's then assume that:

rho_ji*U_ji = rho_i * U_i = m_i (independent from j, the generic face on boundary i)

so that you can take it outside the sum and finally obtain:

m_i = rho_i * U_i = -wi*Min/SUM(A_ji,j=1,Mi) = -(wi/Ai) * Min

We can see that this is correct because now:

Mout = SUM(m_i*Ai,i=1,N) = -Min*SUM(wi,i=1,N) = -Min

assuming, obviously, that you correctly choose the wi so that they sum to 1.

In your very simple case with 1 inlet, 1 outlet and rho = cost = 1, you simply have that:

U_1 = - Min/Aout

where Min is just the integral of the normal to boundary velocity on the inlet boundary and Aout is the area of the outlet boundary. Note that the sign matters here for the massflux. Assuming you just work with velocities, you just get:

U_1 = Min/Aout

[sbaffini]

Note that, according to the way you implement the immersed boundary, you might need to actually take into account the massflux trough the immersed boundary because, yes, some immersed boundary implementations don't actually have a null mass flux trough solid walls. In that case you could actually use the method above to first corrrect the massflux on the immersed walls and later apply it for the outlet.

Also note that you can either apply it directly, as per your original question and as explained above, or as an additive correction. That is, you first compute the outlet massflux with upwind. That won't preserve the mass in general. You compute the mass imbalance and split it among the outlet faces.

[FMDenaro]

Originally Posted by Donton

Hello!!
I am trying to write a C++ code to simulate 2D flow past a fixed circular cylinder using the immersed boundary method with finite volume and staggered grid, fractional step method. I am trying to use the unsteady outflow boundary condition:
du/dt + U_1 * (du/dx) = 0
dv/dt + U_1 * (dv/dx) = 0

Now my question is how to calculate this U_1 so that the mass flow rate through the inlet boundary and the outlet boundary remain the same. Can anyone please suggest any useful paper or book or thesis regarding this outflow boundary condition implementation?

If you want a constraint on the mass I suppose you are simulating an incompressible flow problem. Therefore the divergence-free constraint is in act as continuity equation, right? In such a case, the correct mass constraint requires that the BC at the outflow is correctly implemented into the pressure equation. Personally I don't like this BC.s as it is a condition for inviscid flow. Note also that du/dx=-dv/dy from the divergence free constraint in 2D.
Have also a look to this old post Convective outlet boundary condition for Unsteady flows

[sbaffini]

I would distinguish the definition of U_1 as the global outlet-boundary velocity which will preserve the massflux (but, as I wrote, one can also work locally with a correction to the upwind flux) from the actual boundary condition in use at the boundary (as emerges from the implementation suggested in OF post).

However, in both cases that global velocity (or the local one + correction) is the one that needs to enter into the pressure equation.

[Donton]

Thank you everyone for your suggestions. I have tried as you have suggested. But still there is some issue.
Total mass flux at inlet Min and total mass flux at outlet Mout still is not exactly same.
For a test case when I am specifying some constant velocity of u = 1 and v = 0 both at inlet and outlet then (Min - Mout) is exactly zero and I am getting converging and good results. The residual flux in the computational is alwayas decreasing.
But when I am specifying the velocity only for the inlet and letting the outlet velocities get calculated from the upstream by the unsteady convective velocity boundary condition or by the neumann boundary condition in that case (Min - Mout) is not exactly zero rather of the order of 10^(-7). So in this case my pressure is gradually increasing with iteration and after few iteration pressure reaches to 10^6 order and finally diverging to nan.The residual flux in the computational domain is decreasing initially but when it reaches of the order of 10^(-3) it starts oscillating aperiodically with the same magnitude order and not decreasing.

Anyone has any suggestion? Alos if I am using unsteady convective velocity boundary condition at the outlet then what should be the boundary condition for the pressure correction equation at the outlet.

[FMDenaro]

I think that you set something wrong in the pressure equation when you let the outlet condition free. Check only after 1 time step the divergence-free constraint cell-by-cell and see of you get a higher error close to the outlet.

[sbaffini]

When you use upwind to determine the local outflow velocity, it won't preserve mass in general, as already written. Now, one option is the one you used in your first experiment, using a constant velocity on the whole boundary (in this case just equal to the one at inlet), and it seems to work.

However, at this point, I have a first question. What did you use as b.c. for the pressure equation at outlet in this first test?

Moreover, as already written by Filippo, that unsteady convective bc is purely inviscid. In my code (which, however, is density based and doesn't use a massflow correction at outlets), for example, I found that including the tangential derivatives at the outlet is fundamental to have a smooth exit in a simple channel without pressure spikes at outlet.

A second option is to just extrapolate the outlet velocity from upwind cells (so that it is not anymore constant on the whole outlet) and correct it locally, so that the sum of upwind velocities + corrections exactly equals Min. In order to do this, you have to:

1) Determine Min by integrating the velocity on inlet boundary.
2) Determine Mout* by integrating the upwind determined velocity U_j* (j indicating the j-th face on outlet) on outlet boundary.
3) On each face of the outlet correct U_j* to get the mass preserving U_j via:

U_j = U_j* + (Min-Mout*)/SUM(A_k)

So that

SUM(U_jA_j) = SUM(U_j* A_j + (Min-Mout*)A_j/SUM(A_k)) = Mout* + Min - Mout* = Min

Note again that signs matter here. You have to use this U_j as bc in the pressure equation.

Note also that this has nothing to do with the actual bc you use in the momentum equation, as long as you actually use this U_j (or the global version) as actual massflux at outlet.

[sbaffini]

Note that the local correction above can also be generalized as follows:

U_j = U_j* + (Min-Mout*)w_j/(A_j SUM(w_k))

So that:

SUM(U_jA_j) = SUM(U_j* A_j + (Min-Mout*)w_j/SUM(w_k)) = Mout* + Min - Mout* = Min

where w_j can itself be, for example, equal to A_j, or A_j^2, or anything you think could provide a better distribution of the correction over the faces of the outlet.

[Donton]

Finally my code seems to work with the suggestions you gave. Thank you all.

For my first test case BC s were like:
at inlet u=1, v=0, dp/dx = 0
at top and bottom wall du/dy=0, v=0, dp/dy = 0
at outlet u=1, v=0, dp/dx = 0
and it worked.

For my second case BC s were like:
at inlet u=1, v=0, dp/dx = 0
at top and bottom wall du/dy=0, v=0, dp/dy = 0
at outlet du/dx=1, dv/dx=0, dp/dx = 0
it did not work.

For my third case BC s were like:
at inlet u=1, v=0, dp/dx = 0
at top and bottom wall du/dy=0, v=0, dp/dy = 0
at outlet du/dt + U1 * du/dx = 0, dv/dt + U1 * dv/dx = 0, dp/dx = 0
it did not work.

For my fourth case BC s were like:
at inlet u=1, v=0, dp/dx = 0
at top and bottom wall du/dy=0, v=0, dp/dy = 0
at outlet du/dt + U1 * du/dx = 0, dv/dt + U1 * dv/dx = 0, p = 0
Finally it worked both with the correction and without the correction you mentioned in the last reply.

Can anyone say why for the 2nd and 3rd case the code did not work. Are the boundary conditions for both velocity and pressure correct in those two cases or they are wrong?

[sbaffini]

For second case, du/dx = 1 at outlet I think is quite arbitrary (and wrong).

For the other cases, what I'm concerned with is the pressure gradient at outlet.

What I remember as correct bc for the pressure is:

dp/dn = n . (Un*-Un)/dt

where Un* is the velocity normal to boundary before the pressure correction and Un is the one that should be mass preserving. So you see that you could actually apply the massflow correction just within the pressure equation (the bc actually). Otherwise you can apply it just before the pressure equation, in which case dp/dn = 0 should be correct.

What do you mean by not working? Pure blow up, lack of convergence or convergence toward a (presumably) wrong solution?

Also, what does it practically mean for you to use:

du/dt + U1 * du/dx = 0, dv/dt + U1 * dv/dx = 0

as bc?