- **OpenFOAM Programming & Development**
(*http://www.cfd-online.com/Forums/openfoam-programming-development/*)

- - **add source term in energy equation**
(*http://www.cfd-online.com/Forums/openfoam-programming-development/109026-add-source-term-energy-equation.html*)

add source term in energy equationHi, Foamers,
I added a source term in temperature equation, but the strange thing is temperature becomes unbounded hence unphysical. I have no any idea. The following are relevant codes. The basic idea is : I added a source term in TEqn, and the source term is also dependent on T. When I deal with ddt(alphaTemp), I already stored .oldTime() each time. I cannot see where the problem would be. ----------------- fvScalarMatrix TEqn ( fvm::ddt(cp, T) + fvm::div(phi*fvc::interpolate(cp), T) - fvm::laplacian(lambda/rho, T) == ST ); ------------------- ST = ( - hs*fvc::ddt(alphaTemp) - hs*fvc::div(phi,alphaTemp) ); ------------------- forAll(mesh.cells(), celli) { if (T[celli] < Ts.value()) { alphaTemp[celli] = 0; } else if (T[celli] > Tl.value()) { alphaTemp[celli] = 1; } else { alphaTemp[celli] = (T[celli]-Ts.value())/(Tl.value()-Ts.value()); } }; ---------------------- Any response is welcome and appreciated. :) |

Any idea? It takes me quite a long time to figure out. :(
Many thanks to any comments. :) |

I haven't looked too deeply at your code, and I don't know if there's some fundamental problem... like, say the source term has the wrong sign. I'll assume you know what you are doing, and everything is okay.
What I was going to mention is that numerical models sometimes behave very badly around switching functions, especially if you are taking derivatives of them. So your if/else structure may be causing the problems. You could try making it a smooth function without any conditionals... it will also make your code faster. In fact, you could probably make it a field operation instead of a cell-by-cell loop. Google logistic functions. Not sure if this will help, just an idea. |

Quote:
Thanks for your comments. I tried: alphaTemp = 0.5*Foam::erf(4.0*(T-Tmelt)/(Tl-Ts))+scalar(0.5) which is a field operation, but still didn't figure out my problem. Also, the sign of source term has no problem. |

I believe this is stability problem about solving TEqn with source term.
----- part of fvSolution ------- T { solver PBiCG; preconditioner DILU; tolerance 1e-7; relTol 0.1; } TFinal { solver PBiCG; preconditioner DILU; tolerance 1e-9; relTol 0; } PIMPLE { momentumPredictor no; nOuterCorrectors 3; nCorrectors 3; nNonOrthogonalCorrectors 0; pRefCell 0; pRefValue 0; } ---------------------------------------------- My problem is a 2D cavity melting, initial temperature of the whole domain is 301.45, BC is left side 311.15 fixed, right side 301.45 fixed and up and bottom insulated. ---------------- part of log (temperature) -------------- Time = 0.05 uniform 301.45; PIMPLE: iteration 2 (306.456 303.191 302.055 301.66 301.523 301.475 301.459 301.453 301.451 301.45) PIMPLE: iteration 3 (218.682 189.548 262.537 287.918 296.744 299.814 300.881 301.253 301.384 301.434); Time = 0.109375 (306.456 303.191 302.055 301.66 301.523 301.475 301.459 301.453 301.451 301.45); PIMPLE: iteration 2 (308.363 304.81 302.974 302.112 301.73 301.565 301.497 301.469 301.457 301.452); PIMPLE: iteration 3 (301.785 278.104 243.453 279.972 293.494 298.502 300.358 301.048 301.307 301.415) -------------------------------------------- I also tried to set pimple iterations even as 25, still didn't work... Is this stability problem, if yes, how can I handle it? Many thanks in advance. |

All times are GMT -4. The time now is 00:44. |