CFD Online Discussion Forums - Implementing third order derivatives in OpenFOAM

Hi All,

After successful use of standard solvers of OpenFOAM as well as making small changes to existing solvers, I am now in process of developing a new solver based on the shallowWaterFoam solver from version 1.6 (as a personal research topic). The objective is to include (standard) dispersive Boussinesq Wave (BW) terms before extending the model us such that it can be used to test a newly derived boussinesq approach.

As a starter, I need to add the following two dispersive terms (as presented in Madsen and Sorensen, 1992) to the LHS of the hUEqn as presented in ShallowWaterFoam, leading to the following fvVectorMatrix hUEqn:

fvVectorMatrix hUEqn
(
fvm::ddt(hU)
+ fvm::div(phiv, hU)
-A*fvc::div(fvc::grad(fvc::ddt(hU)))
-B*fvc::grad(fvc::laplacian(h))
);

For sake of stability I did not implement the above, but the below hUEqn, allowing for a time lag for the BW-terms (as is done with turbulence).

fvVectorMatrix hUEqn
(
fvm::ddt(hU)
+ fvm::div(phiv, hU)
+BW_term
);

In which the BW_term is defined after completion of the PISO loop, but within the time loop:
tt == fvc::laplacian(hU);
BW_1 == -A*( fvc::ddt(tt) );
BW_2 == -B*fvc::grad(fvc::laplacian(h));
BW_term == BW_1+BW_2;

This approach does not work due to the fact that I take (i) a third derivative of the variable 'h' {d(h)/dxxx} as well as a second spatial and first time derivative of hU {d(hU)/dxxt}. Variation of interpolation schemes did not help me: I still get a fast growing grassy terms. This is especially the case near the boundaries. Investigations show that taking the laplacian(h) leads to a visually smooth solution but numerically discontinuous solution. This holds also for taking the laplacian from (hU). Taking respectively the grad and ddt of these terms obviously introduces the grassy results.

So I would like to ask the following:
(i) Has anyone encountered this before and how did you solve this problem (I suspect that people dealing with acoustics have the same problem)?
(ii) Is there a way to introduce some kind of smoother on the laplacian results before taking the grad or ddt?
(iii) Would derivation of a dxxx class, following the approach specified in phd-research of Jasak (section 3.3.1) be a solution?
(iv) Am I overlooking something else (probably something obvious, since it is the first time I dive into OpenFOAM like this..)?

Thanks in advance!
Alfons