[CGuaq]

Hello everyone,

I try to implement CRW model on OpenFOAM 9 based on Langevin equation and I need to add a drift term with gradient of sigma (in normal flow direction).
My recurrence equations for turbulent velocity components are:

equationsCRW.png

This is my code, it's not over yet (in particular for backup variables that need to be defined out of the function for initialization).

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.142857;

    const scalar UrelMag = mag(U - Uc - UTurb);
    
    const objectRegistry& obr = this->owner().mesh();
            
    const word turbName = IOobject::groupName
    (
         momentumTransportModel::typeName,
         this->owner().U().group()
    );
    const momentumTransportModel& turbModel = obr.lookupObject<momentumTransportModel>(turbName);
    
    scalar y_(wallDist::New(this->owner().mesh()).y()[celli]);
    
    const tmp<volScalarField> tk = turbModel.k();
    const volScalarField& k = tk();
    
    const tmp<volScalarField> tepsilon = turbModel.epsilon();
    const volScalarField& epsilon = tepsilon();
    
    const tmp<volScalarField> tnu = turbModel.nu();
    const volScalarField& nu = tnu();
       
    dimensionedScalar UTau("UTau",dimensionSet(0, 1, -1, 0, 0, 0, 0),0.32);
    dimensionedScalar tL("tL", dimensionSet(0, 0, 1, 0, 0, 0, 0), 0);
    dimensionedScalar Tp("Tp", dimensionSet(0, 0, 1, 0, 0, 0, 0), 0);
    dimensionedScalar nu_("nu_", dimensionSet(0, 2, -1, 0, 0, 0, 0), 0.000014607);
    dimensionedScalar k_("k_", dimensionSet(0, 2, -2, 0, 0, 0, 0), k[celli]);
    dimensionedScalar epsilon_("epsilon_", dimensionSet(0, 2, -3, 0, 0, 0, 0), epsilon[celli]);
    dimensionedScalar dp ("dp", dimensionSet(0, 1, 0, 0, 0, 0, 0), 1e-08);
    
    
    const scalar yStar = UTau.value()*y_/nu_.value();   

    
    if (yStar >100)
    {
    	tL=cps*k_/epsilon_;
    }
    else
    {
    	if (yStar > 5)
    	{
    		tL = (7.122+0.5731*yStar-0.00129*pow(yStar, 2))*nu_/pow(UTau, 2);
    	}
    	else
    	{
    		tL = 10*nu_/pow(UTau, 2);
    	}
    }
    	   
    const scalar Te = 2*tL.value();  
    
    const scalar Kn = (2.0*6.6e-08)/dp.value(); 
    const scalar Cc = 1 + Kn*(1.257 + 0.4*exp(-1.1/Kn));
    
    Tp = 2000*pow(dp , 2)*Cc/(18*nu_);
    
    const scalar Le = Te*sqrt(2*k[celli]/3);
    const scalar Dcross = Tp.value()*UrelMag;
    
    scalar tTurbLoc;
    
    if (Le < Dcross)
    {
    	tTurbLoc =-Tp.value()*log(1-(Le/(Tp.value()*(UrelMag+ small))));
    }
    else
    {
    	tTurbLoc = Te;
    }
    
    const scalar Stk=Tp.value()/tL.value();
    
   
