|
[Sponsors] |
Solution does not converge with temperature-dependent density - chtMultiRegionFOAM |
|
LinkBack | Thread Tools | Search this Thread | Display Modes |
September 17, 2020, 20:53 |
Solution does not converge with temperature-dependent density - chtMultiRegionFOAM
|
#1 |
New Member
Fabrício
Join Date: Mar 2018
Posts: 4
Rep Power: 8 |
Hello, everyone.
I am trying to validate the melting and solidification model by reproducing the experimental procedure found in (Jones et al., 2006). This model was also used by (Muhammad el al., 2015) and it is attached here. To simulate this experiment, I am using OF8 and the solver chtMultiRegion, in which the solid regions are the acrylic and the polycarbonate, whereas the pcm material is the fluid region. As a first approach, I used constant thermal properties in all materials and the solution easily converged. Next, I used variable thermal properties for the fluid region, except the density. Again, the solutions converged well and some results are shown in the second figure. My problem has started when I assumed a polynomial function or the Boussinesq approximation for the density. The values of p_rgh, p and Co explode and the solution displays horrible results in the fluid region after some iterations, despite not crashing. I suspect that this behavior occurs when the temperature begins to change in the fluid region, since it takes some time to diverge. Trying to solve this problem, I tried the following procedures:
None of these procedures have worked. So, I think that I might be setting the B.C in the fluid zone or using the pRef incorrectly. However, at this point, I am out of solutions in mind and if someone could help me, I would appreciate. I am attaching the temperature, velocity and p_rgh fields after the solution diverges. Best regards, #Code: Code:
PIMPLE: Iteration 100 Solving for fluid region pcm DILUPBiCGStab: Solving for Ux, Initial residual = 0.025274, Final residual = 9.00165e-05, No Iterations 2 DILUPBiCGStab: Solving for Uy, Initial residual = 0.0200058, Final residual = 8.38601e-05, No Iterations 2 DILUPBiCGStab: Solving for Uz, Initial residual = 0.0911808, Final residual = 0.000685502, No Iterations 1 DILUPBiCGStab: Solving for e, Initial residual = 1, Final residual = 0.00887018, No Iterations 2 Min/max T:296.15 343.15 DICPCG: Solving for p_rgh, Initial residual = 0.43215, Final residual = 0.00425036, No Iterations 76 DICPCG: Solving for p_rgh, Initial residual = 0.00382949, Final residual = 7.87853e-05, No Iterations 77 pressureControl: p min -7.88183e+28 diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0 time step continuity errors : sum local = 11326.3, global = 0.0218294, cumulative = 0.0218294 DICPCG: Solving for p_rgh, Initial residual = 0.108544, Final residual = 0.00107419, No Iterations 64 DICPCG: Solving for p_rgh, Initial residual = 0.00104525, Final residual = 9.27925e-05, No Iterations 65 pressureControl: p min -7.88183e+28 diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0 time step continuity errors : sum local = 11884.1, global = 0.0218294, cumulative = 0.0436589 Solving for solid region acrylic DICPCG: Solving for e, Initial residual = 0.0190569, Final residual = 1.83692e-06, No Iterations 1 DICPCG: Solving for e, Initial residual = 1.8329e-06, Final residual = 3.44873e-10, No Iterations 1 Min/max T:296.15 328.937 Solving for solid region polycarbonate DICPCG: Solving for e, Initial residual = 0.00971627, Final residual = 1.11858e-06, No Iterations 1 DICPCG: Solving for e, Initial residual = 1.11916e-06, Final residual = 2.52876e-10, No Iterations 1 Min/max T:296.15 343.15 PIMPLE: Not converged within 100 iterations ExecutionTime = 26.66 s ClockTime = 26 s Region: pcm Courant Number mean: 1.53665e+10 max: 5.8385e+10 Region: acrylic Diffusion Number mean: 0.0100368 max: 11.2252 Region: polycarbonate Diffusion Number mean: 0.0110556 max: 0.0464376 Time = 0.07 0/pcm/T Code:
/*--------------------------------*- C++ -*----------------------------------*\ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org \\ / A nd | Version: 7 \\/ M anipulation | \*---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class volScalarField; location "0/heater"; object T; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // dimensions [ 0 0 0 1 0 0 0 ]; internalField uniform 296.15; boundaryField { top { type zeroGradient; } front { type wedge; } back { type wedge; } pcm_to_acrylic { type compressible::turbulentTemperatureCoupledBaffleMixed; value uniform 296.15; kappaMethod solidThermo; Tnbr T; } pcm_to_polycarbonate { type compressible::turbulentTemperatureCoupledBaffleMixed; value uniform 296.15; kappaMethod solidThermo; Tnbr T; } Code:
/*--------------------------------*- C++ -*----------------------------------*\ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org \\ / A nd | Version: 7 \\/ M anipulation | \*---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class volScalarField; location "0/heater"; object p; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // dimensions [1 -1 -2 0 0 0 0]; internalField uniform 1e5; boundaryField { top { type calculated; value $internalField; } front { type wedge; } back { type wedge; } pcm_to_acrylic { type calculated; value $internalField; } pcm_to_polycarbonate { type calculated; value $internalField; } Code:
========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org \\ / A nd | Version: 8 \\/ M anipulation | \*---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class volScalarField; location "0/fluid"; object p_rgh; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // dimensions [ 1 -1 -2 0 0 0 0 ]; internalField uniform 99500; boundaryField { top { type fixedFluxPressure; value $internalField; } front { type wedge; } back { type wedge; } pcm_to_acrylic { type fixedFluxPressure; value $internalField; } pcm_to_polycarbonate { type fixedFluxPressure; value $internalField; } } // ************************************************************************* // Code:
/*--------------------------------*- C++ -*----------------------------------*\ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org \\ / A nd | Version: 8 \\/ M anipulation | \*---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class volVectorField; location "0/fluid"; object U; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // dimensions [ 0 1 -1 0 0 0 0 ]; internalField uniform (0 0 0); boundaryField { top { type noSlip; } front { type wedge; } back { type wedge; } pcm_to_acrylic { type noSlip; } pcm_to_polycarbonate { type noSlip; } } // ************************************************************************* // Code:
/*--------------------------------*- C++ -*----------------------------------*\ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org \\ / A nd | Version: 7 \\/ M anipulation | \*---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class dictionary; location "constant/steel"; object thermophysicalProperties; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // thermoType { type heRhoThermo; mixture pureMixture; transport polynomial; thermo hPolynomial; equationOfState icoPolynomial; specie specie; energy sensibleInternalEnergy; } dpdt off; mixture { // PCM specie { molWeight 50; } equationOfState { rhoCoeffs<8> (2456.56 -5.22 0 0 0 0 0 0 ); } thermodynamics { Hf 0; Sf 0; CpCoeffs<8> (-3726.924 18.96 0 0 0 0 0 0 ); } transport { muCoeffs<8> (0.11994 -0.0006529 90 0 0 0 0 0); kappaCoeffs<8> (3.726502 -1.108e-2 0 0 0 0 0 0); } } // ************************************************************************* // Code:
/*---------------------------------------------------------------------------*\ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org \\ / A nd | Version: 8 \\/ M anipulation | \*---------------------------------------------------------------------------*/ Build : 8-9b73cf21a682 Exec : checkMesh Date : Sep 17 2020 Time : 20:33:23 Host : "fabricio-A320M-S2H" PID : 25153 I/O : uncollated Case : /home/fabricio/Desktop/OF_cases/Phase_Change/Validation nProcs : 1 sigFpe : Enabling floating point exception trapping (FOAM_SIGFPE). fileModificationChecking : Monitoring run-time modified files using timeStampMaster (fileModificationSkew 10) allowSystemOperations : Allowing user-supplied system call operations // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // Create time Create polyMesh for time = 0 Time = 0 Mesh stats points: 40401 internal points: 0 faces: 80100 internal faces: 39700 cells: 20000 faces per cell: 5.99 boundary patches: 6 point zones: 0 face zones: 0 cell zones: 3 Overall number of cells of each type: hexahedra: 19800 prisms: 200 wedges: 0 pyramids: 0 tet wedges: 0 tetrahedra: 0 polyhedra: 0 Checking topology... Boundary definition OK. Cell to face addressing OK. Point usage OK. Upper triangular ordering OK. Face vertices OK. Number of regions: 1 (OK). Checking patch topology for multiply connected surfaces... Patch Faces Points Surface topology left 0 0 ok (empty) top 100 201 ok (non-closed singly connected) bottom 100 201 ok (non-closed singly connected) right 200 402 ok (non-closed singly connected) back 20000 20301 ok (non-closed singly connected) front 20000 20301 ok (non-closed singly connected) Checking geometry... Overall domain bounding box (0 0 -0.00151577) (0.0347169 0.0658 0.00151577) Mesh has 2 geometric (non-empty/wedge) directions (1 1 0) Mesh has 3 solution (non-empty) directions (1 1 1) Wedge back with angle 2.5 degrees Wedge front with angle 2.5 degrees All edges aligned with or perpendicular to non-empty directions. Boundary openness (1.71133e-14 5.13242e-18 -2.47923e-16) OK. Max cell openness = 3.18472e-16 OK. Max aspect ratio = 1.89533 OK. Minimum face area = 5.2623e-09. Maximum face area = 1.0472e-06. Face area magnitudes OK. Min volume = 1.7313e-12. Max volume = 3.44528e-10. Total volume = 3.46259e-06. Cell volumes OK. Mesh non-orthogonality Max: 0 average: 0 Non-orthogonality check OK. Face pyramids OK. Max skewness = 0.330796 OK. Coupled point location match (average 0) OK. Mesh OK. End Code:
/*--------------------------------*- C++ -*----------------------------------*\ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org \\ / A nd | Version: 7 \\/ M anipulation | \*---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class dictionary; location "system"; object controlDict; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // application chtMultiRegionFoam; startFrom latestTime; startTime 0; stopAt endTime; endTime 4000; deltaT 0.01; writeControl timeStep; writeInterval 1; purgeWrite 0; writeFormat binary; writePrecision 6; writeCompression off; timeFormat general; timePrecision 6; runTimeModifiable true; adjustTimeStep no; maxCo 1; maxDeltaT 0.01; // ************************************************************************* // Code:
*--------------------------------*- C++ -*----------------------------------*\ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | | \\ / O peration | Version: 3.0.0 | | \\ / A nd | Web: www.OpenFOAM.org | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class dictionary; location "system"; object fvSchemes; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // ddtSchemes { default Euler; } gradSchemes { default Gauss linear; } divSchemes { default none; div(phi,U) bounded Gauss linearUpwind grad(U); div(phi,e) bounded Gauss linearUpwind default; div(phi,R) bounded Gauss linearUpwind default; div(phi,K) Gauss linear; div(phi,Ekp) Gauss linear; div(R) Gauss linear; div(phiv,p) Gauss linear; div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear; } laplacianSchemes { default Gauss linear corrected; } interpolationSchemes { default linear; } snGradSchemes { default corrected; } // ************************************************************************* // Code:
/*--------------------------------*- C++ -*----------------------------------*\ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | | \\ / O peration | Version: 3.0.0 | | \\ / A nd | Web: www.OpenFOAM.org | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class dictionary; location "system"; object fvSolution; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // solvers { rho { solver PCG; preconditioner DIC; tolerance 1e-6; relTol 0.1; } rhoFinal { $rho; tolerance 1e-6; relTol 0.1; } p_rgh { solver PCG; preconditioner DIC; tolerance 1e-4; relTol 0.01; maxIter 5000; } p_rghFinal { $p_rgh; relTol 0; } "(U|h|e).*" { solver PBiCGStab; preconditioner DILU; tolerance 1e-7; relTol 0.01; maxIter 5000; } } PIMPLE { momentumPredictor yes; nOuterCorrectors 100; nCorrectors 2; nNonOrthogonalCorrectors 1; pRefCell 0; pRefValue 0; outerCorrectorResidualControl { p_rgh { relTol 0; tolerance 0.01; } } } relaxationFactors { fields { "rho.*" 1; "p_rgh.*" 0.7; } equations { "U.*" 0.7; "(h|e)" 0.7; } } // ************************************************************************* // Jones, Benjamin J., et al. "Experimental and numerical study of melting in a cylinder." International Journal of Heat and Mass Transfer 49.15-16 (2006): 2724-2738. Muhammad, M. D., O. Badr, and H. Yeung. "Validation of a CFD melting and solidification model for phase change in vertical cylinders." Numerical Heat Transfer, Part A: Applications 68.5 (2015): 501-511. |
|
September 20, 2020, 11:45 |
|
#2 |
Senior Member
Joachim Herb
Join Date: Sep 2010
Posts: 650
Rep Power: 22 |
What I have done once, when the polynomial density function resulted in crashs, was to fit or better mis-use the janaf and sutherland methods for cp and mu (for the temperature/pressure range of interest) and used Boussinesq for the equation of state:
Code:
thermoType { type heRhoThermo; mixture pureMixture; transport sutherland; thermo janaf; equationOfState Boussinesq; specie specie; energy sensibleInternalEnergy; } Another way, I got it to work, was to create a new equation of state mode for the liquid including the compressibility. So something like: Code:
makeThermos ( rhoThermo, heRhoThermo, pureMixture, polynomialTransport, sensibleInternalEnergy, hPolynomialThermo, icoPolCompress, specie ); Code:
template<class Specie, int PolySize> class icoPolCompress; Code:
//- Is the equation of state is incompressible i.e. rho != f(p) static const bool incompressible = true; Code:
template<class Specie, int PolySize> inline Foam::scalar Foam::icoPolCompress<Specie, PolySize>::rho ( scalar p, scalar T ) const { scalar rho0 = rhoCoeffs_.value(T); // calc compressibility return rho0*(1 + (p - 1e5)*compr); // apply it to density } |
|
September 20, 2020, 15:48 |
|
#3 |
Senior Member
Join Date: Sep 2013
Posts: 353
Rep Power: 21 |
Your settings seem well suited for the problem. Are you certain that is it not a problem of your polynomials getting out of bound? If for some reason the temperature drops or increases in such a fashion that one of your values becomes negative your simulation is bound to blow up. Might want to add an fvOption to limit the temperature between those ranges, for numerical reasons. From a quick glance your temperature range in which all of your polynomials are positive is quite small.
Why the slight difference in initial values between p and p_rgh? maxCo=1 seems a bit high and you did not specify maxDi. Although this does not seem diffusion bound. Nevertheless maybe 0.5 for both, especially with Euler time stepping. The rest is your divSchemes. Does it work with upwind schemes? You are using bounded with your schemes. This is usually only done with steadyState cases. This shouldn't be there, but might be ok due to the outerCorrectors.... Since your mesh seems perfectly orthogonal, and only slightly skewed you could even consider switching your laplacian/snGradSchemes to uncorrected or orthogonal. Leads to a slightly increased stability and should be fine on such a nice grid. In fvSolution consider setting the relTol to zero for the Final parameters. |
|
September 20, 2020, 20:21 |
|
#4 |
Senior Member
Peter Hess
Join Date: Apr 2011
Location: Austria
Posts: 250
Rep Power: 17 |
What about the gravity?
There is no inlet or oulet bc in the velocity field. That why I guess, you have free convection. The prgh field shows no gravity effect. The temperature shows no gravity effect also. Post the gravity field please. Regards Peter |
|
January 20, 2021, 10:29 |
|
#5 |
New Member
Fabrício
Join Date: Mar 2018
Posts: 4
Rep Power: 8 |
Sorry for the late reply and thanks for helping me.
So, the first mistake that I noticed was that the Boussinesq approximation is already implemented in the source term related to the melting/solidification process. solidificationMeltingSource.C Code:
forAll(cells_, i) { const label celli = cells_[i]; const scalar Vc = V[celli]; const scalar alpha1c = alpha1_[celli]; const scalar S = -Cu_*sqr(1.0 - alpha1c)/(pow3(alpha1c) + q_); const vector Sb = rhoRef_*g*beta_*deltaT_[i]; Sp[celli] += Vc*S; Su[celli] += Vc*Sb; } Code:
thermoType { type heRhoThermo; mixture pureMixture; transport polynomial; thermo hPolynomial; equationOfState icoPolynomial; specie specie; energy sensibleInternalEnergy; } dpdt off; mixture { // PCM specie { molWeight 50; } equationOfState { rhoCoeffs<8> (769 0 0 0 0 0 0 0); } thermodynamics { Hf 0; Sf 0; CpCoeffs<8> (-3726.924 18.96 0 0 0 0 0 0 ); } transport { muCoeffs<8> (0.11994 -0.0006529 9e-7 0 0 0 0 0); kappaCoeffs<8> (2.16893 -5.89361e-3 0 0 0 0 0 0); } } The remaining mistakes were related to some polynomials being out of bound, as stated by Bloerb. I am sharing some files and results just for info, thanks again! fvOptions Code:
options { solidificationMeltingSource1 { type solidificationMeltingSource; active yes; selectionMode all; Tsol 309.55; Tliq 309.65; L 248000; thermoMode thermo; beta 8.161e-4; rhoRef 769; Cu 1.6e+06; q 1.0e-03; } limitT { type limitTemperature; active yes; selectionMode all; min 296.15; max 343.15; } } // ************************************************************************* // Code:
ddtSchemes { default Euler; } gradSchemes { default leastSquares; } divSchemes { default none; div(phi,U) Gauss linearUpwind grad(U); div(phi,e) Gauss linearUpwind default; div(phi,K) Gauss linearUpwind default; div(phiv,p) Gauss linear; div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear; } laplacianSchemes { default Gauss linear corrected; } interpolationSchemes { default linear; } snGradSchemes { default corrected; } |
|
September 8, 2024, 15:14 |
Thermal shrinkage of PCM during solidification
|
#6 |
New Member
JANGA RAKESH KUMAR
Join Date: Aug 2024
Posts: 13
Rep Power: 2 |
Hello Foamers,
I'm also working on a similar problem, but also involves air phase at top of liquid initially, and also densities different in liquid and solid of PCM, so that a shrinkage void forms after solidification. I have tried to setup this using interFoam but not handling, also tried to couple solidification fvOption with interFoam but I ran into errors. Anyone know how to setup this problem??? If anyone has any idea of, how to setup this problem, please explain. Any hep regarding this problem is greatly appreciated. Thanks in advance. -Rakesh |
|
Tags |
chtmultiregionfoam, icopolynomial |
Thread Tools | Search this Thread |
Display Modes | |
|
|
Similar Threads | ||||
Thread | Thread Starter | Forum | Replies | Last Post |
Getting density, temperature, and other cell flow variables via macros in Fluent UDF | ingabobjoe | Fluent UDF and Scheme Programming | 2 | August 15, 2020 22:08 |
Continuity Error for Temperature Dependent Density | TomasDenk | OpenFOAM Running, Solving & CFD | 5 | January 30, 2020 19:53 |
temperature dependent density of water in fluent | sahar.mh | Fluent UDF and Scheme Programming | 0 | November 15, 2019 11:00 |
UDF to Define Temperature Dependent Negative Heat Source | ATIKADAR | Fluent UDF and Scheme Programming | 1 | September 23, 2019 04:52 |
Calculation of the Governing Equations | Mihail | CFX | 7 | September 7, 2014 07:27 |