|
[Sponsors] |
Help!! customize surface tension term in interFoam |
|
LinkBack | Thread Tools | Search this Thread | Display Modes |
February 12, 2016, 01:15 |
Help!! customize surface tension term in interFoam
|
#1 |
New Member
Join Date: May 2014
Posts: 5
Rep Power: 12 |
Hi guys, I am trying to replace the original surface tension term from interFoam using a custom code.
The original surface tension term mathematically can be expressed as Kappa*sigma*grad(alpha1). And interFoam expresses this term as fvc::interpolate(interface.sigmaK())*fvc::snGrad(a lpha1) Based on my understanding, what interFoam does is to calculate the surface tension on face center and then reconstruct to the cell center. Therefore, "fvc::interpolate(interface.sigmaK())*fvc::snGrad( alpha1)" is a scalar rather than a vector. What I am trying to do is to replace the above mentioned term by another function aiming to calculate the force between any cell to the rest of the cells in the domain. Here is my code added to createFileds.H. It complies and runs OK, but no matter how long I run the code, the results of test case does not make sense. BTW, my test case is to place a air rectangle in a liquid box and see if the rectangle can become a circle eventually. Test case does not make sense means that the air rectangle does not change its shape. Code:
volVectorField newSurfaceTension //create volume vector field ( IOobject ( "newSurfaceTension", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::NO_WRITE ), mesh, dimensionedVector("newSurfaceTension",dimensionSet(1,-4,-2,0,0,0,0),Foam::vector(0,0,0)) ); forAll(mesh.C(),celli) //calculate surface tension on celli { forAll(mesh.C(),cellj) // relation between the celli to the rest of the cells in the domain { if( ((alpha1[celli]<0.5)&& (alpha1[cellj]<0.5)) || ((alpha1[celli]>=0.5)&& (alpha1[cellj]>=0.5)) ) //if both less than 0.5 or both greater than 0.5 { s=constant1; } else // if one is below 0.5 and the other is above 0.5 { s=constant2; } newSurfaceTension[celli]=newSurfaceTension[celli]+(mesh.C()[celli]-mesh.C()[cellj])*s*constant3*mag(mesh.C()[celli]-mesh.C()[cellj]); // my formula } } surfaceVectorField newSurfaceTension1("newSurfaceTension1",fvc::interpolate(newSurfaceTension)); //convert from volVectorField to surfaceVectorField const surfaceVectorField& sf=mesh.Sf(); surfaceScalarField newSurfaceTension2("newSurfaceTension2",newSurfaceTension1 & sf); //covert from surfaceVectorField to surfaceScalarField After that I will update surface tension term in UEqn.H and pEqn.H. Here I use UEqn.H as an example. Code:
fvVectorMatrix UEqn ( fvm::ddt(rho, U) + fvm::div(rhoPhi, U) + turbulence->divDevRhoReff(rho, U) ); UEqn.relax(); if (pimple.momentumPredictor()) { solve ( UEqn == fvc::reconstruct ( ( //fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1) newSurfaceTension2 - ghf*fvc::snGrad(rho) - fvc::snGrad(p_rgh) ) * mesh.magSf() ) ); } |
|
Thread Tools | Search this Thread |
Display Modes | |
|
|
Similar Threads | ||||
Thread | Thread Starter | Forum | Replies | Last Post |
Access the surface tension force at a patch | michielm | OpenFOAM Programming & Development | 1 | December 13, 2013 03:50 |
Surface Tension term in two phase flows | Soghrat | Main CFD Forum | 2 | July 7, 2009 10:20 |
Surface tension | asder | CFX | 2 | August 4, 2008 09:41 |
Surface Tension term | Soghrat | Main CFD Forum | 1 | April 24, 2007 06:45 |
Surface tension - Continuum Surface Force model | Robert | Main CFD Forum | 0 | May 2, 2002 07:34 |