CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Running, Solving & CFD

Temperature convergence problem with nonsmooth thermal conductivity

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   November 6, 2017, 18:10
Default Temperature convergence problem with nonsmooth thermal conductivity
  #1
New Member
 
Join Date: Jan 2017
Posts: 7
Rep Power: 8
bmercer is on a distinguished road
Hi all,

I am trying to run a problem where thermal conductivity is a non-smooth function of temperature, where the non-smoothness results from dealing with a phase change material. I am having trouble getting the temperature equation of the solver to converge under certain conditions and would be curious if anyone has any thoughts as to what may be going wrong.

In my model, the conductivity is a linear function of temperature. For regions yet to be melted, the conductivity is k1(T) = a*T + b. Once a region has been melted, the conductivity becomes k2 = kfactor*k1(T), even if it re-solidifies. kfactor is between 10 and 100 for problems I want to run. This creates regions in the mesh where cells located next to each other can have similar temperatures but have a jump in conductivity between 10x and 100x. When I set kfactor to 1, my case runs smoothly with no convergence issues. If I modify my code to ignore the melted/unmelted distinction and let conductivity k = k2(T), I get the same (good) behavior as when I use kfactor = 1. This leads me to believe that I donít have a bug in the implementation of the conductivity function and that the non-smooth jump in conductivity between cells may be the culprit.

Below I include my solver setting files, as well as the last few lines of the OpenFoam equation residual output which includes the initial temperature equation blow-up and subsequent simulation crash. Iíd be grateful if anyone has any ideas on how to resolve the issue of the temperature equation blowing up for kfactor > 10. Iím hopeful that tweaking some solver settings could iron out the convergence issues Ė Iíve made some blind attempts so far, but no luck. Any thoughts are appreciated! Thank you!

fvSchemes:
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         Gauss upwind;
    div(phi,U)      Gauss upwind;
    div((phi*interpolate(cp)),T) Gauss upwind;
}

laplacianSchemes
{
    default                   none; //Gauss linear corrected; 
    laplacian((lambda|rho),T) Gauss linear uncorrected; //corrected 
    laplacian(nu,U)           Gauss linear uncorrected; // corrected
    laplacian((1|A(U)),p_rgh)     Gauss linear uncorrected; // corrected
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         uncorrected; //corrected;
}

fluxRequired
{
  default no;
  p_rgh;
}


// ************************************************************************* //
fvSolution:
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
{

    p_rgh
    {
        solver          GAMG;
        smoother        GaussSeidel;
        tolerance       1e-8;
        relTol          0.01;
    }

    p_rghFinal
    {
        $p_rgh;
        relTol          0;
    }

    T
    {
        solver          PBiCG;
        preconditioner   DILU;
        tolerance        1e-7;
        relTol            0.1;
    }

    TFinal
    {
        solver          PBiCG;
        preconditioner   DILU;
        tolerance        1e-9;
        relTol              0;
    }

    U
    {
        solver          PBiCG;
        preconditioner  DILU;
        tolerance       1e-6;
        relTol          0.1;
    }

    UFinal
    {
        $U;
        relTol          0;
    }
}

PIMPLE
{
    momentumPredictor yes; //yes;
    nOuterCorrectors 1; //1;
    nCorrectors     1;
    nNonOrthogonalCorrectors 0;
    pRefCell 0;
    pRefValue 0;
}


// ************************************************************************* //
OpenFoam solver output:
Code:
...
...
...
Courant Number mean: 0.0006917597 max: 0.27723292
deltaT = 9.9667774e-07
PIMPLE: iteration 1
DILUPBiCG:  Solving for Ux, Initial residual = 0.0047223043, Final residual = 2.1859185e-07, No Iterations 2
DILUPBiCG:  Solving for Uy, Initial residual = 0.0048172404, Final residual = 2.3390652e-07, No Iterations 2
DILUPBiCG:  Solving for Uz, Initial residual = 0.0041326105, Final residual = 1.7710191e-09, No Iterations 3
DILUPBiCG:  Solving for T, Initial residual = 0.00060871525, Final residual = 9.7934538e-10, No Iterations 4
GAMG:  Solving for p_rgh, Initial residual = 0.067143822, Final residual = 8.550094e-09, No Iterations 13
time step continuity errors : sum local = 3.4964727e-13, global = 9.497921e-14, cumulative = 4.9549033e-11
ExecutionTime = 284.84 s  ClockTime = 287 s

Time = 0.000703654485

Interpolating laser coordinate
Interpolating laser coordinate
Vx   = 100
Vy   = 0
Vmag = 100

