CFD Online Discussion Forums

CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   OpenFOAM (http://www.cfd-online.com/Forums/openfoam/)
-   -   help for a piece of code (http://www.cfd-online.com/Forums/openfoam/105584-help-piece-code.html)

bryant_k August 2, 2012 23:43

change values of a surfaceSalarField at boundary
 
Hello everyone
I want to change values of a surfaceSalarField sphi at boundary 'surface'.To other regions,the values keep default.So I write a piece of code :
Code:

sphi=sigma*(interpolate(U) &mesh.Sf())-fvc::snGrad(ha);

label patchI=mesh.boundaryMesh().findPatchID("surface");
forAll(sphi,patchI)
{
sphi[patchI]=sigma*(interpolate(U) &mesh.Sf());
}


But this is wrong.So can you tell me how to modify the above code?

regards!

bryant

wyldckat August 3, 2012 09:22

Hi Bryant, I got your message.

Look at the file "applications/utilities/postProcessing/wall/yPlusLES/yPlusLES.C" for answers ;)
Also look at "applications/utilities/postProcessing/wall/yPlusRAS/yPlusRAS.C" to see a more basic method of calculating yPlus. ;)

Best regards,
Bruno

bryant_k August 3, 2012 10:51

Hi Bruno
I have looked at the files which you suggested.I am confused with the code inside. I am not familiar with the LES or RAS.
Maybe I have not make my points clearly.Firstly the quantity 'sphi' is a surfaceScalarFiled.And then there are two different expressions for it.
Code:

sphi=sigma*(interpolate(U) &mesh.Sf())//this expression is just for the boundary 'surface'


sphi=sigma*(interpolate(U) &mesh.Sf())-fvc::snGrad(ha);//this expression is for other regions

So can you give me some suggestions in details?

Regards!


bryant

wyldckat August 3, 2012 14:03

If you abstract yourself from whatever LES and RAS are, you would see the details in bold and italic:
Code:

const fvPatchList& patches = mesh.boundary();

        const volScalarField nuLam(sgsModel->nu());

        forAll(patches, patchi)
        {
            const fvPatch& currPatch = patches[patchi];

            if (isA<wallFvPatch>(currPatch))
            {
                yPlus.boundaryField()[patchi] =
                    d[patchi]
                  *sqrt
                    (
                        nuEff.boundaryField()[patchi]
                      *mag(U.boundaryField()[patchi].snGrad())
                    )
                  /nuLam.boundaryField()[patchi];
                const scalarField& Yp = yPlus.boundaryField()[patchi];

                Info<< "Patch " << patchi
                    << " named " << currPatch.name()
                    << " y+ : min: " << gMin(Yp) << " max: " << gMax(Yp)
                    << " average: " << gAverage(Yp) << nl << endl;
            }
        }

It even detects if the patch is a wall...

This should translate to your case to be something like:
Code:

sphi=sigma*(interpolate(U) &mesh.Sf())-fvc::snGrad(ha);

label patchI=mesh.boundaryMesh().findPatchID("surface"); forAll(sphi,patchI)
{
  sphi.boundaryField()[patchI] =
    sigma*(U.boundaryField()[patchi] & mesh.Sf());
}

Ah, wait, now I understand your issue with this... "yPlus" is a "volScalarField"...

OK, the closest I could find so far was "applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/turbulenceModel/wallViscosity.H"...
I found it among the results from using this command:
Code:

find . -name "*.[CH]" | xargs grep -sl surfaceScalarField | sed 's=/[a-zA-Z0-9]*\..*==' | uniq | xargs grep -e "forAll(.*patc" -rsl
Basically it:
  1. Searches for files "*.[CH]";
  2. Then among those, the ones that mention "surfaceScalarFields";
  3. From those get the folders they're in;
  4. Make sure the list of folders is unique;
  5. The find in those folders for files that use "forAll" on patches.


Mmm... this is complicated... Anyone else wants to help?

wyldckat August 3, 2012 14:16

OK, I kept digging a bit deeper... in this file "applications/solvers/electromagnetics/magneticFoam/createFields.H", you'll find these pieces of code:
Code:

    surfaceScalarField murf
    (
//...
    );

    surfaceScalarField Mrf
    (
//...
    );

        const labelList& faces = mesh.faceZones()[magnetZonei];

        const scalar muri = magnets[i].mur();
        const scalar Mri = magnets[i].Mr().value();
        const vector& orientationi = magnets[i].orientation();

        const surfaceVectorField& Sf = mesh.Sf();

        forAll(faces, i)
        {
            label facei = faces[i];
            murf[facei] = muri;
            Mrf[facei] = Mri*(orientationi & Sf[facei]);
        }

Which would imply the need to do face-by-face calculations...

bryant_k August 5, 2012 03:09

Hi Bruno
Thank you for your help.
I have modified my code according the reference file.
Code:

sphi=sigma*(interpolate(U) &mesh.Sf())-fvc::snGrad(ha);
forAll(sphi,i)
{
label zonei = mesh.faceZones().findZoneID("surface");
const labelList& faces = mesh.faceZones()[zoneonei];
    forAll(faces, i)
        {
            label facei = faces[i];   
            sphi[facei] = sigma.value()*(U.boundaryField()[facei]& mesh.Sf[facei]);
        }
}

I can complie it fluently.But it can not run correctly("segments fault").
I have a doubt at this line:
Code:

sphi[facei] = sigma.value()*(U.boundaryField()[facei]& mesh.Sf[facei]);
if I write
Code:

sphi[facei] = sigma*(U.boundaryField()[facei]& mesh.Sf[facei]);
it will failed in compiling.
Do you have any idea about this?

regards!

bryant

wyldckat August 5, 2012 08:08

Hi Bryant,

Something's very wrong here:
Code:

forAll(sphi,i)
{
//...
  forAll(faces, i)

Do not use the same variable "i" for iterating in both loops!

And why is "zonei" being found inside the outer loop and not outside?

Best regards,
Bruno


All times are GMT -4. The time now is 10:35.