View Issue Details

IDProjectCategoryView StatusLast Update
0000973OpenFOAMBugpublic2020-01-10 16:18
Reporteruser728Assigned Tohenry  
PrioritynormalSeverityminorReproducibilityalways
Status closedResolutionunable to reproduce 
Summary0000973: waveSurfacePressureFvPatchScalarField.C ImplizitEuler/CrankNicolson wrong?!
DescriptionI tried potentialFreeSurfaceFoam with the tutorial "oscillatingBox".
I tested IE (Implizit Euler ddt scheme) and CN 1.0 (CrankNicolson) ddt scheme with very small time steps. Here, periodial water height is similar over time. That is fine.
Larger timesteps causes a discrepancy in solution.

I think, treatment in boundary condition is wrong.

Original code with path OpenFOAM-2.2.x/src/finiteVolume/fields/fvPatchFields/derived/waveSurfacePressure/waveSurfacePressureFvPatchScalarField.C :

switch (timeScheme)
    {
        case tsEuler:
        case tsCrankNicolson:
        {
            zetap = zeta.oldTime().boundaryField()[patchI] + dZetap;

            break;
        }
        case tsBackward:
        {
            scalar dt0 = db().time().deltaT0Value();

            scalar c = 1.0 + dt/(dt + dt0);
            scalar c00 = dt*dt/(dt0*(dt + dt0));
            scalar c0 = c + c00;

            zetap =
                (
                    c0*zeta.oldTime().boundaryField()[patchI]
                  - c00*zeta.oldTime().oldTime().boundaryField()[patchI]
                  + dZetap
                )/c;

            break;
        }
        default:
        {
            FatalErrorIn
            (
                "waveSurfacePressureFvPatchScalarField<Type>::updateCoeffs()"
            ) << " Unsupported temporal differencing scheme : "
                << ddtSchemeName << nl
                << " on patch " << this->patch().name()
                << " of field " << this->dimensionedInternalField().name()
                << " in file " << this->dimensionedInternalField().objectPath()
                << abort(FatalError);
        }
    }

BC treats IE and CN in a similar way. I do not understand this. I would use something like:

switch (timeScheme)
    {
        case tsEuler:
        {
            
            zetap = zeta.oldTime().boundaryField()[patchI] + dZetap;
            //like old code
            break;
        }
        case tsCrankNicolson:
        {
            
            
            scalar coeff = zeta.mesh().ddtScheme(zeta.name())[1].number();
            vectorField dZetapOld(dt*nf()*phi.oldTime().boundaryField()[patchI]/patch().magSf());
            
            zetap = zeta.oldTime().boundaryField()[patchI] + coeff*dZetapOld/2.0 + (2-coeff)*dZetap/2.0;
//which blends between IE and CN

            //zetap = dZetap*0.5 + 0.5*dZetapOld+ zeta.oldTime().boundaryField()[patchI];//which would be pure CN
            
            
            break;
        }
        case tsBackward:
etc...

where coeff is the blending between IE (coeff==0) and pure CN (coeff==1.0).
(Maybe there is a more elegant way to set scalar coeff. Sorry about that.)

Best regards,

Andy


TagsNo tags attached.

Activities

There are no notes attached to this issue.

Issue History

Date Modified Username Field Change
2013-08-22 14:29 user728 New Issue
2020-01-10 16:18 henry Assigned To => henry
2020-01-10 16:18 henry Status new => closed
2020-01-10 16:18 henry Resolution open => unable to reproduce