Courant Number mean: 0.00069140376 max: 0.27616178
deltaT = 9.9667774e-07
PIMPLE: iteration 1
DILUPBiCG:  Solving for Ux, Initial residual = 0.0043974518, Final residual = 4.3794962e-07, No Iterations 2
DILUPBiCG:  Solving for Uy, Initial residual = 0.0044430772, Final residual = 3.5340057e-07, No Iterations 2
DILUPBiCG:  Solving for Uz, Initial residual = 0.0047550922, Final residual = 1.9589747e-09, No Iterations 3
DILUPBiCG:  Solving for T, Initial residual = 0.00061698156, Final residual = 93345.397, No Iterations 1001
GAMG:  Solving for p_rgh, Initial residual = 0.99749556, Final residual = 6.74928e-09, No Iterations 21
time step continuity errors : sum local = 5.4155474e-11, global = 6.3819725e-14, cumulative = 4.9612853e-11
ExecutionTime = 291.17 s  ClockTime = 293 s

Time = 0.0007046511628

Interpolating laser coordinate
Interpolating laser coordinate
Vx   = 100
Vy   = 0
Vmag = 100

Courant Number mean: 0.00073293333 max: 2.4614658
deltaT = 4.0444894e-07
PIMPLE: iteration 1
DILUPBiCG:  Solving for Ux, Initial residual = 0.99996533, Final residual = 4.1997867e-08, No Iterations 3
DILUPBiCG:  Solving for Uy, Initial residual = 0.99996593, Final residual = 7.471817e-07, No Iterations 3
DILUPBiCG:  Solving for Uz, Initial residual = 0.84149551, Final residual = 4.3093462e-09, No Iterations 3
DILUPBiCG:  Solving for T, Initial residual = 1, Final residual = 5.326721e-10, No Iterations 11
GAMG:  Solving for p_rgh, Initial residual = 1, Final residual = 8.2242381e-09, No Iterations 41
time step continuity errors : sum local = 3187537.2, global = -3187537.1, cumulative = -3187537.1
ExecutionTime = 291.84 s  ClockTime = 294 s

Time = 0.0007050556117

Interpolating laser coordinate
Interpolating laser coordinate
Vx   = 100
Vy   = 0
Vmag = 100

Courant Number mean: 5.0625884e+08 max: 2.3247589e+12
deltaT = 1.7397457e-19
PIMPLE: iteration 1
DILUPBiCG:  Solving for Ux, Initial residual = 0.002170343, Final residual = 1.3269969e-08, No Iterations 2
DILUPBiCG:  Solving for Uy, Initial residual = 0.004592665, Final residual = 2.2649814e-07, No Iterations 2
DILUPBiCG:  Solving for Uz, Initial residual = 0.011799993, Final residual = 4.0934735e-10, No Iterations 3
DILUPBiCG:  Solving for T, Initial residual = 1, Final residual = 4.219571e-10, No Iterations 4
GAMG:  Solving for p_rgh, Initial residual = 1, Final residual = 6.0997958e-09, No Iterations 28
time step continuity errors : sum local = 2.5050765e+21, global = -2.5050765e+21, cumulative = -2.5050765e+21
ExecutionTime = 292.35 s  ClockTime = 294 s

Time = 0.0007050556117

Interpolating laser coordinate
Interpolating laser coordinate
Vx   = 79.769059
Vy   = 0
Vmag = 79.769059

Courant Number mean: 4.0417955e+23 max: 1.9357799e+27
deltaT = 8.9873118e-47
PIMPLE: iteration 1
DILUPBiCG:  Solving for Ux, Initial residual = 0.00087760067, Final residual = 5.1400685e-08, No Iterations 2
DILUPBiCG:  Solving for Uy, Initial residual = 0.0011133698, Final residual = 9.93714e-09, No Iterations 2
DILUPBiCG:  Solving for Uz, Initial residual = 0.013898725, Final residual = 2.1458128e-08, No Iterations 3
[5] #0  [9] [4] #0  Foam::error::printStack(Foam::Ostream&)[8] #0  Foam::error::printStack(Foam::Ostream&)Foam::error::printStack(Foam::Ostream&)#0  Foam::error::printStack(Foam::Ostream&)--------------------------------------------------------------------------
An MPI process has executed an operation involving a call to the
"fork()" system call to create a child process.  Open MPI is currently
operating in a condition that could result in memory corruption or
other system errors; your MPI job may hang, crash, or produce silent
data corruption.  The use of fork() (or system() or other calls that
create child processes) is strongly discouraged.  
...
...
...
bmercer is offline   Reply With Quote

