CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Programming & Development

adding dispersion force term in UEqn leads to oscillations of velocity field

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree1Likes
  • 1 Post By sr_tenedor

Reply
 
LinkBack Thread Tools Display Modes
Old   May 28, 2013, 09:46
Default adding dispersion force term in UEqn leads to oscillations of velocity field
  #1
New Member
 
Fabian Gabel
Join Date: Feb 2013
Location: Darmstadt, Germany
Posts: 5
Rep Power: 4
sr_tenedor is on a distinguished road
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
Attached Images
File Type: jpg pU.jpg (21.7 KB, 22 views)
sr_tenedor is offline   Reply With Quote

Old   June 6, 2013, 09:29
Default Solution
  #2
New Member
 
Fabian Gabel
Join Date: Feb 2013
Location: Darmstadt, Germany
Posts: 5
Rep Power: 4
sr_tenedor is on a distinguished road
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"
luiscardona likes this.
sr_tenedor is offline   Reply With Quote

Old   June 6, 2014, 06:59
Default
  #3
Senior Member
 
M. Montero
Join Date: Mar 2009
Location: Madrid
Posts: 112
Rep Power: 8
be_inspired is on a distinguished road
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.

Pressure Jump Correction. Actuator Disc Model and numerical wiggles
be_inspired is offline   Reply With Quote

Reply

Tags
dispersion force ueqn

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
How to add a source term (body force) to icoFoam? chris Liu OpenFOAM Programming & Development 4 October 27, 2014 16:38
Force can not converge colopolo CFX 13 October 4, 2011 22:03
SIMPLE force term isnbt converging what have I missed brooksmoses OpenFOAM Running, Solving & CFD 2 July 17, 2006 20:25
Problems with SUPG body force term FEM question Main CFD Forum 0 January 21, 2006 18:51
Terrible Mistake In Fluid Dynamics History Abhi Main CFD Forum 12 July 8, 2002 09:11


All times are GMT -4. The time now is 03:37.