[CGuaq]

Hi everyone,

I try to modify on OpenFOAM 9 the stochasticDispersionRAS.C file to take into account anisotropy near the walls in my simulation (with particles). To do that, I want to add 3 functions for UTurb calculation and so I also need to add a "test" for y* (or y+) in each cells of my geometry to know if I am in near wall region (my criteria is y* < 80 to apply anisotropy correction).

I made the following modifications in the "member functions" part of "stochasticDipsersionRAS.C", the rest of the code is exactly the same:

template<class CloudType>
Foam::vector Foam::myStochasticModel<CloudType>::update
(
    const scalar dt,
    const label celli,
    const vector& U,
    const vector& Uc,
    vector& UTurb,
    scalar& tTurb
)
{
    Random& rnd = this->owner().rndGen();

    const scalar cps = 0.16432;

    const scalar k = this->kPtr_->primitiveField()[celli];
    const scalar epsilon = this->epsilonPtr_->primitiveField()[celli] + rootVSmall;

    const scalar UrelMag = mag(U - Uc - UTurb);

    const scalar tTurbLoc = min(k/epsilon, cps*pow(k, 1.5)/epsilon/(UrelMag + small));

    const objectRegistry& obr = this->owner().mesh();
               
    const word turbName = IOobject::groupName
    (
         momentumTransportModel::typeName,
         this->owner().U().group()
    );
    const momentumTransportModel& turbModel = obr.lookupObject<momentumTransportModel>(turbName);
    const scalarField& y_ = turbModel.y()[celli];

    const scalar uTau = 0.32; //value for test;
    const scalar yStar = uTau*y_[celli]/nuw[celli];

    // Parcel is perturbed by the turbulence
    if (dt < tTurbLoc)
    {
        tTurb += dt;

        if (tTurb > tTurbLoc)
        {            
        	tTurb = 0;
        	const scalar sigma = sqrt(2*k/3.0);
        	// Calculate a random direction dir distributed uniformly
               // in spherical coordinates

             	const scalar theta = rnd.sample01<scalar>()*twoPi;
            	const scalar u = 2*rnd.sample01<scalar>() - 1;
            	const scalar a = sqrt(1 - sqr(u));
            	
            	if (yStar < 80.0)
            	{
            		scalar f_x;
            		scalar f_y;
            		scalar f_z;
            		scalar dir_x;
            		scalar dir_y;
            		scalar dir_z;
            		
            		f_x=1+0.285*(yStar+6)*exp(-0.455*pow(yStar+6, 0.53));
            		f_y=1-exp(-0.02*yStar);
            		f_z=sqrt(3-sqr(f_x)-sqr(f_y));
            	       // update the turbulence fluctuating velocity for particle
            			
            		dir_x=f_x*a*cos(theta);
            		dir_y=f_y*a*sin(theta);
            		dir_z=f_z*u;
            		const vector dir(dir_x, dir_y, dir_z);
            		UTurb=sigma*mag(rnd.scalarNormal())*dir;
            	}
            		
            	else
            	{
            		const vector dir(a*cos(theta), a*sin(theta), u);
            		UTurb = sigma*mag(rnd.scalarNormal())*dir;
       	}
    	    }
    	    
    	}
    
    else
    {
        tTurb = great;
        UTurb = Zero;
    }

    return Uc + UTurb;
}

I can compile without error message but when I launch my simulation (on 4 processors) I have the following error message:

Solving 2-D cloud cloud

[1] --> FOAM FATAL ERROR: 
[1] hanging pointer at index 2095 (size 7), cannot dereference
[1] 
[1]     From function const T& Foam::UPtrList<T>::operator[](Foam::label) const [with T = Foam::fvPatchField<double>; Foam::label = int]
[1]     in file /opt/openfoam9/src/OpenFOAM/lnInclude/UPtrListI.H at line 100.
[1] 
FOAM parallel run aborting
[1] 
[0] #0  [3] #[1] #0  [2] #0  0  Foam::error::printStack(Foam::Ostream&)Foam::error::printStack(Foam::Ostream&)Foam::error::printStack(Foam::Ostream&)Foam::error::printStack(Foam::Ostream&) at ??:?
 at ??:?
 at ??:?
 at ??:?

Please can someone help me to understand this error and why I have it ? I am not sure concerning my call for y (distance between a cell and the near wall) and my ystar calculation.
Also I would like to know if it's necessary to add a "forAll" loop in each cells in my code between "if (tTurb > tTurbLoc)" and "if (yStar < 80.0)" conditions.

Thank you in advance for your help !