[anuargimenez]

Hi all,

I am trying to simulate a polymer flow through a single screw extruder, such as the one shown below:


I modified some libraries to correctly introduce the non-Newtonian fluid behavior and ran some simulations with simple geometries to check that everything worked fine. They converge very well and have physical meaning.

Then, I changed to a more complex geometry, as the one shown in the attachment. It is a static single screw extruder, i am not considering the rotation of the screw.
ExtruderGeometry.png

The problem is that I cannot achieve convergence. The enthalpy is going well, however the velocities and pressure remain constant between e-3 and e-4. An image of residuals is attached.
Residuals.jpg

It is a stationary simulation, using rhoSimpleFoam.

My initial conditions are:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2112                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    object      p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [1 -1 -2 0 0 0 0];

internalField   uniform 1e5;

boundaryField
{    
    wall1
    {
        type            fixedFluxPressure;
        value           uniform 0;
    }
    
    wall2
    {
        type            fixedFluxPressure;
        value           uniform 0;
    }
    
    screw
    {
        type            fixedFluxPressure;
        value           uniform 0;  
    }
    
    outlet
    {
        type            fixedValue;
        value           uniform 1e5;
    }
    inlet
    {
        type            zeroGradient;
    }
}


// ************************************************************************* //

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2112                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    object      T;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 0 0 1 0 0 0];

internalField   uniform 405.65; //initial condition

boundaryField
{   
    wall1
    {
	type            	convectiveRadiativeHeatFlux;
	mode            	convectiveRadiative;
	Ta              	constant 513.15;
	h               	uniform 215;			//Heat transfer coefficient    [W/(m²*K)]
	emissivityNozzle	0.95;				//Emissivity nozzle (or brass) [-]
	emissivityFilament	0.95;				//Emissivity filament [-]
	dNozzle		4.0;				//inner diameter nozzle [mm}
	dFilament		1.75;				//filament diameter [mm]
        kappaMethod     	fluidThermo;
        value           	$internalField;		//initial Temperature value
    }

    wall2
    {
        type            fixedValue;
        value           uniform 513.15; // solo calor en la zona 2
    }
    
    screw
    {
        type            zeroGradient;
    }
    
    outlet
    {
        type            zeroGradient;
    }

    inlet
    {
        type            fixedValue;
        value           uniform 298.15; // temperature ambient 25ºC
    }
}


// ************************************************************************* //

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2112                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volVectorField;
    object      U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 1 -1 0 0 0 0];

internalField   uniform (0 0 0);

boundaryField
{  
    wall1
    {
        type            slip;
    }
    
    wall2
    {
        type            noSlip;
    }
    
    screw
    {
        type            slip;
    }
    
    outlet
    {
        type            zeroGradient;
    }

    inlet
    {
        type            flowRateInletVelocity;
        massFlowRate    constant 1e-05; // kg/s
        rho             rho;
        rhoInlet        1000;    // Guess for rho 

    }
}


// ************************************************************************* //

controlDict:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2112                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

application     polymSimpleFoam;

startFrom       latestTime;

startTime       0;

stopAt          endTime;

endTime         0.003;

deltaT          0.00000005; //Time step should be chosen so small that a Courant number is less than 1

writeControl    timeStep;

writeInterval   600; //Number of output folders or steps in ParaView (deltaT*writeInterval = step width of the folders)

purgeWrite      0;

writeFormat     binary;

writePrecision  10;

writeCompression off;

timeFormat      general;

timePrecision   6;

graphFormat     raw;

runTimeModifiable true;

//_______below additional outputs and options
functions //https://www.openfoam.com/documentation/guides/latest/doc/guide-function-objects.html
{

	// Writes the residuals for each iteration step to a file in the postProcessing folder, graphing while solving via the "AllMonitor" script
	residual1 
	{
		type		residuals;
		libs		("libutilityFunctionObjects.so");
		writeControl    timeStep;
		writeInterval   1;
		fields (p U h);
	}

	// The fieldMinMax function object computes the values and locations of field minima and maxima. These are good indicators of calculation performance, e.g. to confirm that predicted results are within expected bounds, or how well a case is converging.
	minMax
	{
		type            volFieldValue;
		libs            ("libfieldFunctionObjects.so");
		fields          (U p);
		executeControl  writeTime;
		writeControl    writeTime;

		writeFields	false;
		log		false;

		regionType	all;
		operation	max;
	}

	writeThermophysicalProperties
	{
		type		writeObjects;
		libs		("libutilityFunctionObjects.so");
		objects	(thermo:alpha thermo:mu);
		//executeControl	timeStep;
		writeControl	writeTime;
		//writeInterval	60;
	}
}

// ************************************************************************* //

fvSchemes:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2112                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

ddtSchemes
{
    default             steadyState;
}

gradSchemes
{
    default         Gauss linear upwind;

    limited         cellLimited Gauss linear 0.1;
    grad(U)         $limited;
}

