[rolando]

Hello Hrvoje,
I tried the version you sent me and experimented with it.
For my testcase I noticed two things.
The testcase is a mesh moving, with a constant velocity and a relative flow in the moving direction.
I use the icoDyMFoam solver, where I activated the fvc::ddtPhiCorr() term in the PISO loop and I use the backward time scheme.

1) If I use a dynamic time-step the time-step adjustment is far more stable than without using this term. But it is still sensitive to the initial time-step (see above).

2) If I use a constant time-step the results are slightly worse than without the fvc::ddtPhiCorr() term. This refers to the solutions of U and p. Phi and meshPhi however is calculated correctly.
This is caused by the method fvcDdtPhiCoeff(), which seems to be a weighting factor in fvc::ddtPhiCorr(). The problem appears, as in this method the difference between the "real fluxes" and the fluxes built by the velocity is slightly different.
To overcome this trouble, I introduced an additional artificial weighting factor, which accounts for the change of phi in the past.
backwardDdtScheme<type>::fvcDdtPhiCorr then reads:

if (mesh().moving())
{
.
.

return tmp<fluxfieldtype>
(
new fluxFieldType
(
ddtIOobject,
rDeltaT
*fvcDdtPhiCoeff(U.oldTime(), phi.oldTime())
*min(mag(phi.oldTime() - phi.oldTime().oldTime()) / mag(phi.oldTime()), 1.0)
*(
(
fvc::interpolate(rA*V0oV)
*coefft0*phi.oldTime()
- fvc::interpolate(rA*V00oV)
*coefft00*phi.oldTime().oldTime()
)
- (
fvc::interpolate
(
rA*
(
coefft0*U.oldTime()*V0oV
- coefft00*U.oldTime().oldTime()*V00oV
)
) & mesh().Sf()
)
)
)
);

.
.
}

For my case this seems to cure the problem and the results of p and U are what they should be.

Rolando

[lr103476]

Hi everybody,

I also observed the problem with backward scheme in combination with a variable timestep. Using a custom dynamicFvMesh library I let the internal mesh points oscillate to create changing cell volumes. In such a way, using a uniform freestream, I could check if the solution is not influenced by the mesh motion and if the Geometric Conservation Law is satisfied.

Summarising, I got no problems with CrankNicholson or Euler using constant and variable timesteps. Using backward I obtained a different from freestream solution when the timestep was varying. With constant timestep the solution was OK with the backward scheme.

Together with a colleague (who knows more about time integration schemes than I do) we found a solution to this problem, which lies in backwardDdtScheme::meshPhi.

I will briefly explain the principle of the solution. The ddt is properly defined as:
phi(n+1) = ((3/2)V(n+1) - 2V(n) + (1/2)V(n-1)) / dt_1

This is solved using the values of the flux phi on timesteps n-1/2 and n+1/2:

phi(n-1/2)=(V(n)-V(n-1))/dt_0
phi(n+1/2)=(V(n+1)-V(n))/dt_1

Here dt_1 is the current and dt_0 the previous timestep. As far a we could see you need to make a correction in meshPhi for the fact that you have to use the dt_0 in the above line. Therefore we multiplied the coefft00 from backwardDdtScheme::meshPhi with dt_0/dt_1. When the timestep is constant, i.e. dt_0=dt_1, this ratio equals 1 and you will get the correct solution.

So we changed:
scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));

Into:
scalar coefft00 = deltaT/(deltaT + deltaT0);

Now I obtain the correct solution using backward differencing on a moving mesh with variable timestep.

Regards, Frank

[henry]

Did you make this change to coefft00 only in backwardDdtScheme<type>::meshPhi or also in all the ddt functions in backwardDdtScheme.C?

Henry

[henry]

I have redone the derivation of the backward scheme for moving meshes and I agree with your change and will include it in the 1.4 release. Thanks a lot for the help in resolving this problem.

Henry

[lr103476]

Allright. From now on we can safely use this second-order backward scheme.

If you still need the answer to your last question, here it is....I only changed the coefft00 in the backwardDdtScheme::meshPhi. In all ddt functions the definition of coefft00 is still OK, since at those places the current timestep dt_1 was used, according to the definition of ddt. Only there were you need the 'n+1/2' and 'n-1/2' values, which is in meshPhi, this definition of coefft00 needs a correction.

Frank

[henry]

Yes I agree.

Henry

[rolando]

Thanks a lot for that post Frank,
I tested it and the timestep oscillations are no longer present.

Rolando