# Addition of body force as function of position in icoFoam solver.

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

 July 14, 2020, 05:23 Addition of body force as function of position in icoFoam solver. #1 New Member   Jack Join Date: Jul 2016 Posts: 8 Rep Power: 10 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: Code: ``` 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).

 September 4, 2020, 04:18 #2 Member   MNM Join Date: Aug 2017 Posts: 69 Rep Power: 8 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 Code: ```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 Code: ```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 Code: ```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" Code: ```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. Divyaprakash, mllokloks and Tibo99 like this.

 Tags body force, coordinates, icofoam problem, source code