View Issue Details
ID | Project | Category | View Status | Date Submitted | Last Update |
---|---|---|---|---|---|
0000973 | OpenFOAM | Bug | public | 2013-08-22 14:29 | 2020-01-10 16:18 |
Reporter | Assigned To | henry | |||
Priority | normal | Severity | minor | Reproducibility | always |
Status | closed | Resolution | unable to reproduce | ||
Summary | 0000973: waveSurfacePressureFvPatchScalarField.C ImplizitEuler/CrankNicolson wrong?! | ||||
Description | I 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 | ||||
Tags | No tags attached. | ||||