[wildfire230]

Hi all,


I am trying to solve a 2D PDE in time h = h(x,y,t) which is given here:


This corresponds to basically a type of blistering film problem. I am able to implement a solver for this by modifying the scalarTransportFoam solver with the following HEqn:

fvScalarMatrix HEqn
(
fvm::ddt(H)
-fvc::laplacian(Foam::pow(H,3.0),fvc::laplacian(fvc::laplacian(H)))
==
fvOptions(H)
);

HEqn.relax();
fvOptions.constrain(HEqn);
HEqn.solve();
fvOptions.correct(H);

Unfortunately, this scheme has extremely bad stability characteristics with explicit discretization and I need timesteps 1e-10 just for stability, which is impractical. Is there any way to rewrite this with implicit discretization? It doesn't seem to like switching any of the fvc to fvm.


Thanks

[HPE]

fvm::laplacian instead of fvc::laplacian? Otherwise, adding as a source term, I suppose.

[wildfire230]

Hi HPE,


Thanks, it seems like you cannot use fvm on any nested operations. For example, if you want the biharmonic operator, it seems like it will only compile as fvc::laplacian(fvc::laplacian(H)). It throws errors whenever I replace either of these fvc with fvm. Can you please expand on your suggestion to replace this as a source term? Basically implement the entire equation as a source term? dh/dt = F and then implement that as a calculated source term with fvOptions, or something like that? I'm not sure if that will solve the implicit vs. explicit calculation.


Thanks again

[HPE]

Can you construct the inner laplacian outside of fvMatrix, make the outer laplacian fvm::, and pass it as a vol?Field into the outer laplacian?

If H is a volScalarField, laplacian of laplacian of H will be again a volScalarField. But this new H will be H*, and will break down the fvMatrix.

Hmm.. So the main problem is H* rather than the inner laplacian oper. Let me think further.

[HPE]

That term is currently a source term. Can you put it to rhs, with fvm::SuSp(your term), and test the stability again?

[HPE]

You may also want to implement your special laplacianLaplacian operator directly instead of executing laplacian twice, which would be expensive if you concern over the performance.

[wildfire230]

Hi HPE,


Thanks again for the advice. For fvm::SuSp, is that something that I can use directly in the modified solver? For example, if I had the biharmonic operator, would that be implemented as:


fvm::SuSp(fvc::laplacian(fvc::laplacian(H))


I may be misunderstanding. For directly implementing something like laplacianLaplacian as a new operator, is that relatively straightforward? To be honest, I have never messed with working with these operators and I am not exactly sure where they are even implemented in OpenFOAM's directory structure.


Thanks again.

[HPE]

About the new operator: just leave it then. Your sims will be considerably fast enough.

So out of my head something like this:

tmp<volScalarField> HStar
(
fvc::laplacian
(
pow3(H)*fvc::laplacian(fvc::laplacian(H))
)
);

In fvMatrix' rhs:

fvm::SuSp(HStar/H, H);

or

tmp<volScalarField> HStar
(
pow3(H)*fvc::laplacian(fvc::laplacian(H))
);

in fvMatrix's left hand side

-fvm::laplacian(Hstar/H, H);

Division and multiplication by H is an unfortunate trick, but it is being used in OpenFOAM anyways.

So, I havent compiled them as I have been writing on the phone (as always), please let me know if it would be any help.

[HPE]

Oh, and dont forget to bound H at the constructor stage as well as after fvMatrix to avoid the division by zero error:

bound(H, HMin); where you define HMin somewhere like:

const dimensionedScalar HMin("HMin", H.dimensions(), SMALL);

[wildfire230]

Hi HPE,


Thanks again, here is a quick update on your recommendations. Your first suggestion of :


tmp<volScalarField> HStar
(
fvc::laplacian
(
pow3(H)*fvc::laplacian(fvc::laplacian(H))
)
);

In fvMatrix' rhs:

fvm::SuSp(HStar/H, H);


DOES compile, although it seems like it does not really offer any improvement to the stability. I'm not sure if it is still basically calculating everything explicitly and so failing for that reason or something more subtle. Your second suggestion of:


tmp<volScalarField> HStar
(
pow3(H)*fvc::laplacian(fvc::laplacian(H))
);

in fvMatrix's left hand side

-fvm::laplacian(Hstar/H, H);


does not actually compile for me. It gives me a "no matching function call error" for ‘Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> >::tmp(Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> >, Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> >)’


I'm not sure if you would have any other suggestions? In any case thanks a lot already.

[HPE]

Let me look at them when I sit in front of a PC, then, since I was always writing via my phone (always delaying good quality answers). The second option's problem is curious.

The first option should have helped, but may be the division causes a new issue.

Which numerical scheme set do you use?

[wildfire230]

Thanks, currently I have the solver setup like this:

tmp<volScalarField> HStar
(
fvc::laplacian(Foam:ow(H,3.0),fvc::laplacian(fvc ::laplacian(H)))
);

fvScalarMatrix HEqn
(
fvm::ddt(H)
- UnitCheck * fvm::SuSp(HStar/H, H)
==
fvOptions(H)
);



And the schemes I am using are:
ddtSchemes: backward
gradSchemes: Gauss linear
divSchemes: none
laplacianSchemes: Gauss linear corrected
interpolationSchemes: linear
snGradSchemes: corrected

[wildfire230]

Also just FYI, it seems like the only schemes actually being used are ddtSchemes and laplacianSchemes, the rest I can all set to default none and everything runs.

[LuckyTran]

Originally Posted by wildfire230

Also just FYI, it seems like the only schemes actually being used are ddtSchemes and laplacianSchemes, the rest I can all set to default none and everything runs.

That is correct. Your equation has one ddt and one bilaplacian (which you are currently tinkering with in a laplacian approach). The other schemes are not needed and can be set to anything. You might need the in the future depending on how your implementation finally ends up working.

[wildfire230]

Does anyone have an idea of how I could implement an actual bilaplacian operator into openfoam so that I could directly use implicit calculation of that operator?