/////////////////////backup variables to put out of the function////////////////////////////////////////////////
    scalar sigma_xn=1e-12; 
    scalar sigma_yn=1e-12; 
    scalar sigma_zn=1e-12; 
    	
    scalar sigma_n=1e-12; 
    	
    scalar u_xn=0; 
    scalar u_yn=0; 
    scalar u_zn=0; 
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////


    if (dt < tTurbLoc)
    {
        	scalar lambda1 = mag(rnd.scalarNormal());
            	scalar lambda2 = mag(rnd.scalarNormal());
            	scalar lambda3 = mag(rnd.scalarNormal());
            	
            	const scalar theta = rnd.sample01<scalar>()*twoPi;
            	const scalar u = 2*rnd.sample01<scalar>() - 1;
            	const scalar a = sqrt(1 - sqr(u));
            	
            	if (yStar < 100.0)
            	{
            		          		
            		const scalar sigma_x = (0.4*yStar)/(1+0.0239*pow(yStar, 1.496))*UTau.value();
            		const volScalarField sigma_y = (0.0116*pow(yStar, 2))/(1+0.203*yStar+0.0014*pow(yStar, 2.421))*UTau*nu/nu;
            		const scalar sigma_z = (0.19*yStar)/(1+0.0361*pow(yStar, 1.322))*UTau.value();
            		
            		const volVectorField gradSigmaVol = fvc::grad(sigma_y);
            		
            		const volScalarField gradSigmaY = gradSigmaVol.component(vector::Y);
            		
            		const scalar u_x = (sigma_x/sigma_xn)*u_xn*exp(-dt/tL.value())+sigma_x*sqrt(1-exp(-2*dt/tL.value()))*lambda1;
            		const volScalarField u_y = (sigma_y/sigma_yn)*u_yn*exp(-dt/tL.value())+sigma_y*sqrt(1-exp(-2*dt/tL.value()))*lambda2+(tL*sigma_y)/(1+Stk)*(1-exp(-dt/tL.value()))*gradSigmaY;
            		const scalar u_z = (sigma_z/sigma_zn)*u_zn*exp(-dt/tL.value())+sigma_z*sqrt(1-exp(-2*dt/tL.value()))*lambda3;
               	           			
            		const vector dir(u_x*a*cos(theta), u_y[celli]*a*sin(theta), u_z*u);    		
            		UTurb = dir;
            	
            		u_xn = u_x;
            		u_yn = u_y[celli];
            		u_zn = u_z;
            	
            		sigma_xn = sigma_x;
            		sigma_yn = sigma_y[celli];
            		sigma_zn = sigma_z;    		
            	}
            		
            	else
            	{
            		const scalar sigma = sqrt(2*k[celli]/3.0);
            		const volScalarField sigmaVol = sqrt(2*k/3.0);
            		const volVectorField gradSigmaVol = fvc::grad(sigmaVol);
            		
            		const volScalarField gradSigmaY = gradSigmaVol.component(vector::Y); 
            		
            		
            		const scalar u_x = (sigma/sigma_n)*u_xn*exp(-dt/tL.value())+sigma*sqrt(1-exp(-2*dt/tL.value()))*lambda1;
            		const volScalarField u_y = (sigmaVol/sigma_n)*u_yn*exp(-dt/tL.value())+sigmaVol*sqrt(1-exp(-2*dt/tL.value()))*lambda2+(tL*sigmaVol)/(1+Stk)*(1-exp(-dt/tL.value()))*gradSigmaY;
            		const scalar u_z = (sigma/sigma_n)*u_zn*exp(-dt/tL.value())+sigma*sqrt(1-exp(-2*dt/tL.value()))*lambda3;
            		
               		const vector dir(u_x*a*cos(theta), u_y[celli]*a*sin(theta), u_z*u);    		
            		UTurb = dir;
            	
            		u_xn = u_x;
            		u_yn = u_y[celli];
            		u_zn = u_z;
            	
            		sigma_n = sigma;        		
            		
       	}   
    }
    
    else
    {
        UTurb = Zero;
        u_xn = 0;
        u_yn = 0;
        u_zn = 0;
            	
        sigma_xn = 1e-12;
        sigma_yn = 1e-12;
        sigma_zn = 1e-12;
        
        sigma_n =1e-12;
    }

    return Uc + UTurb;
}

My code works with 1 processor but when I launch in parallel I get the following error message when particles are injected in the flow:
error_message.png

This error disappeared when I delete gradient term, I don't really understand why. Can someone help me to solve this problem ?
Thank you very much.

[CGuaq]

I resolved my problem with gradient, but now I have a last problem to save values of UTurb and sigma components for the next time step.

I created 2 volVectorFiels in the constructor (u_n to save UTurb components values and sigma_n for sigma components) initialized to Zero:

sigma_n
    (
	IOobject
	(
		"sigma_n",
		this->owner_.mesh().time().timeName(),
		this->owner_.mesh(),
		IOobject::NO_READ,
		IOobject::AUTO_WRITE
	),
	this->owner_.mesh(),
	dimensionedVector("sigma_n", dimless, vector(10e-12, 10e-12, 10e-12))
    ),

    u_n
    (
	IOobject
	(
		"u_n",
		this->owner_.mesh().time().timeName(),
		this->owner_.mesh(),
		IOobject::NO_READ,
		IOobject::AUTO_WRITE
	),
	this->owner_.mesh(),
	dimensionedVector("u_n", dimless, vector(0, 0, 0))
    )

Now my problem is to modified this volVectorFields in each cells and for each time steps. I use the following lines to access to values in cells in volVectorFields:

const volVectorField& sigma_n = obr.lookupObject<volVectorField>("sigma_n");
volVectorField& sigma_nMod = const_cast<volVectorField&>(sigma_n);
    
const volVectorField& u_n = obr.lookupObject<volVectorField>("u_n");
volVectorField& u_nMod = const_cast<volVectorField&>(u_n);

And futher I modify the values in each cells with:

u_nMod.primitiveFieldRef()[celli].component(0) = UTurb.component(0);
u_nMod.primitiveFieldRef()[celli].component(1) = UTurb.component(1);
u_nMod.primitiveFieldRef()[celli].component(2) = UTurb.component(2);
            	
sigma_nMod.primitiveFieldRef()[celli].component(0) = sigma_x;
sigma_nMod.primitiveFieldRef()[celli].component(1) = sigma_y;
sigma_nMod.primitiveFieldRef()[celli].component(2) = sigma_z;

I can compile my code and lauch simulation without errors but it doesn't seem to work ...
Did you have another solution to modified volVectorField or volScalarField cell by cell at each time step ? Or maybe help me to find the error in my code ? I'am beginner in programming under OpenFOAM so any idea could be helpful for me !
Thank you very much.

[csimr]

hello,


I have the same problem with VolVectorField. there is someone who has worked with this model before ?

[sinvastil]

Hi Charles,

I wonder if you have solved your problem and if you achieved satisfactory results using the CRW model?

One more question, from which source have you taken the CRW integrated equations shown in the image of the first post?


Thank you.


Regards

[CGuaq]

Hi sinvastil,

The equations comes from the article of Ahmadi & Mofakham intitled "Improved discrete random walk stochastic model for simulating particle dispersion and deposition in inhomogeneous turbulent flow", you can also find some informations in the other articles from them. But this post is old and I modified this equations by considering the results from the article "stochastic modeling of particle diffusion in a turbulent boundary layer". I obtained very good results with this new model for academic cases but also in more complex flows.
Finally,what problem are you interested in?

Regards

[sinvastil]

Hi Charles,

At the moment, I am just exploring the possibility to implement it in OpenFOAM because of the reported limitations of the DRW. Basically, I am looking at the equations, literature, and look for people's opinion about the accuracy of the model for more complex scenarios. So I am glad you got good results for complex flows.


Next, since to my knowledge the code is not available online, I will evaluate the steps I need for the implementation in OpenFOAM. Therefore, it is likely that I will ask some support in the near future.

Thank you.

Regards

[CGuaq]

Hi sinvastil,

Indeed this model is quit difficult to implement depending on your knowledge on OpenFOAM. I will publish an article in the next few weeks on this subject and the application to more complex flows.
After that, I will share my code via gitHub, it is far from perfect and could be improved on certain points. If you need you can use it or participate in its improvement.

[sinvastil]

Hi Charles,

In terms of programming in OpenFOAM I am still a beginner never programmed such model, only small code modications.

That would be very nice, in that way I can hopefully contribute to the improvement of the code instead of writing it from scratch and struggle with the same things that you already did.

I look forward to read your article and test the code. We keep in touch and good luck for the coming weeks.

Regards.