divSchemes
{
    default         none;

    div(phi,U)      bounded Gauss linearUpwind grad(U);

    energy          bounded Gauss linearUpwind limited;

    div(phi,e)      bounded Gauss linearUpwind limited;
    div(phi,Ekp)    bounded Gauss linearUpwind limited;
    div(phi,h)      bounded Gauss linearUpwind limited;
    div(phi,K)      bounded Gauss linearUpwind limited;

    div(((rho*nuEff)*dev2(T(grad(U)))))    Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         corrected;
}


// ************************************************************************* //

fvSolution:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2112                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

solvers
{
    p
    {
        solver          	GAMG;
        preconditioner  	GAMG;
        tolerance       	1e-07;
        relTol          	0.1;
	smoother		DICGaussSeidel;
	cacheAgglomeration	true;
	nCellsInCoarsestLevel	145;
	agglomerator		faceAreaPair;
	mergeLevels		1;
	directSolveCoarsest	false;
    }

    U
    {
        solver         smoothSolver;
	smoother	GaussSeidel;
        nSweeps		1;
        tolerance       1e-08;
        relTol          0.01;
    }

    "(e|h)"
    {
        solver         smoothSolver;
	smoother	GaussSeidel;
        nSweeps		2;
        tolerance       1e-09;
        relTol          0.01;
    }
}

SIMPLE
{
    nNonOrthogonalCorrectors 2;
    pMinFactor      0.99;
    pMaxFactor      200;

    transonic       no;
    consistent      yes;

    rhoMax          rhoMax [ 1 -3 0 0 0 0 0 ] 1.2;
    rhoMin          rhoMin [ 1 -3 0 0 0 0 0 ] 0.9;
    
    residualControl
    {
        p               1e-6;
        U               1e-5;
        h		 1e-7;
    }
}

relaxationFactors
{
    U               0.6;
    p               0.5;
    rho             0.01;
    h               0.5;
}


// ************************************************************************* //

I already try to refine the mesh, but the residuals remain almost the same.
It would be a great help if you could give me some advice to improve convergence.

Thanks,
Anuar

[Pappelau]

Hello Anuar,
from shortly locking on your case :

Your control dict "looks" like a transient simualtion but u defined steadystate in fvSchemes....
this is pure optics


From your resiudal plot your case is converging very well.. they not need to be very low for good convergence. More important are your physical parameters like temperature distribution and flowfield speed. Never only watch residuals they just show you the direction if something is odd.



But...

IF you use consistent setting with simple (SIMPLEC) your under relaxation should be near 1 for pressure and rho and 0.95 for other convectives (https://www.openfoam.com/documentati...orithm-control) Chapter 6.3.2

[anuargimenez]

Originally Posted by Pappelau

Hello Anuar,
from shortly locking on your case :

Your control dict "looks" like a transient simualtion but u defined steadystate in fvSchemes....
this is pure optics


From your resiudal plot your case is converging very well.. they not need to be very low for good convergence. More important are your physical parameters like temperature distribution and flowfield speed. Never only watch residuals they just show you the direction if something is odd.



But...

IF you use consistent setting with simple (SIMPLEC) your under relaxation should be near 1 for pressure and rho and 0.95 for other convectives (https://www.openfoam.com/documentati...orithm-control) Chapter 6.3.2

Hi Pappelau,

thank you so much for your quick response.

I read somewhere on the internet, that to improve convergence it is convenient to decrease the number of Courant, so I set a very small timestep. However, as you say, since it is a stationary simulation, without time derivative, it is useless. I have corrected this by setting dt = 1.

I have no experience with the SIMPLEC algorithm, your answer was very helpful. I have changed the under relaxation coefficients to:

relaxationFactors
{
    U               0.9;
    p               1;
    rho             1;
    h               0.9;
}

After these changes, here is the residuals plot:
Residualsv2.jpg

It behaves a little unstable at the beginning and then converges with a little high residual. I have to mention that small changes in the relaxationFactors values highly affect the residuals plot. I don't know if I could consider this a good convergence. Maybe a change in fvSchemes would help to improve this.

Regarding the physical meaning of the simulation with these parameters, the temperature seems ok. However, the velocity (and therefore pressure) do not make me feel confident (which does not mean that they are wrong).

The boundary condition for inlet velocity is as follows:

    inlet
    {
        type                flowRateInletVelocity;
        massFlowRate        3e-05;// kg/s
        extrapolateProfile  yes;
        rho                 rho;
        rhoInlet            1040;
        value               uniform (0 0 0);    
    }

And the posprocessing is shown below:

Temperature:
Post_T.png

Pressure:
Post_p.png

Velocity:
Post_U.png

I would greatly appreciate your opinion.
Thanks,
Anuar

[Pappelau]

Good Morning,
from my Point of View the resiudals are just fine, they show a nice convergence .
the Pictures for Temperature seems very plausible also for velocity. The pressure goes high due to expansion(heating)(of your polymer), I assume ? (Never worked with that specific solver nor the underlying physics)

As you mentioned your fvSchemes u can always take a look at:

https://cfd.direct/openfoam/user-guide/v9-fvschemes/

Here u find good information on "practical" values