Old   November 7, 2017, 14:08
Default
  #2
MSF
New Member
 
Join Date: Apr 2014
Location: Germany
Posts: 24
Rep Power: 11
MSF is on a distinguished road
Hi

There is a paper from Voller saying that you should use harmonic interpolation for non constant thermal conductivity. (TREATMENT OF DISCONTINUOUS THERMAL CONDUCTIVITY IN CONTROL-VOLUME SOLUTIONS OF PHASE-CHANGE PROBLEMS
V. R. Voller & C. R. Swaminathan)
But I don't think that this is your main problem. With a non constant thermal conductivity you are introducing a non linearity into your problem and you have to use iterations to resolve it. I don't know how you are resolving the phase change but you can update the thermal conductivity in that iteration loop.

Best Moritz
MSF is offline   Reply With Quote

Old   November 8, 2017, 18:26
Default
  #3
New Member
 
Join Date: Jan 2017
Posts: 7
Rep Power: 8
bmercer is on a distinguished road
Hi Moritz,

That is a great tip and using a harmonic interpolation scheme for the laplacian does seem to solve my convergence issues. Thank you! Also just a note that my solver is already set up to iterate to handle the nonlinearities in the temperature equation, with temperature and conductivity updates in the appropriate places.

Your answer also led me to to a better understanding of how the laplacian term and non-constant conductivity is handled in OpenFoam/FVM versus how I would actually like it to behave. This may be the wrong subforum to discuss this, but I'll elaborate here:

Let's assume my model consists of phases 1 and 2. As previously described, each phase has conductivity as a linear function of temperature but k2(T) = 100*k1(T). For my specific problem, heat flows well (higher conductivity) between a 1-2 and 2-2 interface but not between a 1-1 interface. If you do harmonic interpolation you can see that it actually gives an interpolated value of a conductivity much closer to the low k1 value than the k2 value at a 1-2 interface, which is not the result I would like for my model. So I would actually like to define conductivity between a 1-1 interface as k1(T) and between a 1-2 and 2-2 interface as k2(T), with T taking the interpolated temperature value between cells. How to do this in OpenFoam is another question altogether, but I would like to ask, is there anything about that approach that is numerically flawed and would lead non-convergence of the solver?
bmercer is offline   Reply With Quote

Old   November 13, 2017, 04:56
Default
  #4
MSF
New Member
 
Join Date: Apr 2014
Location: Germany
Posts: 24
Rep Power: 11
MSF is on a distinguished road
Hi

I think it is possible to define the surfaceScalarField for the laplacian directly. So you could interpolate your phase fraction (let us assume its called alpha) on the cell faces. If alpha is zero you prescribe thermal conductivity 1 to that face. If alpha is greater than zero thermal conductivity 2. I not 100% sure if you can do it directly on the surfaceScalarField with pos/neg function but I am sure that you can do it with a forAll loop.

Best Moritz
MSF is offline   Reply With Quote

Old   November 14, 2017, 13:25
Default
  #5
New Member
 
Join Date: Jan 2017
Posts: 7
Rep Power: 8
bmercer is on a distinguished road
Hi Moritz,

Thanks again for your thoughts. The programmer's manual confirms that you can supply a surfaceScalarField as the coefficient argument to fvm::laplacian, so your method for prescribing thermal conductivity should work for me. I will try it out and see it how it behaves.

Also, I just wanted to note that I was able to get my model to work with linear interpolation of the diffusion coefficient when I reduced the time step. So to summarize for any other interested readers, I was able to avoid non-convergence of my solver by
1) Using harmonic interpolation of the diffusion coefficient (specified in fvSchemes). This solved my convergence issue, although I realized the resulting solution didn't behave how I wanted it to for my specific model.
2) Retaining linear interpolation of the coefficient but reducing the timestep for my transient simulation. A larger reduction was required when I used boundary conditions which induced sharper thermal gradients.

Thanks very much for your replies, Moritz, as I learned a lot and also solved my issue.
bmercer is offline   Reply With Quote

Reply

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
[openSmoke] libOpenSMOKE Tobi OpenFOAM Community Contributions 562 January 25, 2023 09:21
Setting the height of the stream in the free channel kevinmccartin CFX 12 October 13, 2022 21:43
Problem in setting Boundary Condition Madhatter92 CFX 12 January 12, 2016 04:39
Indoor Air Temperature distribution - convergence problem Shogan FLUENT 4 May 8, 2014 02:07
The unit of thermal conductivity when adding temperature equation to icoFoam? mxylondon OpenFOAM 2 February 26, 2013 05:26


All times are GMT -4. The time now is 20:54.