CFD Online Discussion Forums

CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   OpenFOAM Programming & Development (http://www.cfd-online.com/Forums/openfoam-programming-development/)
-   -   adding dispersion force term in UEqn leads to oscillations of velocity field (http://www.cfd-online.com/Forums/openfoam-programming-development/118472-adding-dispersion-force-term-ueqn-leads-oscillations-velocity-field.html)

sr_tenedor May 28, 2013 09:46

adding dispersion force term in UEqn leads to oscillations of velocity field
 
1 Attachment(s)
Hello Foamers,

I'm trying to implement an additional force in my solver, to represent dispersion forces between a wall and the fluid particles near the wall. The modelling approach for the magnitude of this force is as follows:

f = -3A / (delta^4)

where delta is the normal distance to the wall.

For my implementation I use a steady state SIMPLE procedure

Code:

  #include "FAdh.H"   

    fvVectorMatrix UEqn
    (
        fvm::div(phi, U)
      - fvm::laplacian(nu, U)
      ==
        //fAdh*vector(0,1,0)
        fvc::reconstruct
        (
            fvc::interpolate(fAdh) * (vector(0.0,1.0,0.0) & mesh.Sf())
        )
    );


    UEqn.relax();
   

    if (momentumPredictor)
    {
        scalar eqnResidual = solve(UEqn == -fvc::grad(p)).initialResidual();
        maxResidual = max(maxResidual, eqnResidual);
    }

Attached are the simulation results of a 1-d case. Herein the bottom boundary represents the wall. As one can see, even if the solver indicates convergence of the results, there remain some oscillations in the velocity field. This leads to an U which is not obeying continuity. I already tried to modify relaxation factors and convergence criteria without success.

Maybe someone of you has a clue about what part of my implementation is causing this physically incorrect velocity field?

Kind regards,
Fabian

sr_tenedor June 6, 2013 09:29

Solution
 
It seems, that I forgot about calculating an additional flux. Now my implementation is based on the modelling of the gravity term in this wiki-entry. In case someone is interested in the final implementation (comments are welcome !):

Code:

    #include "FAdh.H"

    fvVectorMatrix UEqn
    (
        fvm::div(phi, U)
      - fvm::laplacian(nu, U)
    );

    UEqn.relax();

    if (momentumPredictor)
    {
        scalar eqnResidual = solve(UEqn ==
                                    fvc::reconstruct
                                    (
                                        fvc::interpolate(fAdh) * (vector(0.0,1.0,0.0) & mesh.Sf())
                                      - fvc::snGrad(p)*mesh.magSf()
                                    )
                                  ).initialResidual();
        maxResidual = max(maxResidual, eqnResidual);
    }

    p.boundaryField().updateCoeffs();

    volScalarField rUA("rUA", 1.0/UEqn.A());
    surfaceScalarField rUAf("(1|A(U))", fvc::interpolate(rUA));

    U = rUA*UEqn.H();

    surfaceScalarField phiU
    (
        (fvc::interpolate(U) & mesh.Sf())
        +fvc::ddtPhiCorr(rUA, U, phi)
    );

    phi = phiU  + rUAf*fvc::interpolate(fAdh)*(vector(0.0,1.0,0.0) & mesh.Sf());

    // Non-orthogonal pressure corrector loop
    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
    {
        fvScalarMatrix pEqn
        (
            fvm::laplacian(rUAf, p) == fvc::div(phi)
        );

        pEqn.setReference(pRefCell, pRefValue);

        if (nonOrth == 0)
        {
            scalar eqnResidual = pEqn.solve().initialResidual();
            maxResidual = max(maxResidual, eqnResidual);
        }
        else
        {
            pEqn.solve();
        }

        if (nonOrth == nNonOrthCorr)
        {
            phi -= pEqn.flux();
        }
    }

    U += rUA*fvc::reconstruct((phi - phiU)/rUAf);
    U.correctBoundaryConditions();

    #include "continuityErrs.H"


be_inspired June 6, 2014 06:59

In case the force is a vectorField, How could that be implemented?

I am trying to include, OF2.3.0, a vector force calculated over cell centers.

http://www.cfd-online.com/Forums/ope...tml#post495673


All times are GMT -4. The time now is 18:11.