CFD Online Discussion Forums - Problems with applying thoms formula

Hello everyone. I'm working on a solver that solves vorticiy-stream function equation. One of the things I wanted to apply, to one of my case, was 0 velocity at the boundary. To achieve that, I wasn't simply applying 0 velocity at the boundary, since I'm not using primitive variables. My idea was to, make such boundary condition that will create vortex itself, that in result of "mixing" will create a 0-velocity slice near the wall.

So this is a code of my Thoms formula, atleast this is what i think its called.

Code:

const fvPatch& boundaryPatch = patch(); // this patch

            // make a reference to the mesh so that I can access variables

            const fvMesh& mesh = patch().boundaryMesh().mesh();

            const volScalarField& psi = db().lookupObject<volScalarField>("psi");

            const volScalarField& OM = db().lookupObject<volScalarField>("OM");

            label inletPatchID = mesh.boundaryMesh().findPatchID("bottom");

            label inletPatchID2 = mesh.boundaryMesh().findPatchID("left");

            // Get reference to boundary field value

            const scalarField& psib = psi.boundaryField()[inletPatchID];  

            const scalarField& OMb = OM.boundaryField()[inletPatchID];

            

            // get coordinate for cell centres

            const fvPatchVectorField& centre = mesh.C().boundaryField()[inletPatchID];

            const fvPatchVectorField& centreleft = mesh.C().boundaryField()[inletPatchID2];

            //const fvPatchVectorField& centre = mesh.C().boundaryField()[inletPatchID];

            

            //area of patch

            scalarField areasc = mesh.magSf().boundaryField()[inletPatchID];

            scalarField OMwall = OMb;

            scalar area = sum(areasc);

            

             scalar d;

         

          // THOMAS FORMULA      

                        forAll (psib, cellI)

                        {

                                d = mag(mesh.C()[cellI].component(vector::Y) - centre[cellI].component(vector::Y));

                                OMwall [cellI] = -2.0*psib[cellI]/(d*d);

                        }

                        Info << (OMwall);



         operator==(OMwall);

Generally all i want to achieve is the vorticity at the wall that will be equal to:

$\omega$ = $\frac{2*\psi}{\Delta y^2}$