CFD Online Logo CFD Online URL
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Programming & Development

Translation of equation into OpenFoam

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

LinkBack Thread Tools Search this Thread Display Modes
Old   January 16, 2023, 22:52
Default Translation of equation into OpenFoam
New Member
Michael Coe
Join Date: Jan 2021
Posts: 4
Rep Power: 4
mco is on a distinguished road
Hi everyone,

I've consulted the programming guide for OpenFoam and I'm working on translating an equation for a decay constant to put into a solver. This is a solver for cyclic heat transfer with constant wall temperature.

Suppose I have a field \theta that is a normalized temperature field.

The equation to find the decay constant, \theta, is:

\lambda = - \frac{1}{L} ln\left(1 - \frac{\alpha \int_{wall} \left(\frac{\partial \theta}{\partial n} \right) ds}{\int_{in} \left( u_x \theta + \alpha \frac{\partial \theta}{\partial x}\right)_{in}dy} \right)

here, u_x is the velocity in the direction of flow, n is the normal to the wall, \alpha is the thermal diffusivity of the fluid, and L is the length of the domain.

The numerator represents the heat flux leaving the control volume over the wall surface, and the denominator represents the net heat flux at the inlet boundary condition.

My solution to this is the following code in a custom solver:

// Calculate decay constant lambda

volScalarField alphaEff("alphaEff",turbulence->nu()/Pr + turbulence->nut()/Prt);

vector flowDir = (Ubar.value()/mag(Ubar.value()));
label inletID = mesh.boundaryMesh().findPatchID("inlet");
label wallsID = mesh.boundaryMesh().findPatchID("walls");

if(runTime.value() >= 20)
    // Calculate decay constant lambda Numerator
    dimensionedScalar lambdaNum
        dimensionSet(0, 3, -1, 1, 0, 0, 0),

     // Calculate decay constant lambda Denomenator
    dimensionedScalar lambdaDen
        dimensionSet(0, 3, -1, 1, 0, 0, 0),
        gSum((theta.boundaryField()[inletID]*(flowDir & U.boundaryField()[inletID])*mesh.boundary()[inletID].magSf()) + 
        (alphaEff.boundaryField()[inletID]*(flowDir & fvc::grad(theta)().boundaryField()[inletID])*mesh.boundary()[inletID].magSf()))
    lambda = -(1/Lc)*Foam::log(1 - (lambdaNum/lambdaDen));

Info << "Lambda: " << lambda.value() << endl;

All the constants are defined in the createFields.H file. The solver waits until time steps to let the cyclic pressure field stabilize before calculating \lambda_L

I suspect that it might be something to do with the second term in lambdaDen. Specifically how I take the gradient of the field value at the inlet.

Any help with this would be greatly appreciated.

Best Regards,

mco is offline   Reply With Quote


gradient adaption, periodic heat transfer, solver compilation

Thread Tools Search this Thread
Search this Thread:

Advanced Search
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 Off
Pingbacks are On
Refbacks are On

Similar Threads
Thread Thread Starter Forum Replies Last Post
Porous Modeling of Energy equation in OpenFOAM mohammad_kordo OpenFOAM Running, Solving & CFD 9 November 22, 2020 07:18
How to contribute to the community of OpenFOAM users and to the OpenFOAM technology wyldckat OpenFOAM 17 November 10, 2017 15:54
OpenFOAM v3.0.1 Training, London, Houston, Berlin, Jan-Mar 2016 OpenFOAM Announcements from Other Sources 0 January 5, 2016 03:18
OpenFOAM Training, London, Chicago, Munich, Sep-Oct 2015 OpenFOAM Announcements from Other Sources 2 August 31, 2015 13:36
Solving Advection-reaction Equation using OpenFoam arijith OpenFOAM Programming & Development 1 October 13, 2014 10:02

All times are GMT -4. The time now is 16:09.