|
[Sponsors] |
Modifying solidificationMeltingSource for Mushy zone phasechange |
|
LinkBack | Thread Tools | Search this Thread | Display Modes |
January 30, 2017, 10:07 |
Modifying solidificationMeltingSource for Mushy zone phasechange
|
#1 |
Member
a
Join Date: Oct 2014
Posts: 49
Rep Power: 12 |
Hi Foam experts,
I was trying to modify the solidificationMeltingSource (Isothermal ohase change source) for the mushy ploblems in OpenFoam version 3.0. I was successful till defining the liquid fraction using Tsolidus and Tliquidus. But some how I am unable to get the "Phi field" in order to add the term " fvc::div(phi, alpha1_))" in solidificationMeltingSourceTemplates.C. Code:
\*---------------------------------------------------------------------------*/ #include "fvMatrices.H" #include "fvcDdt.H" #include "fvcDiv.H" // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // template<class RhoFieldType> void Foam::fv::mushysolidificationMeltingSource::apply ( const RhoFieldType& rho, fvMatrix<scalar>& eqn ) { if (debug) { Info<< type() << ": applying source to " << eqn.psi().name() << endl; } const volScalarField Cp(this->Cp()); update(Cp); dimensionedScalar L("L", dimEnergy/dimMass, L_); // contributions added to rhs of solver equation if (eqn.psi().dimensions() == dimTemperature) { // isothermal phase change - only include time derivative eqn -= L/Cp*(fvc::ddt(rho, alpha1_) + fvc::div(phi, alpha1_)); //mushy phase change // eqn -= L/Cp*(fvc::ddt(rho, alpha1_)); } else { // isothermal phase change - only include time derivative eqn -= L*(fvc::ddt(rho, alpha1_) + fvc::div(phi, alpha1_)); //mushy phase change // eqn -= L*(fvc::ddt(rho, alpha1_)); } } // ************************************************************************* // Code:
In file included from sources/derived/mushysolidificationMeltingSource/mushysolidificationMeltingSource.H:269:0, from sources/derived/mushysolidificationMeltingSource/mushysolidificationMeltingSource.C:26: sources/derived/mushysolidificationMeltingSource/mushysolidificationMeltingSourceTemplates.C: In member function ‘void Foam::fv::mushysolidificationMeltingSource::apply(const RhoFieldType&, Foam::fvMatrix<double>&)’: sources/derived/mushysolidificationMeltingSource/mushysolidificationMeltingSourceTemplates.C:56:57: error: ‘phi’ was not declared in this scope eqn -= L/Cp*(fvc::ddt(rho, alpha1_) + fvc::div(phi, alpha1_)); //mushy phase change ^ sources/derived/mushysolidificationMeltingSource/mushysolidificationMeltingSourceTemplates.C:62:54: error: ‘phi’ was not declared in this scope eqn -= L*(fvc::ddt(rho, alpha1_) + fvc::div(phi, alpha1_)); //mushy phase change ^ sources/derived/mushysolidificationMeltingSource/mushysolidificationMeltingSource.C: In member function ‘void Foam::fv::mushysolidificationMeltingSource::update(const volScalarField&)’: sources/derived/mushysolidificationMeltingSource/mushysolidificationMeltingSource.C:170:16: warning: unused variable ‘Cpc’ [-Wunused-variable] scalar Cpc = Cp[cellI]; ^ make: *** [Make/linuxGccDPInt32Opt/sources/derived/mushysolidificationMeltingSource/mushysolidificationMeltingSource.o] Error 1 Thanks in advance. kindly help |
|
January 30, 2017, 12:46 |
|
#2 |
Senior Member
Chris Sideroff
Join Date: Mar 2009
Location: Ottawa, ON, CAN
Posts: 434
Rep Power: 22 |
The flux, phi, is not stored in that class (or parent classes) so you need to get a reference to it through the registry. Look in one of the other fvOptions sources that requires the flux for an example. The mesh is accessible so this should work:
Code:
const surfaceScalarField& phi = mesh_.lookupObject<surfaceScalarField>(phiName_); |
|
January 31, 2017, 05:46 |
|
#3 | |
Member
a
Join Date: Oct 2014
Posts: 49
Rep Power: 12 |
Quote:
|
||
June 26, 2017, 03:19 |
|
#4 |
New Member
diwakar
Join Date: Sep 2016
Posts: 11
Rep Power: 10 |
hi cfd@kgp
were you able to run make the solver for solidificationMelting using fvOptions for mushy zone?? |
|
July 1, 2017, 14:44 |
|
#5 |
New Member
diwakar
Join Date: Sep 2016
Posts: 11
Rep Power: 10 |
Anyways I got it working for mushy zone ( both linear and schiel law case)
|
|
September 6, 2017, 11:15 |
|
#6 | |
New Member
Maurício Guilherme Alves dos Reis
Join Date: Feb 2015
Posts: 10
Rep Power: 11 |
Quote:
I'm using OpenFOAM 4.1, I have made the proposed modifications, but the following error occurs during compilation: error: ‘div’ is not a member of ‘Foam::fvc’ eqn -= L/Cp*(fvc::ddt(rho, alpha1_) + fvc::div(phi, alpha1_)); can anybody help me? |
||
September 14, 2017, 06:09 |
|
#7 |
Member
Bruno Lebon
Join Date: Dec 2012
Posts: 33
Rep Power: 13 |
Did you include fvcDiv.H?
My mushyZoneSourceTemplates.C file looks like this and it compiles fine: Code:
/*---------------------------------------------------------------------------*\ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | \\ / A nd | Copyright (C) 2014-2015 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License This file is part of OpenFOAM. OpenFOAM is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. OpenFOAM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>. \*---------------------------------------------------------------------------*/ #include "fvMatrices.H" #include "fvcDdt.H" #include "fvcDiv.H" // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // template<class RhoFieldType> void Foam::fv::mushyZoneSource::apply ( const RhoFieldType& rho, fvMatrix<scalar>& eqn ) { if (debug) { Info<< type() << ": applying source to " << eqn.psi().name() << endl; } const surfaceScalarField& phi = mesh_.lookupObject<surfaceScalarField>(phiName_); const volScalarField Cp(this->Cp()); update(); dimensionedScalar L("L", dimEnergy/dimMass, L_); // contributions added to rhs of solver equation if (eqn.psi().dimensions() == dimTemperature) { // isothermal phase change - only include time derivative eqn -= L/Cp*(fvc::ddt(rho, alpha1_) + fvc::div(phi, alpha1_)); // eqn -= L/Cp*(fvc::ddt(rho, alpha1_)); } else { // isothermal phase change - only include time derivative eqn -= L*(fvc::ddt(rho, alpha1_) + fvc::div(phi, alpha1_)); // eqn -= L*(fvc::ddt(rho, alpha1_)); } } // ************************************************************************* // |
|
September 14, 2017, 06:59 |
|
#8 | |
New Member
diwakar
Join Date: Sep 2016
Posts: 11
Rep Power: 10 |
Quote:
-Diwakar Janghel |
||
September 14, 2017, 07:03 |
|
#9 |
Member
Bruno Lebon
Join Date: Dec 2012
Posts: 33
Rep Power: 13 |
Thanks, I was replying to MauricioReis's error ... but since you are here, did you ever try to use a tabulated fraction liquid:
Code:
void Foam::fv::mushyZoneSource::update() { if (curTimeIndex_ == mesh_.time().timeIndex()) { return; } if (debug) { Info<< type() << ": " << name_ << " - updating phase indicator" << endl; } // update old time alpha1 field alpha1_.oldTime(); const volScalarField& T = mesh_.lookupObject<volScalarField>(TName_); interpolationTable<scalar> fraction_curve("constant/fL"); forAll(cells_, i) { label celli = cells_[i]; scalar Tc = T[celli]; scalar alpha1New; if (Tc < Tsolidus_) { alpha1New = 0.0; } else if (Tc < Tliquidus_) { alpha1New = fraction_curve(Tc); } else { alpha1New = 1.0; } alpha1_[celli] = max(0, min(alpha1New, 1)); } alpha1_.correctBoundaryConditions(); curTimeIndex_ = mesh_.time().timeIndex(); } |
|
September 14, 2017, 07:13 |
|
#10 | |
New Member
diwakar
Join Date: Sep 2016
Posts: 11
Rep Power: 10 |
Quote:
i did not get your question? are you talking about schiel law? |
||
September 14, 2017, 07:23 |
|
#11 |
Member
Bruno Lebon
Join Date: Dec 2012
Posts: 33
Rep Power: 13 |
Yes, I implemented the fraction solid dependence on temperature using a temperature table ... which was itself calculated assuming Scheil solidification. Did you ever try this before? It converges fine ... to the wrong temperature profile.
|
|
September 14, 2017, 07:25 |
|
#12 |
New Member
diwakar
Join Date: Sep 2016
Posts: 11
Rep Power: 10 |
Can you share the formula which you are using to relate solid fraction with temperature..?
|
|
September 14, 2017, 07:35 |
|
#13 |
Member
Bruno Lebon
Join Date: Dec 2012
Posts: 33
Rep Power: 13 |
Not a formula, but a table:
Code:
( ( 0 0) (873 0) (883.3 0.05) (884.7 0.1) (885.4 0.2) (886.8 0.3) (887.5 0.4) (891.6 0.5) (900.5 0.6) (907.4 0.7) (912.2 0.8) (913.6 0.9) (915.0 1.0) (9999 1.0) ) |
|
September 14, 2017, 07:54 |
|
#14 |
New Member
diwakar
Join Date: Sep 2016
Posts: 11
Rep Power: 10 |
||
September 14, 2017, 09:59 |
|
#15 | |
New Member
Maurício Guilherme Alves dos Reis
Join Date: Feb 2015
Posts: 10
Rep Power: 11 |
Quote:
|
||
November 15, 2017, 10:48 |
|
#16 |
Senior Member
Join Date: Sep 2013
Posts: 353
Rep Power: 21 |
Can someone explain to me why alpha is updated in the following fashion in the original solidificationMeltingSource? I understand that C*(Tc - Tmelt_) is basically the energy above the melt temperature. We devide by L to get a fraction....why is this needed though? Why not simple if T<Tmelt alpha=0 if T>Tmelt alpha=1...
Code:
scalar alpha1New = alpha1_[celli] + relax_*Cpc*(Tc - Tmelt_)/L_; Code:
alpha1 = 0.5*Foam::erf(4.0*(T-Tmelt)/(Tl-Ts))+scalar(0.5); Last edited by Bloerb; November 15, 2017 at 17:31. |
|
April 12, 2018, 11:53 |
|
#17 | |
Member
a
Join Date: Oct 2014
Posts: 49
Rep Power: 12 |
Quote:
Bloerb/Stefan I am curious how did you implemented it using functions, please let us know your progress. meanwhile, I had send you a messege too, but not sure about its delivery so thought to get your response here. Thanks in advance |
||
May 10, 2018, 05:46 |
|
#18 | |
Member
Bruno Lebon
Join Date: Dec 2012
Posts: 33
Rep Power: 13 |
Quote:
When you have a mushy zone, you need to replace Tmelt_ but the temperature corresponding to alpha1_[celli] (either using a formula -- lever rule, Scheil or whatever --) or a lookup table. |
||
May 10, 2018, 12:13 |
|
#19 |
Member
Tarang
Join Date: Feb 2011
Location: Delhi, India
Posts: 47
Rep Power: 15 |
Hi,
I have made some changes in solidificationMeltingSource and I use it with tinkered buoyantPimpleSolver which I named as meltFoam. I basically solve energy equation before NS. Further, I solve Energy equation 'n' (nEnergyCorrector in PIMPLE dict) times so that energy equation converges reasonably before I move to next PIMPLE corrector. I am attaching both things here, feel free to comment or suggest improvements. |
|
January 25, 2021, 05:33 |
Asking for which OF version?
|
#20 | |
Member
Join Date: Nov 2020
Posts: 53
Rep Power: 5 |
Quote:
|
||
Tags |
melting openfoam, solidification, solidification/melting |
Thread Tools | Search this Thread |
Display Modes | |
|
|
Similar Threads | ||||
Thread | Thread Starter | Forum | Replies | Last Post |
[Commercial meshers] Mesh conversion problem (fluent3DMeshToFoam) | Aadhavan | OpenFOAM Meshing & Mesh Conversion | 2 | March 8, 2018 02:47 |
Possible Bug in pimpleFoam (or createPatch) (or fluent3DMeshToFoam) | cfdonline2mohsen | OpenFOAM | 3 | October 21, 2013 10:28 |
[Commercial meshers] fluentMeshToFoam multidomain mesh conversion problem | Attesz | OpenFOAM Meshing & Mesh Conversion | 12 | May 2, 2013 11:52 |
Problem in running ICEM grid in Openfoam | Tarak | OpenFOAM | 6 | September 9, 2011 18:51 |
Problem in IMPORT of ICEM input file in FLUENT | csvirume | FLUENT | 2 | September 9, 2009 02:08 |