[FluidFox]

Hi All,

I am new to OpenFOAM, but have experience with CFD.
I am moving over from an axisymmetric 2D code of my own to run 3D simulations in OpenFoam.
I have happily followed a number of tutorials and have now set up my own testcases to run.

What I want to do now is add in a body force.
My test case is a cylindrical cavity, and I want to add a force that will drive rotation of the fluid (e_theta direction in polar coordinates)*.

I have followed a few tutorials and have made a copy of icoFoam as my_icoFoam and successfully run wmake and executed. I have read the source code and it mostly makes sense on a high level, but I am having trouble on where to start to implement a force as a function of position.

Based on other entries I need to modify the code to look something like this:

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

but I am unsure how to define `bodyforce`.

My main question at the moment is how I write a function/variable in the appropriate form to be included in the equation above. From what I have read I can extract cartesian components from the mesh as mesh.C().component(vector::Y), but I am unsure how to use these. Do I need to define an additional field that represents my body force?

My other main question is about defining vectors. Do I have to implement this in cartesian components for OpenFoam, or is it somehow possible to define the body force vectors in cylindrical polar coordinate system.

I have tried to read documentation and other forum posts before asking this, but since I am new to OF I may have missed something or not understood, so apologies if that is the case. Any help or pointing in the direction of a useful tutorial would be much appreciated.Thanks!


*For reference, the simplest force I am looking to implement is a constant times radial distance (k*r) pointing in the direction perpendicular to r and z (clockwise from above).

[SHUBHAM9595]

Hey Jack,

1) As you are solving for a volVectorField (U), you first need to initialize your source term. One of the ways to do it is to define it in createFields.H

Info<< "Reading custom source CS\n" << endl;
volVectorField CS
(
    IOobject
    (
        "CS",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    mesh
);

2) You can define your force term (K*R) in UEqn.H as u already mentioned using the mesh coordinates

volScalarField xx = mesh.C().component(vector::X);
volScalarField yy = mesh.C().component(vector::Y);
volScalarField rr  = sqrt((xx*xx)+(yy*yy));
volScalarField CS1 = K1*rr; // here K1 is const

3) Then u need to assign above volScalarField to the previously initialised volVectorField

CS.replace(vector::X, CS1);
CS.replace(vector::Y, CS1);
CS.replace(vector::Z, CS1);

4) Finally, u have to add CS in ur UEqn.H, now as for the present source term "Kr", its independent of the field variable which we are solving (U), thus u can add it as an EXPLICIT source term. Otherwise u have to add it via "fvm::Sp"

fvVectorMatrix UEqn
        (
            fvm::ddt(U)
          + fvm::div(phi, U)
          - (1/rho)*(CS) 
          - fvm::laplacian(nu, U)
        );

Make sure you use the right dimensions while defining the constant K1, otherwise it'll show u an error while solving.