[olesen]

I've select rhoSimpleFoam as a starting point for using openFOAM for standard industrial applications - steady-state, k/eps, moderate compressibility - thus my apologies for posting with presumably banal questions.

Backtracking from my original geometry, I am now using blockMesh to create a very simple test case.
I assume that encountered stability problems arise from incorrect configuration, but I am uncertain as to what I should be changing.

The test geometry is a simple square duct (50mmx50mm) with a u-bend. The inlet is L/D = 1, the bend is R/D = 1, and the outlet is L/D = 10.

I've used surfaceNormalFixedValue for the inlet velocity and fixedValue (1.1e5) for the outlet pressure.

With U=15m/s, I obtain a physically reasonable result (relaxation parameters 0.3/0.7). At a slightly higher inlet velocity (U=20m/s), the solution seems to diverge rapidly regardless of the relaxation parameters (lowest tried was 0.05/0.3).

I really cannot figure out if the difficulty arises from the physics (outlet too close) or from incorrect specification of the BCs and/or schemes. Ideally the inlet would be specified as an integral mass-flux, which would also alleviate the divergence problem, but I cannot imagine that this is the sole factor here.

If anyone has a few minutes to look at the problem and see what is going awry, I would be greatly appreciate it. I could also offer the problem for a potential rhoSimpleFoam tutorial.

Thanks,
/mark

[olesen]

As a followup to my own post.
Thanks to Joern Beilke for suggesting initialization with potentialFoam and providing some settings for fvSolution that had worked for him.

Initializing with potentialFoam did help for lower (fixedValue) inlet velocities, but problems persist at higher flow rates. Underrelaxing rho (0.1) and p (approx. 0.05) seemed to delay the divergence problem, but not prevent it.

With help of "computeMassFlux.H" (found on the forum), the mass flow can be seen to steadily increase. The obvious explanation would be that increases in the density at the inlet result in a correspondingly increased mass flow at the outlet (via adjustPhi).

Using my own "integralInlet" boundary condition to force a fixed inlet mass flow, the solver runs reasonably stably for moderate Mach numbers.

I am still, however, puzzled about the pressure correction. In the pEqn.h we have the following lines:

AU = UEqn().A();
...
pEqn( fvm::laplacian(rho/AU, p) == fvc::div(phi) );

I can't figure out where the density corrections are hidden - cf. Eq. 10.12 Ferziger & Peric. Or am I missing something quite obvious?

/mark