[Venky_94]

Hey I'm trying to compile a new viscosity model (the equation is attached in the image).
I'm modifying the cross power law model and changing the calcNu() function as below.

Foam::tmp<Foam::volScalarField>
Foam::viscosityModels::APL::calcNu() const
{
	if (strainRate() <= 0) 
    {
	return pow(pow(mu0_,p_)+pow(K_*pow(scalar(0.00000000001),(n_-1)),p_),pow(p_,-1))/density_;
	}
    else 
    {
	return pow(pow(mu0_,p_)+pow(K_*pow(dimensionedScalar("one", dimTime, 1.0)*strainRate(),(n_-1)),p_),pow(p_,-1))/density_;
	}
}

My aim is to avoid unrealistic values during startup when velocity is zero, which is why I want to implement an if loop where the strain rate is replaced by a very small value during startup after which my actual equation should kick in.

The issue I'm facing is that I don't understand how to implement the if loop for the volVectorField of strainRate and I'm getting the below error. Could someone suggest a workaround or way to do this?

viscosityModels/APL/APL.C:56:19: error: no match for ‘operator<=’ (operand types are ‘Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> >’ and ‘int’)
if (strainRate() <= 0)
~~~~~~~~~~~~~^~~~

[SvenBent]

Hi.

From the error you can see that you are trying compare an entire field to a single integer value.
One way to solve the problem is to use neg0 and pos to differentiate between the negative and positive strain rate.

Foam::viscosityModels::APL::calcNu() const
{
	const volScalarField& sRate = strainRate();

        return(
            neg0(sRate) * (
            pow(pow(mu0_,p_)+pow(K_*pow(scalar(0.00000000001),(n_-1)),p_),pow(p_,-1))/density_
                        ) *
            pos(sRate)* ( 
            pow(pow(mu0_,p_)+pow(K_*pow(dimensionedScalar("one", dimTime, 1.0)*sRate,(n_-1)),p_),pow(p_,-1))/density_
                        )
            );
}

Best regards
Mathias

Originally Posted by Venky_94

Hey I'm trying to compile a new viscosity model (the equation is attached in the image).
I'm modifying the cross power law model and changing the calcNu() function as below.

Foam::tmp<Foam::volScalarField>
Foam::viscosityModels::APL::calcNu() const
{
	if (strainRate() <= 0) 
    {
	return pow(pow(mu0_,p_)+pow(K_*pow(scalar(0.00000000001),(n_-1)),p_),pow(p_,-1))/density_;
	}
    else 
    {
	return pow(pow(mu0_,p_)+pow(K_*pow(dimensionedScalar("one", dimTime, 1.0)*strainRate(),(n_-1)),p_),pow(p_,-1))/density_;
	}
}

My aim is to avoid unrealistic values during startup when velocity is zero, which is why I want to implement an if loop where the strain rate is replaced by a very small value during startup after which my actual equation should kick in.

The issue I'm facing is that I don't understand how to implement the if loop for the volVectorField of strainRate and I'm getting the below error. Could someone suggest a workaround or way to do this?

[Venky_94]

Originally Posted by SvenBent

Hi.

From the error you can see that you are trying compare an entire field to a single integer value.
One way to solve the problem is to use neg0 and pos to differentiate between the negative and positive strain rate.

Foam::viscosityModels::APL::calcNu() const
{
	const volScalarField& sRate = strainRate();

        return(
            neg0(sRate) * (
            pow(pow(mu0_,p_)+pow(K_*pow(scalar(0.00000000001),(n_-1)),p_),pow(p_,-1))/density_
                        ) *
            pos(sRate)* ( 
            pow(pow(mu0_,p_)+pow(K_*pow(dimensionedScalar("one", dimTime, 1.0)*sRate,(n_-1)),p_),pow(p_,-1))/density_
                        )
            );
}

Best regards
Mathias

Thanks Mathias. I tried out the change you suggested and my simulation does not go beyond the '0' time step. However it helped me understand the issue very clearly though. When the initial velocity field is set as uniform(0 0 0), the strain rate is 0 at some points, where the viscosity values blow up to infinity because of the asymptotic part of the equation. I'm able to get the simulation running as expected if the initial velocity field is set to anything but 0 [for example, uniform(0.000000001 0 0)].

I compiled the code as below and identified that the issue is when the strain rate values are 0. All I need right now is a way to replace the strain rate by a very small positive value (say 1.0E-17) in the equation wherever it is 0.

Foam::tmp<Foam::volScalarField>
Foam::viscosityModels::APL::calcNu() const
{
	return pow(pow(mu0_,p_)+pow(K_*pow(dimensionedScalar("one", dimTime, 1.0)*strainRate(),(n_-1)),p_),pow(p_,-1))/density_;
}

I'd be grateful if you could suggest a way to do that.

[SvenBent]

You should be able to use:

return pow(pow(mu0_,p_)+pow(K_*pow(dimensionedScalar("one", dimTime, 1.0)*max(strainRate(),vSmall),(n_-1)),p_),pow(p_,-1))/density_;

[Venky_94]

Originally Posted by SvenBent

You should be able to use:

return pow(pow(mu0_,p_)+pow(K_*pow(dimensionedScalar("one", dimTime, 1.0)*max(strainRate(),vSmall),(n_-1)),p_),pow(p_,-1))/density_;

Thanks a lot. I just had to replace vSmall with a dimensionedscalar with the units of strainrate and it's running perfectly.