4 * Created on: 1 janv. 2015
5 * Author: Michael Ndjinga
7 #include "DriftModel.hxx"
11 DriftModel::DriftModel(pressureEstimate pEstimate, int dim, bool useDellacherieEOS){
15 _dragCoeffs=vector<double>(2,0);
19 if( pEstimate==around1bar300K){//EOS at 1 bar and 373K
20 cout<<"Fluid is water-Gas mixture around saturation point 1 bar and 373 K (100°C)"<<endl;
21 *_runLogFile<<"Fluid is water-Gas mixture around saturation point 1 bar and 373 K (100°C)"<<endl;
22 _Tsat=373;//saturation temperature at 1 bar
23 double esatv=2.505e6;//Gas internal energy at saturation at 1 bar
24 double esatl=4.174e5;//water internal energy at saturation at 1 bar
25 double cv_v=1555;//Gas specific heat capacity at saturation at 1 bar
26 double cv_l=3770;//water specific heat capacity at saturation at 1 bar
27 double gamma_v=1.337;//Gas heat capacity ratio at saturation at 1 bar
28 double rho_sat_l=958;//water density at saturation at 1 bar
29 double sound_speed_l=1543;//water sound speed at saturation at 1 bar
30 _fluides[0] = new StiffenedGas(gamma_v,cv_v,_Tsat,esatv); //ideal gas law for Gas at pressure 1 bar and temperature 100°C, gamma=1.34
31 _fluides[1] = new StiffenedGas(rho_sat_l,1e5,_Tsat,esatl,sound_speed_l,cv_l); //stiffened gas law for water at pressure 1 bar and temperature 100°C
32 _hsatl=4.175e5;//water enthalpy at saturation at 1 bar
33 _hsatv=2.675e6;//Gas enthalpy at saturation at 1 bar
35 _useDellacherieEOS=false;
37 else{//EOS at 155 bars and 618K
38 cout<<"Fluid is water-Gas mixture around saturation point 155 bars and 618 K (345°C)"<<endl;
39 *_runLogFile<<"Fluid is water-Gas mixture around saturation point 155 bars and 618 K (345°C)"<<endl;
40 _useDellacherieEOS=useDellacherieEOS;
43 _Tsat=656;//saturation temperature used in Dellacherie EOS
44 _hsatl=1.633e6;//water enthalpy at saturation at 155 bars
45 _hsatv=3.006e6;//Gas enthalpy at saturation at 155 bars
46 _fluides[0] = new StiffenedGasDellacherie(1.43,0 ,2.030255e6 ,1040.14); //stiffened gas law for Gas from S. Dellacherie
47 _fluides[1] = new StiffenedGasDellacherie(2.35,1e9,-1.167056e6,1816.2); //stiffened gas law for water from S. Dellacherie
51 double esatv=2.444e6;//Gas internal energy at saturation at 155 bar
52 double esatl=1.604e6;//water internal energy at saturation at 155 bar
53 double sound_speed_v=433;//Gas sound speed at saturation at 155 bar
54 double sound_speed_l=621;//water sound speed at saturation at 155 bar
55 double cv_v=3633;//Gas specific heat capacity at saturation at 155 bar
56 double cv_l=3100;//water specific heat capacity at saturation at 155 bar
57 double rho_sat_v=102;//Gas density at saturation at 155 bar
58 double rho_sat_l=594;//water density at saturation at 155 bar
59 _Tsat=618;//saturation temperature at 155 bars
60 _hsatl=1.63e6;//water enthalpy at saturation at 155 bars
61 _hsatv=2.6e6;//Gas enthalpy at saturation at 155 bars
62 _fluides[0] = new StiffenedGas(rho_sat_v,1.55e7,_Tsat,esatv, sound_speed_v,cv_v); //stiffened gas law for Gas at pressure 155 bar and temperature 345°C
63 _fluides[1] = new StiffenedGas(rho_sat_l,1.55e7,_Tsat,esatl, sound_speed_l,cv_l); //stiffened gas law for water at pressure 155 bar
66 _latentHeat=_hsatv-_hsatl;
67 cout<<"Liquid saturation enthalpy "<< _hsatl<<" J/Kg"<<endl;
68 *_runLogFile<<"Liquid saturation enthalpy "<< _hsatl<<" J/Kg"<<endl;
69 cout<<"Vapour saturation enthalpy "<< _hsatv<<" J/Kg"<<endl;
70 *_runLogFile<<"Vapour saturation enthalpy "<< _hsatv<<" J/Kg"<<endl;
71 cout<<"Latent heat "<< _latentHeat<<endl;
72 *_runLogFile<<"Latent heat "<< _latentHeat<<endl;
75 void DriftModel::initialize(){
76 cout<<"Initialising the drift model"<<endl;
77 *_runLogFile<<"Initialising the drift model"<<endl;
79 _Uroe = new double[_nVar];
80 _gravite = vector<double>(_nVar,0);//Not to be confused with _GravityField3d (size _Ndim). _gravite (size _Nvar) is usefull for dealing with source term and implicitation of gravity vector
81 for(int i=0; i<_Ndim; i++)
82 _gravite[i+2]=_GravityField3d[i];
84 _GravityImplicitationMatrix = new PetscScalar[_nVar*_nVar];
87 _Vitesse=Field("Velocity",CELLS,_mesh,3);//Forcement en dimension 3 (3 composantes) pour le posttraitement des lignes de courant
91 _VoidFraction=Field("Void fraction",CELLS,_mesh,1);
92 _Enthalpy=Field("Enthalpy",CELLS,_mesh,1);
93 _Concentration=Field("Concentration",CELLS,_mesh,1);
94 _Pressure=Field("Pressure",CELLS,_mesh,1);
95 _mixtureDensity=Field("Mixt density",CELLS,_mesh,1);
96 _Temperature=Field("Temperature",CELLS,_mesh,1);
97 _DensiteLiquide=Field("Liquid density",CELLS,_mesh,1);
98 _DensiteVapeur=Field("Steam density",CELLS,_mesh,1);
99 _EnthalpieLiquide=Field("Liquid enthalpy",CELLS,_mesh,1);
100 _EnthalpieVapeur=Field("Steam enthalpy",CELLS,_mesh,1);
101 _VitesseX=Field("Velocity x",CELLS,_mesh,1);
104 _VitesseY=Field("Velocity y",CELLS,_mesh,1);
106 _VitesseZ=Field("Velocity z",CELLS,_mesh,1);
110 if(_entropicCorrection)
111 _entropicShift=vector<double>(3);//at most 3 distinct eigenvalues
113 ProblemFluid::initialize();
116 bool DriftModel::iterateTimeStep(bool &converged)
118 if(_timeScheme == Explicit || !_usePrimitiveVarsInNewton)
119 return ProblemFluid::iterateTimeStep(converged);
124 if(_NEWTON_its>0){//Pas besoin de computeTimeStep à la première iteration de Newton
126 computeTimeStep(stop);
128 if(stop){//Le compute time step ne s'est pas bien passé
129 cout<<"ComputeTimeStep failed"<<endl;
134 computeNewtonVariation();
136 //converged=convergence des iterations
137 if(_timeScheme == Explicit)
139 else{//Implicit scheme
141 KSPGetIterationNumber(_ksp, &_PetscIts);
142 if( _MaxIterLinearSolver < _PetscIts)//save the maximum number of iterations needed during the newton scheme
143 _MaxIterLinearSolver = _PetscIts;
144 if(_PetscIts>=_maxPetscIts)//solving the linear system failed
146 cout<<"Systeme lineaire : pas de convergence de Petsc. Itérations maximales "<<_maxPetscIts<<" atteintes"<<endl;
147 //*_runLogFileogFile<<"Systeme lineaire : pas de convergence de Petsc. Itérations maximales "<<_maxPetscIts<<" atteintes"<<endl;
151 else{//solving the linear system succeeded
152 //Calcul de la variation relative Uk+1-Uk
156 for(int j=1; j<=_Nmailles; j++)
158 for(int k=0; k<_nVar; k++)
161 VecGetValues(_newtonVariation, 1, &I, &dx);
162 VecGetValues(_primitiveVars, 1, &I, &x);
163 if (fabs(x)*fabs(x)< _precision)
165 if(_erreur_rel < fabs(dx))
166 _erreur_rel = fabs(dx);
168 else if(_erreur_rel < fabs(dx/x))
169 _erreur_rel = fabs(dx/x);
173 converged = _erreur_rel <= _precision_Newton;
176 double relaxation=1;//Vk+1=Vk+relaxation*deltaV
178 VecAXPY(_primitiveVars, relaxation, _newtonVariation);
180 //mise a jour du champ primitif
181 updateConservatives();
183 if(_nbPhases==2 && fabs(_err_press_max) > _precision)//la pression n'a pu être calculée en diphasique à partir des variables conservatives
185 cout<<"Warning consToPrim: nbiter max atteint, erreur relative pression= "<<_err_press_max<<" precision= " <<_precision<<endl;
186 *_runLogFile<<"Warning consToPrim: nbiter max atteint, erreur relative pression= "<<_err_press_max<<" precision= " <<_precision<<endl;
192 cout<<"Vecteur Vkp1-Vk "<<endl;
193 VecView(_newtonVariation, PETSC_VIEWER_STDOUT_SELF);
194 cout << "Nouvel etat courant Vk de l'iteration Newton: " << endl;
195 VecView(_primitiveVars, PETSC_VIEWER_STDOUT_SELF);
198 if(_nbPhases==2 && _nbTimeStep%_freqSave ==0){
199 if(_minm1<-_precision || _minm2<-_precision)
201 cout<<"!!!!!!!!! WARNING masse partielle negative sur " << _nbMaillesNeg << " faces, min m1= "<< _minm1 << " , minm2= "<< _minm2<< " precision "<<_precision<<endl;
202 *_runLogFile<<"!!!!!!!!! WARNING masse partielle negative sur " << _nbMaillesNeg << " faces, min m1= "<< _minm1 << " , minm2= "<< _minm2<< " precision "<<_precision<<endl;
206 cout << "!!!!!!!!!!!!!!!!!!!!!!!! Complex eigenvalues on " << _nbVpCplx << " cells, max imag= " << _part_imag_max << endl;
207 *_runLogFile << "!!!!!!!!!!!!!!!!!!!!!!!! Complex eigenvalues on " << _nbVpCplx << " cells, max imag= " << _part_imag_max << endl;
219 void DriftModel::computeNewtonVariation()
221 if(_timeScheme==Explicit || !_usePrimitiveVarsInNewton)
222 ProblemFluid::computeNewtonVariation();
227 cout<<"Vecteur courant Vk "<<endl;
228 VecView(_primitiveVars,PETSC_VIEWER_STDOUT_SELF);
231 if(_timeScheme == Explicit)
233 VecCopy(_b,_newtonVariation);
234 VecScale(_newtonVariation, _dt);
235 if(_verbose && _nbTimeStep%_freqSave ==0)
237 cout<<"Vecteur _newtonVariation =_b*dt"<<endl;
238 VecView(_newtonVariation,PETSC_VIEWER_STDOUT_SELF);
246 cout << "Matrice du système linéaire avant contribution delta t" << endl;
247 MatView(_A,PETSC_VIEWER_STDOUT_SELF);
249 cout << "Second membre du système linéaire avant contribution delta t" << endl;
250 VecView(_b, PETSC_VIEWER_STDOUT_SELF);
253 MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
255 VecAXPY(_b, 1/_dt, _old);
256 VecAXPY(_b, -1/_dt, _conservativeVars);
258 for(int imaille = 0; imaille<_Nmailles; imaille++){
259 _idm[0] = _nVar*imaille;
260 for(int k=1; k<_nVar; k++)
261 _idm[k] = _idm[k-1] + 1;
262 VecGetValues(_primitiveVars, _nVar, _idm, _Vi);
263 primToConsJacobianMatrix(_Vi);
264 for(int k=0; k<_nVar*_nVar; k++)
265 _primToConsJacoMat[k]*=1/_dt;
266 MatSetValuesBlocked(_A, 1, &imaille, 1, &imaille, _primToConsJacoMat, ADD_VALUES);
268 MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY);
270 #if PETSC_VERSION_GREATER_3_5
271 KSPSetOperators(_ksp, _A, _A);
273 KSPSetOperators(_ksp, _A, _A,SAME_NONZERO_PATTERN);
278 cout << "Matrice du système linéaire après contribution delta t" << endl;
279 MatView(_A,PETSC_VIEWER_STDOUT_SELF);
281 cout << "Second membre du système linéaire après contribution delta t" << endl;
282 VecView(_b, PETSC_VIEWER_STDOUT_SELF);
287 KSPSetComputeEigenvalues(_ksp,PETSC_TRUE);
290 KSPSolve(_ksp, _b, _newtonVariation);
294 computeScaling(_maxvp);
296 VecAssemblyBegin(_vecScaling);
297 VecAssemblyBegin(_invVecScaling);
298 for(int imaille = 0; imaille<_Nmailles; imaille++)
301 VecSetValuesBlocked(_vecScaling,1 , &indice, _blockDiag, INSERT_VALUES);
302 VecSetValuesBlocked(_invVecScaling,1,&indice,_invBlockDiag, INSERT_VALUES);
304 VecAssemblyEnd(_vecScaling);
305 VecAssemblyEnd(_invVecScaling);
309 cout << "Matrice avant le preconditionneur des vecteurs propres " << endl;
310 MatView(_A,PETSC_VIEWER_STDOUT_SELF);
312 cout << "Second membre avant le preconditionneur des vecteurs propres " << endl;
313 VecView(_b, PETSC_VIEWER_STDOUT_SELF);
316 MatDiagonalScale(_A,_vecScaling,_invVecScaling);
319 cout << "Matrice apres le preconditionneur des vecteurs propres " << endl;
320 MatView(_A,PETSC_VIEWER_STDOUT_SELF);
323 VecCopy(_b,_bScaling);
324 VecPointwiseMult(_b,_vecScaling,_bScaling);
327 cout << "Produit du second membre par le preconditionneur bloc diagonal a gauche" << endl;
328 VecView(_b, PETSC_VIEWER_STDOUT_SELF);
332 KSPSolve(_ksp,_b, _bScaling);
333 VecPointwiseMult(_newtonVariation,_invVecScaling,_bScaling);
337 cout << "solution du systeme lineaire local:" << endl;
338 VecView(_newtonVariation, PETSC_VIEWER_STDOUT_SELF);
345 void DriftModel::convectionState( const long &i, const long &j, const bool &IsBord){
346 // sortie: etat de Roe rho, cm, v, H,T
349 for(int k=1; k<_nVar; k++)
350 _idm[k] = _idm[k-1] + 1;
351 VecGetValues(_conservativeVars, _nVar, _idm, _Ui);
354 for(int k=1; k<_nVar; k++)
355 _idm[k] = _idm[k-1] + 1;
357 VecGetValues(_Uext, _nVar, _idm, _Uj);
359 VecGetValues(_conservativeVars, _nVar, _idm, _Uj);
360 if(_verbose && _nbTimeStep%_freqSave ==0)
362 cout<<"DriftModel::convectionState Left state cell " << i<< ": "<<endl;
363 for(int k =0; k<_nVar; k++)
365 cout<<"DriftModel::convectionState Right state cell " << j<< ": "<<endl;
366 for(int k =0; k<_nVar; k++)
369 if(_Ui[0]<0||_Uj[0]<0)
371 cout<<"!!!!!!!!!!!!!!!!!!!!!!!!densite de melange negative, arret de calcul!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl;
372 *_runLogFile<<"!!!!!!!!!!!!!!!!!!!!!!!!densite de melange negative, arret de calcul!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl;
373 _runLogFile->close();
374 throw CdmathException("densite negative, arret de calcul");
376 PetscScalar ri, rj, xi, xj, pi, pj;
378 ri = sqrt(_Ui[0]);//racine carre de phi_i rho_i
381 _Uroe[0] = ri*rj/sqrt(_porosityi*_porosityj); //moyenne geometrique des densites
382 if(_verbose && _nbTimeStep%_freqSave ==0)
383 cout << "Densité moyenne Roe gauche " << i << ": " << ri*ri << ", droite " << j << ": " << rj*rj << "->" << _Uroe[0] << endl;
384 xi = _Ui[1]/_Ui[0];//cm gauche
385 xj = _Uj[1]/_Uj[0];//cm droite
387 _Uroe[1] = (xi*ri + xj*rj)/(ri + rj);//moyenne de Roe des concentrations
388 if(_verbose && _nbTimeStep%_freqSave ==0)
389 cout << "Concentration de Roe gauche " << i << ": " << xi << ", droite " << j << ": " << xj << "->" << _Uroe[1] << endl;
390 for(int k=0;k<_Ndim;k++){
391 xi = _Ui[k+2];//phi rho u gauche
392 xj = _Uj[k+2];//phi rho u droite
393 _Uroe[2+k] = (xi/ri + xj/rj)/(ri + rj);
394 //"moyenne" des vitesses
395 if(_verbose && _nbTimeStep%_freqSave ==0)
396 cout << "Vitesse de Roe composante "<< k<<" gauche " << i << ": " << xi/(ri*ri) << ", droite " << j << ": " << xj/(rj*rj) << "->" << _Uroe[k+2] << endl;
398 // H = (rho E + p)/rho
399 xi = _Ui[_nVar-1];//phi rho E
401 Ii = i*_nVar+1; // correct Kieu
402 VecGetValues(_primitiveVars, 1, &Ii, &pi);// _primitiveVars pour p
405 consToPrim(_Uj,_Vj,_porosityj);
410 Ii = j*_nVar+1; // correct Kieu
411 VecGetValues(_primitiveVars, 1, &Ii, &pj);
413 xi = (xi + pi)/_Ui[0];
414 xj = (xj + pj)/_Uj[0];
415 _Uroe[_nVar-1] = (ri*xi + rj*xj)/(ri + rj);
416 if(_verbose && _nbTimeStep%_freqSave ==0)
417 cout << "Enthalpie totale de Roe H gauche " << i << ": " << xi << ", droite " << j << ": " << xj << "->" << _Uroe[_nVar-1] << endl;
418 // Moyenne de Roe de Tg et Td
419 Ii = i*_nVar+_nVar-1;
420 VecGetValues(_primitiveVars, 1, &Ii, &xi);// _primitiveVars pour T
423 //consToPrim(_Uj,_Vj,_porosityj);//Fonction déjà appelée
428 Ii = j*_nVar+_nVar-1;
429 VecGetValues(_primitiveVars, 1, &Ii, &xj);
431 if(_verbose && _nbTimeStep%_freqSave ==0)
433 cout<<"Convection interfacial state"<<endl;
434 for(int k=0;k<_nVar;k++)
435 cout<< _Uroe[k]<<" , "<<endl;
439 void DriftModel::diffusionStateAndMatrices(const long &i,const long &j, const bool &IsBord){
440 //sortie: matrices et etat de diffusion (rho, rho cm, q, T)
443 for(int k=1; k<_nVar; k++)
444 _idm[k] = _idm[k-1] + 1;
446 VecGetValues(_conservativeVars, _nVar, _idm, _Ui);
448 for(int k=1; k<_nVar; k++)
449 _idm[k] = _idm[k-1] + 1;
452 VecGetValues(_Uextdiff, _nVar, _idm, _Uj);
454 VecGetValues(_conservativeVars, _nVar, _idm, _Uj);
456 if(_verbose && _nbTimeStep%_freqSave ==0)
458 cout << "DriftModel::diffusionStateAndMatrices cellule gauche" << i << endl;
460 for(int q=0; q<_nVar; q++)
461 cout << _Ui[q] << "\t";
463 cout << "DriftModel::diffusionStateAndMatrices cellule droite" << j << endl;
465 for(int q=0; q<_nVar; q++)
466 cout << _Uj[q] << "\t";
470 for(int k=0; k<_nVar; k++)
471 _Udiff[k] = (_Ui[k]/_porosityi+_Uj[k]/_porosityj)/2;
472 if(_verbose && _nbTimeStep%_freqSave ==0)
474 cout << "DriftModel::diffusionStateAndMatrices conservative diffusion state" << endl;
476 for(int q=0; q<_nVar; q++)
477 cout << _Udiff[q] << "\t";
479 cout << "porosite gauche= "<<_porosityi<< ", porosite droite= "<<_porosityj<<endl;
482 consToPrim(_Udiff,_phi,1);
483 _Udiff[_nVar-1]=_phi[_nVar-1];
485 double Tm=_phi[_nVar-1];
486 double RhomEm=_Udiff[_nVar-1];
489 if(_timeScheme==Implicit)
492 for (int l = 0; l<_Ndim;l++)
493 q_2+=_Udiff[l+2]*_Udiff[l+2];
494 double pression=_phi[1];
495 double m_v=_Udiff[1];
496 double m_l=_Udiff[0]-_Udiff[1];
497 double rho_v=_fluides[0]->getDensity(pression,Tm);
498 double rho_l=_fluides[1]->getDensity(pression,Tm);
499 double alpha_v=m_v/rho_v,alpha_l=1-alpha_v;
501 for(int l=0; l<_nVar*_nVar;l++)
503 double mu = alpha_v*_fluides[0]->getViscosity(Tm)+alpha_l*_fluides[1]->getViscosity(Tm);
504 for(int l=2;l<(_nVar-1);l++)
506 _Diffusion[l*_nVar] = mu*_Udiff[l]/(_Udiff[0]*_Udiff[0]);
507 _Diffusion[l*_nVar+l] = -mu/_Udiff[0];
509 double lambda = alpha_v*_fluides[0]->getConductivity(Tm)+alpha_l*_fluides[1]->getConductivity(Tm);
510 double C_v= alpha_v*_fluides[0]->constante("cv");
511 double C_l= alpha_l*_fluides[1]->constante("cv");
512 double ev0=_fluides[0]->getInternalEnergy(0,rho_v);//Corriger
513 double el0=_fluides[1]->getInternalEnergy(0,rho_l);//Corriger
514 double Rhomem=RhomEm-0.5*q_2/(_Udiff[0]*_Udiff[0]);
515 int q = (_nVar-1)*_nVar;
516 //Formules a verifier
517 _Diffusion[q]=lambda*((Rhomem-m_v*ev0-m_l*el0)*C_l/((m_v*C_v+m_l*C_l)*(m_v*C_v+m_l*C_l))+el0/(m_v*C_v+m_l*C_l)-q_2/(2*_Udiff[0]*_Udiff[0]*(m_v*C_v+m_l*C_l)));
518 _Diffusion[q+1]=lambda*((Rhomem-m_v*ev0-m_l*el0)*(C_v-C_l)/((m_v*C_v+m_l*C_l)*(m_v*C_v+m_l*C_l))+(ev0+el0)/(m_v*C_v+m_l*C_l));
519 for(int k=2;k<(_nVar-1);k++)
521 _Diffusion[q+k]= lambda*_Udiff[k]/(_Udiff[0]*(m_v*C_v+m_l*C_l));
523 _Diffusion[_nVar*_nVar-1]=-lambda/(m_v*C_v+m_l*C_l);
525 if(_verbose && _nbTimeStep%_freqSave ==0)
527 cout << "Matrice de diffusion D, pour le couple (" << i << "," << j<< "):" << endl;
528 displayMatrix( _Diffusion,_nVar," Matrice de diffusion ");
533 void DriftModel::setBoundaryState(string nameOfGroup, const int &j,double *normale){
535 double v2=0, q_n=0;//q_n=quantité de mouvement normale à la face interne
537 for(k=1; k<_nVar; k++)
538 _idm[k] = _idm[k-1] + 1;
539 VecGetValues(_conservativeVars, _nVar, _idm, _externalStates);//On initialise l'état fantôme avec l'état interne
540 for(k=0; k<_Ndim; k++)
541 q_n+=_externalStates[(k+2)]*normale[k];
543 double porosityj=_porosityField(j);
545 if(_verbose && _nbTimeStep%_freqSave ==0)
547 cout << "setBoundaryState for group "<< nameOfGroup<< ", inner cell j= "<<j << " face unit normal vector "<<endl;
548 for(k=0; k<_Ndim; k++){
549 cout<<normale[k]<<", ";
554 if (_limitField[nameOfGroup].bcType==Wall || _limitField[nameOfGroup].bcType==InnerWall){
555 //Pour la convection, inversion du sens de la vitesse normale
556 for(k=0; k<_Ndim; k++)
557 _externalStates[(k+2)]-= 2*q_n*normale[k];
560 for(k=1; k<_nVar; k++)
561 _idm[k] = _idm[k-1] + 1;
563 VecAssemblyBegin(_Uext);
564 VecSetValues(_Uext, _nVar, _idm, _externalStates, INSERT_VALUES);
565 VecAssemblyEnd(_Uext);
567 //Pour la diffusion, paroi à vitesse et temperature imposees
569 for(k=1; k<_nVar; k++)
570 _idm[k] = _idm[k-1] + 1;
571 VecGetValues(_primitiveVars, _nVar, _idm, _Vj);
572 double concentration=_Vj[0];
573 double pression=_Vj[1];
574 double T=_limitField[nameOfGroup].T;
575 double rho_v=_fluides[0]->getDensity(pression,T);
576 double rho_l=_fluides[1]->getDensity(pression,T);
577 if(fabs(concentration*rho_l+(1-concentration)*rho_v)<_precision)
579 cout<<"rhov= "<<rho_v<<", rhol= "<<rho_l<<endl;
580 cout<<"concentration*rho_l+(1-concentration)*rho_v= "<<concentration*rho_l+(1-concentration)*rho_v<<endl;
581 *_runLogFile<<"concentration*rho_l+(1-concentration)*rho_v= "<<concentration*rho_l+(1-concentration)*rho_v<<endl;
582 _runLogFile->close();
583 throw CdmathException("DriftModel::setBoundaryState: Inlet, impossible to compute mixture density, division by zero");
586 _externalStates[0]=porosityj*rho_v*rho_l/(concentration*rho_l+(1-concentration)*rho_v);
587 _externalStates[1]=concentration*_externalStates[0];
589 _externalStates[2]=_externalStates[0]*_limitField[nameOfGroup].v_x[0];
590 v2 +=_limitField[nameOfGroup].v_x[0]*_limitField[nameOfGroup].v_x[0];
593 v2 +=_limitField[nameOfGroup].v_y[0]*_limitField[nameOfGroup].v_y[0];
594 _externalStates[3]=_externalStates[0]*_limitField[nameOfGroup].v_y[0];
597 _externalStates[4]=_externalStates[0]*_limitField[nameOfGroup].v_z[0];
598 v2 +=_limitField[nameOfGroup].v_z[0]*_limitField[nameOfGroup].v_z[0];
601 _externalStates[_nVar-1] = _externalStates[1]*_fluides[0]->getInternalEnergy(_limitField[nameOfGroup].T,rho_v)
602 +(_externalStates[0]-_externalStates[1])*_fluides[1]->getInternalEnergy(_limitField[nameOfGroup].T,rho_l) + _externalStates[0]*v2/2;
604 for(k=1; k<_nVar; k++)
605 _idm[k] = _idm[k-1] + 1;
606 VecAssemblyBegin(_Uextdiff);
607 VecSetValues(_Uextdiff, _nVar, _idm, _externalStates, INSERT_VALUES);
608 VecAssemblyEnd(_Uextdiff);
610 else if (_limitField[nameOfGroup].bcType==Neumann){
612 for(k=1; k<_nVar; k++)
613 _idm[k] = _idm[k-1] + 1;
615 VecAssemblyBegin(_Uext);
616 VecSetValues(_Uext, _nVar, _idm, _externalStates, INSERT_VALUES);
617 VecAssemblyEnd(_Uext);
619 VecAssemblyBegin(_Uextdiff);
620 VecSetValues(_Uextdiff, _nVar, _idm, _externalStates, INSERT_VALUES);
621 VecAssemblyEnd(_Uextdiff);
623 else if (_limitField[nameOfGroup].bcType==Inlet){
626 VecGetValues(_primitiveVars, _nVar, _idm, _Vj);
627 double concentration=_limitField[nameOfGroup].conc;
628 double pression=_Vj[1];
629 double T=_limitField[nameOfGroup].T;
630 double rho_v=_fluides[0]->getDensity(pression,T);
631 double rho_l=_fluides[1]->getDensity(pression,T);
632 if(fabs(concentration*rho_l+(1-concentration)*rho_v)<_precision)
634 cout<<"rhov= "<<rho_v<<", rhol= "<<rho_l<<endl;
635 cout<<"concentration*rho_l+(1-concentration)*rho_v= "<<concentration*rho_l+(1-concentration)*rho_v<<endl;
636 *_runLogFile<<"concentration*rho_l+(1-concentration)*rho_v= "<<concentration*rho_l+(1-concentration)*rho_v<<endl;
637 _runLogFile->close();
638 throw CdmathException("DriftModel::setBoundaryState: Inlet, impossible to compute mixture density, division by zero");
641 _externalStates[0]=porosityj*rho_v*rho_l/(concentration*rho_l+(1-concentration)*rho_v);
642 _externalStates[1]=concentration*_externalStates[0];
643 _externalStates[2]=_externalStates[0]*(_limitField[nameOfGroup].v_x[0]);
644 v2 +=(_limitField[nameOfGroup].v_x[0])*(_limitField[nameOfGroup].v_x[0]);
647 v2 +=_limitField[nameOfGroup].v_y[0]*_limitField[nameOfGroup].v_y[0];
648 _externalStates[3]=_externalStates[0]*_limitField[nameOfGroup].v_y[0];
651 _externalStates[4]=_externalStates[0]*_limitField[nameOfGroup].v_z[0];
652 v2 +=_limitField[nameOfGroup].v_z[0]*_limitField[nameOfGroup].v_z[0];
655 _externalStates[_nVar-1] = _externalStates[1]*_fluides[0]->getInternalEnergy(T,rho_v)+(_externalStates[0]-_externalStates[1])*_fluides[1]->getInternalEnergy(T,rho_l) + _externalStates[0]*v2/2;
657 else if(_nbTimeStep%_freqSave ==0)
658 cout<< "Warning : fluid going out through inlet boundary "<<nameOfGroup<<". Applying Neumann boundary condition"<<endl;
661 for(k=1; k<_nVar; k++)
662 _idm[k] = _idm[k-1] + 1;
664 VecAssemblyBegin(_Uext);
665 VecAssemblyBegin(_Uextdiff);
666 VecSetValues(_Uext, _nVar, _idm, _externalStates, INSERT_VALUES);
667 VecSetValues(_Uextdiff, _nVar, _idm, _externalStates, INSERT_VALUES);
668 VecAssemblyEnd(_Uext);
669 VecAssemblyEnd(_Uextdiff);
671 else if (_limitField[nameOfGroup].bcType==InletPressure){
673 //Computation of the hydrostatic contribution : scalar product between gravity vector and position vector
674 Cell Cj=_mesh.getCell(j);
675 double hydroPress=Cj.x()*_GravityField3d[0];
677 hydroPress+=Cj.y()*_GravityField3d[1];
679 hydroPress+=Cj.z()*_GravityField3d[2];
681 hydroPress*=_externalStates[0]/porosityj;//multiplication by rho
683 //Building the external state
684 VecGetValues(_primitiveVars, _nVar, _idm, _Vj);
685 double concentration, Tm;
687 concentration=_limitField[nameOfGroup].conc;
688 Tm=_limitField[nameOfGroup].T;
691 if(_nbTimeStep%_freqSave ==0)
692 cout<< "Warning : fluid going out through inletPressure boundary "<<nameOfGroup<<". Applying Neumann boundary condition for concentration, velocity and temperature"<<endl;
693 concentration=_Vj[0];
697 double pression=_limitField[nameOfGroup].p + hydroPress;
698 double rho_v=_fluides[0]->getDensity(pression, Tm);
699 double rho_l=_fluides[1]->getDensity(pression, Tm);
700 if(fabs(concentration*rho_l+(1-concentration)*rho_v)<_precision)
702 cout<<"rhov= "<<rho_v<<", rhol= "<<rho_l<<endl;
703 cout<<"concentration*rho_l+(1-concentration)*rho_v= "<<concentration*rho_l+(1-concentration)*rho_v<<endl;
704 *_runLogFile<<"concentration*rho_l+(1-concentration)*rho_v= "<<concentration*rho_l+(1-concentration)*rho_v<<endl;
705 _runLogFile->close();
706 throw CdmathException("DriftModel::jacobian: Inlet, impossible to compute mixture density, division by zero");
709 _externalStates[0]=porosityj*rho_v*rho_l/(concentration*rho_l+(1-concentration)*rho_v);
710 _externalStates[1]=concentration*_externalStates[0];
711 double mv=_externalStates[1], ml=_externalStates[0]-_externalStates[1];
712 for(k=0; k<_Ndim; k++)
714 v2+=_Vj[k+2]*_Vj[k+2];
715 _externalStates[k+2]=_externalStates[0]*_Vj[(k+2)] ;
717 _externalStates[_nVar-1] = mv*_fluides[0]->getInternalEnergy(Tm,rho_v)+ml*_fluides[1]->getInternalEnergy(Tm,rho_l) +_externalStates[0]* v2/2;
720 for(k=1; k<_nVar; k++)
721 _idm[k] = _idm[k-1] + 1;
722 VecAssemblyBegin(_Uext);
723 VecAssemblyBegin(_Uextdiff);
724 VecSetValues(_Uext, _nVar, _idm, _externalStates, INSERT_VALUES);
725 VecSetValues(_Uextdiff, _nVar, _idm, _externalStates, INSERT_VALUES);
726 VecAssemblyEnd(_Uext);
727 VecAssemblyEnd(_Uextdiff);
729 else if (_limitField[nameOfGroup].bcType==Outlet){
730 if(q_n<=0 && _nbTimeStep%_freqSave ==0)
731 cout<< "Warning : fluid going in through outlet boundary "<<nameOfGroup<<". Applying Neumann boundary condition for concentration, velocity and temperature"<<endl;
733 //Computation of the hydrostatic contribution : scalar product between gravity vector and position vector
734 Cell Cj=_mesh.getCell(j);
735 double hydroPress=(Cj.x()-_gravityReferencePoint[0])*_GravityField3d[0];
737 hydroPress+=(Cj.y()-_gravityReferencePoint[1])*_GravityField3d[1];
739 hydroPress+=(Cj.z()-_gravityReferencePoint[2])*_GravityField3d[2];
741 hydroPress*=_externalStates[0]/porosityj;//multiplication by rho
743 //Building the external state
744 VecGetValues(_primitiveVars, _nVar, _idm, _Vj);
746 double concentration=_Vj[0];
747 double pression=_limitField[nameOfGroup].p+hydroPress;
748 double Tm=_Vj[_nVar-1];
749 double rho_v=_fluides[0]->getDensity(pression, Tm);
750 double rho_l=_fluides[1]->getDensity(pression, Tm);
751 if(fabs(concentration*rho_l+(1-concentration)*rho_v)<_precision)
753 cout<<"rhov= "<<rho_v<<", rhol= "<<rho_l<<endl;
754 cout<<"concentration*rho_l+(1-concentration)*rho_v= "<<concentration*rho_l+(1-concentration)*rho_v<<endl;
755 *_runLogFile<<"concentration*rho_l+(1-concentration)*rho_v= "<<concentration*rho_l+(1-concentration)*rho_v<<endl;
756 _runLogFile->close();
757 throw CdmathException("DriftModel::jacobian: Inlet, impossible to compute mixture density, division by zero");
760 _externalStates[0]=porosityj*rho_v*rho_l/(concentration*rho_l+(1-concentration)*rho_v);
761 _externalStates[1]=concentration*_externalStates[0];
762 double mv=_externalStates[1], ml=_externalStates[0]-_externalStates[1];
763 for(k=0; k<_Ndim; k++)
765 v2+=_Vj[k+2]*_Vj[k+2];
766 _externalStates[k+2]=_externalStates[0]*_Vj[(k+2)] ;
768 _externalStates[_nVar-1] = mv*_fluides[0]->getInternalEnergy(Tm,rho_v)+ml*_fluides[1]->getInternalEnergy(Tm,rho_l) +_externalStates[0]* v2/2;
771 for(k=1; k<_nVar; k++)
772 _idm[k] = _idm[k-1] + 1;
773 VecAssemblyBegin(_Uext);
774 VecAssemblyBegin(_Uextdiff);
775 VecSetValues(_Uext, _nVar, _idm, _externalStates, INSERT_VALUES);
776 VecSetValues(_Uextdiff, _nVar, _idm, _externalStates, INSERT_VALUES);
777 VecAssemblyEnd(_Uext);
778 VecAssemblyEnd(_Uextdiff);
781 cout<<"!!!!!!!!!!!!!!! Error DriftModel::setBoundaryState !!!!!!!!!!"<<endl;
782 cout<<"!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Boundary condition not set for boundary named "<<nameOfGroup<< ", _limitField[nameOfGroup].bcType= "<<_limitField[nameOfGroup].bcType<<endl;
783 cout<<"Accepted boundary condition are Neumann, Wall, InnerWall, Inlet, and Outlet"<<endl;
784 *_runLogFile<<"Boundary condition not set for boundary named "<<nameOfGroup<<"Accepted boundary condition are Neumann,Wall, InnerWall, Inlet, and Outlet"<<endl;
785 _runLogFile->close();
786 throw CdmathException("Unknown boundary condition");
790 void DriftModel::convectionMatrices()
792 //entree: URoe = rho, cm, u, H
793 //sortie: matrices Roe+ et Roe- +Roe si scheme centre
795 if(_verbose && _nbTimeStep%_freqSave ==0)
796 cout<<"DriftModel::convectionMatrices()"<<endl;
798 double umn=0, u_2=0; //valeur de u.normale et |u|²
799 for(int i=0;i<_Ndim;i++)
801 u_2 += _Uroe[2+i]*_Uroe[2+i];
802 umn += _Uroe[2+i]*_vec_normal[i];//vitesse normale
805 vector<complex<double> > vp_dist(3);
807 if(_spaceScheme==staggered && _nonLinearFormulation==VFFC)//special case where we need no Roe matrix
809 staggeredVFFCMatricesConservativeVariables(umn);//Computation of classical upwinding matrices
810 if(_timeScheme==Implicit && _usePrimitiveVarsInNewton)//For use in implicit matrix
811 staggeredVFFCMatricesPrimitiveVariables(umn);
813 else//In this case we build a Roe matrix
815 double rhom=_Uroe[0];
817 double Hm=_Uroe[_nVar-1];
818 double hm=Hm-0.5*u_2;
819 double m_v=cm*rhom, m_l=rhom-m_v;
821 Vector vitesse(_Ndim);
822 for(int idim=0;idim<_Ndim;idim++)
823 vitesse[idim]=_Uroe[2+idim];
825 if(cm<_precision)//pure liquid
826 Tm=_fluides[1]->getTemperatureFromEnthalpy(hm,rhom);
827 else if(cm>1-_precision)
828 Tm=_fluides[0]->getTemperatureFromEnthalpy(hm,rhom);
829 else//Hypothèse de saturation
832 double pression= getMixturePressure(cm, rhom, Tm);
834 if(_verbose && _nbTimeStep%_freqSave ==0)
835 cout<<"Roe state rhom="<<rhom<<" cm= "<<cm<<" Hm= "<<Hm<<" Tm= "<<Tm<<" pression "<<pression<<endl;
837 getMixturePressureDerivatives( m_v, m_l, pression, Tm);//EOS is involved to express pressure jump and sound speed
838 if(_kappa*hm+_khi+cm*_ksi<0){
839 *_runLogFile<<"Calcul matrice de Roe: vitesse du son complexe"<<endl;
840 _runLogFile->close();
841 throw CdmathException("Calcul matrice de Roe: vitesse du son complexe");
843 double am=sqrt(_kappa*hm+_khi+cm*_ksi);//vitesse du son du melange
844 if(_verbose && _nbTimeStep%_freqSave ==0)
845 cout<<"DriftModel::convectionMatrices, sound speed am= "<<am<<endl;
847 //On remplit la matrice de Roe
849 RoeMatrixConservativeVariables( cm, umn, ecin, Hm,vitesse);
851 //On remplit les valeurs propres
856 _maxvploc=fabs(umn)+am;
860 /* Construction des matrices de decentrement */
861 if(_spaceScheme== centered){
862 if(_entropicCorrection)
864 *_runLogFile<<"DriftModel::convectionMatrices: entropy schemes not yet available for drift model"<<endl;
865 _runLogFile->close();
866 throw CdmathException("DriftModel::convectionMatrices: entropy schemes not yet available for drift model");
869 for(int i=0; i<_nVar*_nVar;i++)
871 if(_timeScheme==Implicit && _usePrimitiveVarsInNewton)
872 for(int i=0; i<_nVar*_nVar;i++)
873 _absAroeImplicit[i] = 0;
875 else if( _spaceScheme ==staggered){
876 //Calcul de décentrement de type décalé pour formulation Roe
877 staggeredRoeUpwindingMatrixConservativeVariables( cm, umn, ecin, Hm, vitesse);
878 //staggeredRoeUpwindingMatrixPrimitiveVariables( cm, umn, ecin, Hm, vitesse);
880 else if(_spaceScheme == upwind || _spaceScheme ==pressureCorrection || _spaceScheme ==lowMach ){
881 if(_entropicCorrection)
882 entropicShift(_vec_normal);
884 _entropicShift=vector<double>(3,0);//at most 3 distinct eigenvalues
886 vector< complex< double > > y (3,0);
888 for( int i=0 ; i<3 ; i++)
889 y[i] = Poly.abs_generalise(vp_dist[i])+1*_entropicShift[i];
890 Poly.abs_par_interp_directe(3,vp_dist, _Aroe, _nVar,_precision, _absAroe,y);
892 if( _spaceScheme ==pressureCorrection)
894 for( int i=0 ; i<_Ndim ; i++)
895 for( int j=0 ; j<_Ndim ; j++)
896 _absAroe[(2+i)*_nVar+2+j]-=(vp_dist[2].real()-vp_dist[0].real())/2*_vec_normal[i]*_vec_normal[j];
898 else if( _spaceScheme ==lowMach){
899 double M=sqrt(u_2)/am;
900 for( int i=0 ; i<_Ndim ; i++)
901 for( int j=0 ; j<_Ndim ; j++)
902 _absAroe[(2+i)*_nVar+2+j]-=(1-M)*(vp_dist[2].real()-vp_dist[0].real())/2*_vec_normal[i]*_vec_normal[j];
907 *_runLogFile<<"DriftModel::convectionMatrices: scheme not treated"<<endl;
908 _runLogFile->close();
909 throw CdmathException("DriftModel::convectionMatrices: scheme not treated");
912 for(int i=0; i<_nVar*_nVar;i++)
914 _AroeMinus[i] = (_Aroe[i]-_absAroe[i])/2;
915 _AroePlus[i] = (_Aroe[i]+_absAroe[i])/2;
917 if(_timeScheme==Implicit)
919 if(_usePrimitiveVarsInNewton)//Implicitation using primitive variables
922 _Vij[1]=pression;//pressure
923 _Vij[_nVar-1]=Tm;//Temperature
924 for(int idim=0;idim<_Ndim; idim++)
925 _Vij[2+idim]=_Uroe[2+idim];
926 primToConsJacobianMatrix(_Vij);
928 Poly.matrixProduct(_AroeMinus, _nVar, _nVar, _primToConsJacoMat, _nVar, _nVar, _AroeMinusImplicit);
929 Poly.matrixProduct(_AroePlus, _nVar, _nVar, _primToConsJacoMat, _nVar, _nVar, _AroePlusImplicit);
932 for(int i=0; i<_nVar*_nVar;i++)
934 _AroeMinusImplicit[i] = _AroeMinus[i];
935 _AroePlusImplicit[i] = _AroePlus[i];
938 if(_verbose && _nbTimeStep%_freqSave ==0)
940 displayMatrix(_Aroe, _nVar,"Matrice de Roe");
941 displayMatrix(_absAroe, _nVar,"Valeur absolue matrice de Roe");
942 displayMatrix(_AroeMinus, _nVar,"Matrice _AroeMinus");
943 displayMatrix(_AroePlus, _nVar,"Matrice _AroePlus");
947 if(_verbose && _nbTimeStep%_freqSave ==0 && _timeScheme==Implicit)
949 displayMatrix(_AroeMinusImplicit, _nVar,"Matrice _AroeMinusImplicit");
950 displayMatrix(_AroePlusImplicit, _nVar,"Matrice _AroePlusImplicit");
953 /*******Calcul de la matrice signe pour VFFC, VFRoe et décentrement des termes source***/
954 if(_entropicCorrection)
956 InvMatriceRoe( vp_dist);
958 Poly.matrixProduct(_absAroe, _nVar, _nVar, _invAroe, _nVar, _nVar, _signAroe);
960 else if (_spaceScheme==upwind || (_spaceScheme ==lowMach) || (_spaceScheme ==pressureCorrection))//upwind sans entropic
961 SigneMatriceRoe( vp_dist);
962 else if(_spaceScheme== centered)//centre sans entropic
963 for(int i=0; i<_nVar*_nVar;i++)
965 else if(_spaceScheme ==staggered )
974 for(int i=0; i<_nVar*_nVar;i++)
976 _signAroe[0] = signu;
977 _signAroe[1+_nVar] = signu;
978 for(int i=2; i<_nVar-1;i++)
979 _signAroe[i*_nVar+i] = -signu;
980 _signAroe[_nVar*(_nVar-1)+_nVar-1] = signu;
984 *_runLogFile<<"DriftModel::convectionMatrices: well balanced option not treated"<<endl;
985 _runLogFile->close();
986 throw CdmathException("DriftModel::convectionMatrices: well balanced option not treated");
989 if(_verbose && _nbTimeStep%_freqSave ==0 && _timeScheme==Implicit)
990 displayMatrix(_signAroe, _nVar,"Signe de la matrice de Roe");
993 void DriftModel::addDiffusionToSecondMember
998 double Tm=_Udiff[_nVar-1];
999 double lambdal=_fluides[1]->getConductivity(Tm);
1000 double lambdav = _fluides[0]->getConductivity(Tm);
1001 double mu_l = _fluides[1]->getViscosity(Tm);
1002 double mu_v = _fluides[0]->getViscosity(Tm);
1004 if(mu_v==0 && mu_l ==0 && lambdav==0 && lambdal==0 && _heatTransfertCoeff==0)
1007 //extraction des valeurs
1008 for(int k=0; k<_nVar; k++)
1009 _idm[k] = _nVar*i + k;
1011 VecGetValues(_primitiveVars, _nVar, _idm, _Vi);
1012 if (_verbose && _nbTimeStep%_freqSave ==0)
1014 cout << "Contribution diffusion: variables primitives maille " << i<<endl;
1015 for(int q=0; q<_nVar; q++)
1016 cout << _Vi[q] << endl;
1021 for(int k=0; k<_nVar; k++)
1022 _idn[k] = _nVar*j + k;
1024 VecGetValues(_primitiveVars, _nVar, _idn, _Vj);
1028 lambdal=max(lambdal,_heatTransfertCoeff);//wall nucleate boing -> larger heat transfer
1029 for(int k=0; k<_nVar; k++)
1032 VecGetValues(_Uextdiff, _nVar, _idn, _Uj);
1033 consToPrim(_Uj,_Vj,1);
1035 if (_verbose && _nbTimeStep%_freqSave ==0)
1037 cout << "Contribution diffusion: variables primitives maille " <<j <<endl;
1038 for(int q=0; q<_nVar; q++)
1039 cout << _Vj[q] << endl;
1042 double pression=(_Vi[1]+_Vj[1])/2;//ameliorer car traitement different pour pression et temperature
1043 double m_v=_Udiff[1];
1044 double rho_v=_fluides[0]->getDensity(pression,Tm);
1045 double rho_l=_fluides[1]->getDensity(pression,Tm);
1046 double alpha_v=m_v/rho_v,alpha_l=1-alpha_v;
1047 double mu = alpha_v*mu_v+alpha_l*mu_l;
1048 double lambda = alpha_v*lambdav+alpha_l*lambdal;
1050 //pas de diffusion sur la masse totale
1052 //on n'a pas encore mis la contribution sur la masse
1054 //contribution visqueuse sur la quantite de mouvement
1055 for(int k=2; k<_nVar-1; k++)
1056 _phi[k] = _inv_dxi*2/(1/_inv_dxi+1/_inv_dxj)*mu*(_porosityj*_Vj[k] - _porosityi*_Vi[k]);
1058 //contribution visqueuse sur l'energie
1059 _phi[_nVar-1] = _inv_dxi*2/(1/_inv_dxi+1/_inv_dxj)*lambda*(_porosityj*_Vj[_nVar-1] - _porosityi*_Vi[_nVar-1]);
1062 VecSetValuesBlocked(_b, 1, _idm, _phi, ADD_VALUES);
1064 if(_verbose && _nbTimeStep%_freqSave ==0)
1066 cout << "Contribution diffusion au 2nd membre pour la maille " << i << ": "<<endl;
1067 for(int q=0; q<_nVar; q++)
1068 cout << _phi[q] << endl;
1074 //On change de signe pour l'autre contribution
1075 for(int k=0; k<_nVar; k++)
1076 _phi[k] *= -_inv_dxj/_inv_dxi;
1079 VecSetValuesBlocked(_b, 1, _idn, _phi, ADD_VALUES);
1081 if(_verbose && _nbTimeStep%_freqSave ==0)
1083 cout << "Contribution diffusion au 2nd membre pour la maille " << j << ": "<<endl;
1084 for(int q=0; q<_nVar; q++)
1085 cout << _phi[q] << endl;
1092 void DriftModel::sourceVector(PetscScalar * Si,PetscScalar * Ui,PetscScalar * Vi, int i)
1094 double phirho=Ui[0],phim1=Ui[1],phim2=phirho-phim1,phirhoE=Ui[_nVar-1], cv=Vi[0], P=Vi[1], T=Vi[_nVar-1];
1096 for(int k=0; k<_Ndim; k++)
1097 norm_u+=Vi[2+k]*Vi[2+k];
1098 norm_u=sqrt(norm_u);
1099 double h=(phirhoE-0.5*phirho*norm_u*norm_u+P*_porosityField(i))/phirho;//e+p/rho
1102 //if(T>_Tsat && cv<1-_precision)
1103 if(_hsatv>h && h>_hsatl && cv<1-_precision)//heated boiling//
1104 Si[1]=_heatPowerField(i)/_latentHeat*_porosityField(i);//phi*gamma
1105 else if (P<_Psat && cv<1-_precision)//flash boiling
1106 Si[1]=-_dHsatl_over_dp*_dp_over_dt(i)/_latentHeat;
1109 for(int k=2; k<_nVar-1; k++)
1110 Si[k] =_gravite[k]*phirho-(phim1*_dragCoeffs[0]+phim2*_dragCoeffs[1])*norm_u*Vi[k];
1112 Si[_nVar-1]=_heatPowerField(i);
1114 for(int k=0; k<_Ndim; k++)
1115 Si[_nVar-1] +=(_GravityField3d[k]*phirho-(phim1*_dragCoeffs[0]+phim2*_dragCoeffs[1])*norm_u*Vi[2+k])*Vi[2+k];
1117 if(_timeScheme==Implicit)
1119 for(int k=0; k<_nVar*_nVar;k++)
1120 _GravityImplicitationMatrix[k] = 0;
1121 if(!_usePrimitiveVarsInNewton)
1122 for(int k=0; k<_nVar;k++)
1123 _GravityImplicitationMatrix[k*_nVar]=-_gravite[k];
1126 getDensityDerivatives( cv, P, T ,norm_u*norm_u);
1127 for(int k=0; k<_nVar;k++)
1129 _GravityImplicitationMatrix[k*_nVar+0] =-_gravite[k]*_drho_sur_dcv;
1130 _GravityImplicitationMatrix[k*_nVar+1] =-_gravite[k]*_drho_sur_dp;
1131 _GravityImplicitationMatrix[k*_nVar+_nVar-1]=-_gravite[k]*_drho_sur_dT;
1136 if(_verbose && _nbTimeStep%_freqSave ==0)
1138 cout<<"DriftModel::sourceVector"<<endl;
1140 for(int k=0;k<_nVar;k++)
1144 for(int k=0;k<_nVar;k++)
1148 for(int k=0;k<_nVar;k++)
1151 if(_timeScheme==Implicit)
1152 displayMatrix(_GravityImplicitationMatrix, _nVar, "Gravity implicitation matrix");
1156 void DriftModel::pressureLossVector(PetscScalar * pressureLoss, double K, PetscScalar * Ui, PetscScalar * Vi, PetscScalar * Uj, PetscScalar * Vj)
1158 double norm_u=0, u_n=0, phirho;
1159 for(int i=0;i<_Ndim;i++)
1160 u_n += _Uroe[2+i]*_vec_normal[i];
1165 for(int i=0;i<_Ndim;i++)
1166 norm_u += Vi[2+i]*Vi[2+i];
1167 norm_u=sqrt(norm_u);
1169 for(int i=0;i<_Ndim;i++)
1170 pressureLoss[2+i]=-K*phirho*norm_u*Vi[2+i];
1173 for(int i=0;i<_Ndim;i++)
1174 norm_u += Vj[2+i]*Vj[2+i];
1175 norm_u=sqrt(norm_u);
1177 for(int i=0;i<_Ndim;i++)
1178 pressureLoss[2+i]=-K*phirho*norm_u*Vj[2+i];
1180 pressureLoss[_nVar-1]=-K*phirho*norm_u*norm_u*norm_u;
1183 if(_verbose && _nbTimeStep%_freqSave ==0)
1185 cout<<"DriftModel::pressureLossVector K= "<<K<<endl;
1186 cout<<"pressure loss vector phirho= "<< phirho<<" norm_u= "<<norm_u<<endl;
1188 for(int k=0;k<_nVar;k++)
1192 for(int k=0;k<_nVar;k++)
1196 for(int k=0;k<_nVar;k++)
1200 for(int k=0;k<_nVar;k++)
1203 cout<<"pressureLoss="<<endl;
1204 for(int k=0;k<_nVar;k++)
1205 cout<<pressureLoss[k]<<", ";
1209 void DriftModel::porosityGradientSourceVector()
1211 double u_ni=0, u_nj=0, rhoi,rhoj, pi=_Vi[1], pj=_Vj[1],pij;
1212 for(int i=0;i<_Ndim;i++) {
1213 u_ni += _Vi[2+i]*_vec_normal[i];
1214 u_nj += _Vj[2+i]*_vec_normal[i];
1216 _porosityGradientSourceVector[0]=0;
1217 _porosityGradientSourceVector[1]=0;
1218 rhoj=_Uj[0]/_porosityj;
1219 rhoi=_Ui[0]/_porosityi;
1220 pij=(pi+pj)/2+rhoi*rhoj/2/(rhoi+rhoj)*(u_ni-u_nj)*(u_ni-u_nj);
1221 for(int i=0;i<_Ndim;i++)
1222 _porosityGradientSourceVector[2+i]=pij*(_porosityi-_porosityj)*2/(1/_inv_dxi+1/_inv_dxj);
1223 _porosityGradientSourceVector[_nVar-1]=0;
1226 void DriftModel::jacobian(const int &j, string nameOfGroup,double * normale)
1228 if(_verbose && _nbTimeStep%_freqSave ==0)
1229 cout<<"Jacobienne condition limite convection bord "<< nameOfGroup<<endl;
1232 for(k=0; k<_nVar*_nVar;k++)
1233 _Jcb[k] = 0;//No implicitation at this stage
1236 for(k=1; k<_nVar; k++)
1237 _idm[k] = _idm[k-1] + 1;
1238 VecGetValues(_conservativeVars, _nVar, _idm, _externalStates);
1239 double q_n=0;//quantité de mouvement normale à la paroi
1240 for(k=0; k<_Ndim; k++)
1241 q_n+=_externalStates[(k+1)]*normale[k];
1243 // loop of boundary types
1244 if (_limitField[nameOfGroup].bcType==Wall || _limitField[nameOfGroup].bcType==InnerWall)
1246 for(k=0; k<_nVar;k++)
1247 _Jcb[k*_nVar + k] = 1;
1248 for(k=2; k<_nVar-1;k++)
1249 for(int l=2; l<_nVar-1;l++)
1250 _Jcb[k*_nVar + l] -= 2*normale[k-2]*normale[l-2];
1252 else if (_limitField[nameOfGroup].bcType==Inlet)
1254 return;//For the moment no implicitation, should debug inlet
1256 //Boundary state quantities
1257 double v[_Ndim], ve[_Ndim], v2, ve2;
1258 VecGetValues(_primitiveVars, _nVar, _idm, _Vj);
1259 double concentration=_limitField[nameOfGroup].conc;
1260 double pression=_Vj[1];
1261 double T=_limitField[nameOfGroup].T;
1262 double rho_v=_fluides[0]->getDensity(pression,T);
1263 double rho_l=_fluides[1]->getDensity(pression,T);
1264 if(fabs(concentration*rho_l+(1-concentration)*rho_v)<_precision)
1266 cout<<"rhov= "<<rho_v<<", rhol= "<<rho_l<<endl;
1267 cout<<"concentration*rho_l+(1-concentration)*rho_v= "<<concentration*rho_l+(1-concentration)*rho_v<<endl;
1268 *_runLogFile<<"concentration*rho_l+(1-concentration)*rho_v= "<<concentration*rho_l+(1-concentration)*rho_v<<endl;
1269 _runLogFile->close();
1270 throw CdmathException("DriftModel::jacobian: Inlet, impossible to compute mixture density, division by zero");
1273 double rho_e=rho_v*rho_l/(concentration*rho_l+(1-concentration)*rho_v);
1274 double e_v=_fluides[0]->getInternalEnergy(T,rho_v);
1275 double e_l=_fluides[1]->getInternalEnergy(T,rho_l);
1276 ve[0] = _limitField[nameOfGroup].v_x[0];
1281 ve[1] = _limitField[nameOfGroup].v_y[0];
1287 ve[2] = _limitField[nameOfGroup].v_z[0];
1292 double E_e = concentration*e_v+(1-concentration)*e_l + ve2/2;
1294 //Pressure differential
1295 double gamma_v =_fluides[0]->constante("gamma");
1296 double gamma_l =_fluides[1]->constante("gamma");
1297 double omega=concentration/((gamma_v-1)*rho_v*rho_v*e_v)+(1-concentration)/((gamma_l-1)*rho_l*rho_l*e_l);
1298 double rhom=_externalStates[0];
1299 double m_v=concentration*rhom, m_l=rhom-m_v;
1300 getMixturePressureDerivatives( m_v, m_l, pression, T);
1302 double CoeffCol1 = rho_e*rho_e*omega*(_khi+_kappa*v2/2);
1303 double CoeffCol2 = rho_e*rho_e*omega*_ksi;
1304 double CoeffCol3et4 = rho_e*rho_e*omega*_kappa;
1309 for(k=0;k<_Ndim;k++)
1310 _Jcb[2+k]=-CoeffCol3et4*v[k];
1311 _Jcb[_nVar-1]=CoeffCol3et4;
1313 for(k=0;k<_nVar;k++)
1314 _Jcb[_nVar+k]=_Jcb[k]*concentration;
1316 for(int l=0;l<_Ndim;l++)
1317 for(k=0;k<_nVar;k++)
1318 _Jcb[(2+l)*_nVar+k]=_Jcb[k]*ve[l];
1320 for(k=0;k<_nVar;k++)
1321 _Jcb[(_nVar-1)*_nVar+k]=_Jcb[k]*E_e;
1324 for(k=0;k<_nVar;k++)
1327 else if (_limitField[nameOfGroup].bcType==InletPressure && q_n<0)
1329 return;//For the moment no implicitation, should debug inletPressure
1330 //Boundary state quantities
1331 double v[_Ndim], v2=0;
1332 VecGetValues(_primitiveVars, _nVar, _idm, _Vj);
1333 double concentration=_limitField[nameOfGroup].conc;
1334 double pression=_limitField[nameOfGroup].p;
1335 double T=_limitField[nameOfGroup].T;
1336 double rho_v=_fluides[0]->getDensity(pression,T);
1337 double rho_l=_fluides[1]->getDensity(pression,T);
1338 if(fabs(concentration*rho_l+(1-concentration)*rho_v)<_precision)
1340 cout<<"rhov= "<<rho_v<<", rhol= "<<rho_l<<endl;
1341 cout<<"concentration*rho_l+(1-concentration)*rho_v= "<<concentration*rho_l+(1-concentration)*rho_v<<endl;
1342 *_runLogFile<<"concentration*rho_l+(1-concentration)*rho_v= "<<concentration*rho_l+(1-concentration)*rho_v<<endl;
1343 _runLogFile->close();
1344 throw CdmathException("DriftModel::jacobian: InletPressure, impossible to compute mixture density, division by zero");
1347 double rho_ext=rho_v*rho_l/(concentration*rho_l+(1-concentration)*rho_v);
1348 double rho_int=_externalStates[0];
1349 double ratio_densite=rho_ext/rho_int;
1350 for(k=0;k<_Ndim;k++){
1355 for(int l=0;l<_Ndim;l++){
1356 _Jcb[(2+l)*_nVar]=-ratio_densite*v[l];
1357 _Jcb[(2+l)*_nVar+2+l]=ratio_densite;
1360 _Jcb[(_nVar-1)*_nVar]=-ratio_densite*v2;
1361 for(int l=0;l<_Ndim;l++)
1362 _Jcb[(_nVar-1)*_nVar+2+l]=ratio_densite*v[l];
1365 else if (_limitField[nameOfGroup].bcType==Outlet || (_limitField[nameOfGroup].bcType==InletPressure && q_n>=0)){
1366 return;//For the moment no implicitation, should debug inletPressure
1367 //Boundary state quantities
1368 double v[_Ndim], v2;
1369 VecGetValues(_primitiveVars, _nVar, _idm, _Vj);
1370 double concentration=_Vj[0];
1371 double pression=_limitField[nameOfGroup].p;
1372 double T=_Vj[_nVar-1];
1373 double rho_v=_fluides[0]->getDensity(pression,T);
1374 double rho_l=_fluides[1]->getDensity(pression,T);
1375 if(fabs(concentration*rho_l+(1-concentration)*rho_v)<_precision)
1377 cout<<"rhov= "<<rho_v<<", rhol= "<<rho_l<<endl;
1378 cout<<"concentration*rho_l+(1-concentration)*rho_v= "<<concentration*rho_l+(1-concentration)*rho_v<<endl;
1379 *_runLogFile<<"concentration*rho_l+(1-concentration)*rho_v= "<<concentration*rho_l+(1-concentration)*rho_v<<endl;
1380 _runLogFile->close();
1381 throw CdmathException("DriftModel::jacobian: Outlet, impossible to compute mixture density, division by zero");
1384 double rho_ext=rho_v*rho_l/(concentration*rho_l+(1-concentration)*rho_v);
1385 double rho_int=_externalStates[0];
1386 double e_v=_fluides[0]->getInternalEnergy(T,rho_v);
1387 double e_l=_fluides[1]->getInternalEnergy(T,rho_l);
1388 double ratio_densite=rho_ext/rho_int;
1389 for(k=0;k<_Ndim;k++){
1393 double E_m = concentration*e_v+(1-concentration)*e_l + v2/2;//total energy
1395 //Temperature differential
1396 double C_v= _fluides[0]->constante("cv");
1397 double C_l= _fluides[1]->constante("cv");
1398 double ev0=_fluides[0]->getInternalEnergy(0,rho_v);//Corriger
1399 double el0=_fluides[1]->getInternalEnergy(0,rho_l);//Corriger
1401 double omega=concentration*C_v/(rho_v*e_v)+(1-concentration)*C_l/(rho_l*e_l);
1402 double rhomem=_externalStates[0]*(concentration*e_v+(1-concentration)*e_l);
1403 double m_v=concentration*rho_ext, m_l=rho_ext-m_v;
1404 double alpha_v=m_v/rho_v,alpha_l=1-alpha_v;
1405 double denom=m_v *C_v+m_l* C_l;
1407 double khi=(m_v*(ev0*C_l-el0*C_v)-rhomem*C_l)/(denom*denom);
1408 double ksi=(rho_ext*(el0*C_v-ev0*C_l)+rhomem*(C_l-C_v))/(denom*denom);
1409 double kappa=1/denom;
1411 double CoeffCol1 = rho_int*rho_int*omega*(khi+kappa*v2/2)+ratio_densite;//The +ratio_densite helps filling the lines other than total mass
1412 double CoeffCol2 = rho_int*rho_int*omega*ksi;
1413 double CoeffCol3et4 = rho_int*rho_int*omega*kappa;
1418 for(k=0;k<_Ndim;k++)
1419 _Jcb[2+k]=CoeffCol3et4*v[k];
1420 _Jcb[_nVar-1]=-CoeffCol3et4;
1422 for(k=0;k<_nVar;k++)
1423 _Jcb[_nVar+k]=_Jcb[k]*concentration;
1425 for(int l=0;l<_Ndim;l++)
1426 for(k=0;k<_nVar;k++)
1427 _Jcb[(2+l)*_nVar+k]=_Jcb[k]*v[l];
1429 for(k=0;k<_nVar;k++)
1430 _Jcb[(_nVar-1)*_nVar+k]=_Jcb[k]*E_m;
1432 //adding the remaining diagonal term
1433 for(k=0;k<_nVar;k++)
1434 _Jcb[k*_nVar+k]+=ratio_densite;
1437 else if (_limitField[nameOfGroup].bcType!=Neumann)
1439 cout << "group named "<<nameOfGroup << " : unknown boundary condition" << endl;
1440 _runLogFile->close();
1441 throw CdmathException("DriftModel::jacobian: The boundary condition is not recognised: neither inlet, inltPressure, outlet, wall, InnerWall, nor neumann");
1445 //calcule la jacobienne pour une CL de diffusion
1446 void DriftModel::jacobianDiff(const int &j, string nameOfGroup)
1448 if(_verbose && _nbTimeStep%_freqSave ==0)
1449 cout<<"Jacobienne condition limite diffusion bord "<< nameOfGroup<<endl;
1452 double v2=0,cd = 1,cn=0,p0, gamma;
1455 for(k=0; k<_nVar*_nVar;k++)
1458 if (_limitField[nameOfGroup].bcType==Wall || _limitField[nameOfGroup].bcType==InnerWall){
1460 else if (_limitField[nameOfGroup].bcType==Outlet || _limitField[nameOfGroup].bcType==Neumann
1461 ||_limitField[nameOfGroup].bcType==Inlet || _limitField[nameOfGroup].bcType==InletPressure)
1463 for(k=0;k<_nVar;k++)
1464 _JcbDiff[k*_nVar+k]=1;
1467 cout << "group named "<<nameOfGroup << " : unknown boundary condition" << endl;
1468 _runLogFile->close();
1469 throw CdmathException("DriftModel::jacobianDiff: This boundary condition is not recognised");
1474 void DriftModel::computeScaling(double maxvp)
1476 // if(_spaceScheme!=staggered)
1482 for(int q=2; q<_nVar-1; q++)
1484 _blockDiag[q]=1./maxvp;//
1485 _invBlockDiag[q]= maxvp;//1.;//
1487 _blockDiag[_nVar - 1]=1/(maxvp*maxvp);//1
1488 _invBlockDiag[_nVar - 1]= 1./_blockDiag[_nVar - 1] ;// 1.;//
1491 else{//_spaceScheme==staggered
1492 _blockDiag[0]=maxvp*maxvp;
1493 _invBlockDiag[0]=1/_blockDiag[0];
1494 _blockDiag[1]=maxvp*maxvp;
1495 _invBlockDiag[1]=1/_blockDiag[1];
1496 for(int q=2; q<_nVar-1; q++)
1499 _invBlockDiag[q]= 1;//1.;//
1501 _blockDiag[_nVar - 1]=1;//1
1502 _invBlockDiag[_nVar - 1]= 1./_blockDiag[_nVar - 1] ;// 1.;//
1506 Vector DriftModel::computeExtendedPrimState(double *V)
1508 Vector Vext(7+2*_Ndim);
1510 double C=V[0], P=V[1], T=V[_nVar-1];
1511 double rho_v=_fluides[0]->getDensity(P,T);
1512 double rho_l=_fluides[1]->getDensity(P,T);
1513 double e_v=_fluides[0]->getInternalEnergy(T,rho_v);
1514 double e_l=_fluides[1]->getInternalEnergy(T,rho_l);
1515 if(fabs(rho_l*C+rho_v*(1-C))<_precision)
1517 cout<<"rhov= "<<rho_v<<", rhol= "<<rho_l<<", concentration= "<<C<<endl;
1518 *_runLogFile<<"DriftModel::computeExtendedPrimState: impossible to compute void fraction, division by zero"<<endl;
1519 _runLogFile->close();
1520 throw CdmathException("DriftModel::computeExtendedPrimState: impossible to compute void fraction, division by zero");
1523 _externalStates[0]=rho_v*rho_l/(C*rho_l+(1-C)*rho_v);
1524 double alpha_v=rho_l*C/(rho_l*C+rho_v*(1-C)), alpha_l=1-alpha_v;
1525 double rho_m=alpha_v*rho_v+alpha_l*rho_l;
1526 double h_m=(alpha_v*rho_v*e_v+alpha_l*rho_l*e_l+P)/rho_m;
1528 for(int i=0;i<_Ndim;i++)
1530 Vector u_r=relative_velocity(C, Vit,rho_m);
1534 for(int i=0;i<_Ndim;i++)
1537 Vext(3+_Ndim)=alpha_v;
1538 for(int i=0;i<_Ndim;i++)
1539 Vext(4+_Ndim+i)=u_r(i);
1540 Vext((4+2*_Ndim))=rho_v;
1541 Vext((5+2*_Ndim))=rho_l;
1542 Vext(6+2*_Ndim)=h_m;
1548 void DriftModel::primToCons(const double *P, const int &i, double *W, const int &j){
1549 double concentration=P[i*_nVar];
1550 double pression=P[i*_nVar+1];
1551 double temperature=P[i*_nVar+_nVar-1];
1552 double rho_v=_fluides[0]->getDensity(pression,temperature);
1553 double rho_l=_fluides[1]->getDensity(pression,temperature);
1554 double e_v = _fluides[0]->getInternalEnergy(temperature,rho_v);
1555 double e_l = _fluides[1]->getInternalEnergy(temperature,rho_l);
1556 if(fabs(concentration*rho_l+(1-concentration)*rho_v)<_precision)
1558 cout<<"rhov= "<<rho_v<<", rhol= "<<rho_l<<endl;
1559 cout<<"concentration*rho_l+(1-concentration)*rho_v= "<<concentration*rho_l+(1-concentration)*rho_v<<endl;
1560 *_runLogFile<<"concentration*rho_l+(1-concentration)*rho_v= "<<concentration*rho_l+(1-concentration)*rho_v<<endl;
1561 *_runLogFile<<"DriftModel::primToCons: impossible to compute mixture density, division by zero"<<endl;
1562 _runLogFile->close();
1563 throw CdmathException("DriftModel::primToCons: impossible to compute mixture density, division by zero");
1565 W[j*(_nVar)] =(rho_v*rho_l/(concentration*rho_l+(1-concentration)*rho_v))*_porosityField(j);//phi*rho
1566 W[j*(_nVar)+1] =W[j*(_nVar)] *concentration;//phi *rho*c_v
1567 for(int k=0; k<_Ndim; k++)
1568 W[j*_nVar+(k+2)] = W[j*(_nVar)] *P[i*_nVar+(k+2)];//phi*rho*u
1570 W[j*_nVar+_nVar-1] = W[j*(_nVar)+1]* e_v + (W[j*(_nVar)]-W[j*(_nVar)+1])* e_l;//phi rhom e_m
1571 for(int k=0; k<_Ndim; k++)
1572 W[j*_nVar+_nVar-1] += W[j*_nVar]*P[i*_nVar+(k+2)]*P[i*_nVar+(k+2)]*0.5;//phi rhom e_m + phi rho u*u
1574 if(_verbose && _nbTimeStep%_freqSave ==0){
1575 cout<<"DriftModel::primToCons: T="<<temperature<<", pression= "<<pression<<endl;
1576 cout<<"rhov= "<<rho_v<<", rhol= "<<rho_l<<", rhom= "<<W[j*(_nVar)]<<" e_v="<<e_v<<" e_l="<<e_l<<endl;
1579 void DriftModel::primToConsJacobianMatrix(double *V)
1581 double concentration=V[0];
1582 double pression=V[1];
1583 double temperature=V[_nVar-1];
1584 double vitesse[_Ndim];
1585 for(int idim=0;idim<_Ndim;idim++)
1586 vitesse[idim]=V[2+idim];
1588 for(int idim=0;idim<_Ndim;idim++)
1589 v2+=vitesse[idim]*vitesse[idim];
1591 double rho_v=_fluides[0]->getDensity(pression,temperature);
1592 double rho_l=_fluides[1]->getDensity(pression,temperature);
1593 double gamma_v=_fluides[0]->constante("gamma");
1594 double gamma_l=_fluides[1]->constante("gamma");
1595 double q_v=_fluides[0]->constante("q");
1596 double q_l=_fluides[1]->constante("q");
1598 double rho=concentration*rho_v+(1-concentration)*rho_l;;
1600 for(int k=0;k<_nVar*_nVar; k++)
1601 _primToConsJacoMat[k]=0;
1603 if( !_useDellacherieEOS)
1605 StiffenedGas* fluide0=dynamic_cast<StiffenedGas*>(_fluides[0]);
1606 StiffenedGas* fluide1=dynamic_cast<StiffenedGas*>(_fluides[1]);
1607 double e_v = fluide0->getInternalEnergy(temperature);
1608 double e_l = fluide0->getInternalEnergy(temperature);
1609 double cv_v=fluide0->constante("cv");
1610 double cv_l=fluide1->constante("cv");
1611 double e=concentration*e_v+(1-concentration)*e_l;
1614 /******* Masse totale **********/
1615 _primToConsJacoMat[0]=-rho*rho*(1/rho_v-1/rho_l);
1616 _primToConsJacoMat[1]=rho*rho*(
1617 concentration /(rho_v*rho_v*(gamma_v-1)*(e_v-q_v))
1618 +(1-concentration)/(rho_l*rho_l*(gamma_l-1)*(e_l-q_l))
1620 _primToConsJacoMat[_nVar-1]=-rho*rho*(
1621 cv_v* concentration /(rho_v*(e_v-q_v))
1622 +cv_l*(1-concentration)/(rho_l*(e_l-q_l))
1625 /******* Masse partielle **********/
1626 _primToConsJacoMat[_nVar]=rho-concentration*rho*rho*(1/rho_v-1/rho_l);
1627 _primToConsJacoMat[_nVar+1]=concentration*rho*rho*(
1628 concentration /(rho_v*rho_v*(gamma_v-1)*(e_v-q_v))
1629 +(1-concentration)/(rho_l*rho_l*(gamma_l-1)*(e_l-q_l))
1631 _primToConsJacoMat[_nVar+_nVar-1]=-concentration*rho*rho*(
1632 cv_v* concentration /(rho_v*(e_v-q_v))
1633 +cv_l*(1-concentration)/(rho_l*(e_l-q_l))
1635 /******* Quantité de mouvement **********/
1636 for(int idim=0;idim<_Ndim;idim++)
1638 _primToConsJacoMat[2*_nVar+idim*_nVar]=-vitesse[idim]*rho*rho*(1/rho_v-1/rho_l);
1639 _primToConsJacoMat[2*_nVar+idim*_nVar+1]=vitesse[idim]*rho*rho*(
1640 concentration /(rho_v*rho_v*(gamma_v-1)*(e_v-q_v))
1641 +(1-concentration)/(rho_l*rho_l*(gamma_l-1)*(e_l-q_l))
1643 _primToConsJacoMat[2*_nVar+idim*_nVar+2+idim]=rho;
1644 _primToConsJacoMat[2*_nVar+idim*_nVar+_nVar-1]=-vitesse[idim]*rho*rho*(
1645 cv_v* concentration /(rho_v*(e_v-q_v))
1646 +cv_l*(1-concentration)/(rho_l*(e_l-q_l))
1649 /******* Energie totale **********/
1650 _primToConsJacoMat[(_nVar-1)*_nVar]=rho*(e_v-e_l)-E*rho*rho*(1/rho_v-1/rho_l);
1651 _primToConsJacoMat[(_nVar-1)*_nVar+1]=E*rho*rho*(
1652 concentration /(rho_v*rho_v*(gamma_v-1)*(e_v-q_v))
1653 +(1-concentration)/(rho_l*rho_l*(gamma_l-1)*(e_l-q_l))
1655 for(int idim=0;idim<_Ndim;idim++)
1656 _primToConsJacoMat[(_nVar-1)*_nVar+2+idim]=rho*vitesse[idim];
1657 _primToConsJacoMat[(_nVar-1)*_nVar+_nVar-1]=rho*(cv_v*concentration + cv_l*(1-concentration))
1658 -rho*rho*E*( cv_v* concentration /(rho_v*(e_v-q_v))
1659 +cv_l*(1-concentration)/(rho_l*(e_l-q_l)));
1661 else if(_useDellacherieEOS)
1663 StiffenedGasDellacherie* fluide0=dynamic_cast<StiffenedGasDellacherie*>(_fluides[0]);
1664 StiffenedGasDellacherie* fluide1=dynamic_cast<StiffenedGasDellacherie*>(_fluides[1]);
1665 double h_v=fluide0->getEnthalpy(temperature);
1666 double h_l=fluide1->getEnthalpy(temperature);
1667 double h=concentration*h_v+(1-concentration)*h_l;
1669 double cp_v=fluide0->constante("cp");
1670 double cp_l=fluide1->constante("cp");
1672 /******* Masse totale **********/
1673 _primToConsJacoMat[0]=-rho*rho*(1/rho_v-1/rho_l);
1674 _primToConsJacoMat[1]=rho*rho*(
1675 gamma_v* concentration /(rho_v*rho_v*(gamma_v-1)*(h_v-q_v))
1676 +gamma_l*(1-concentration)/(rho_l*rho_l*(gamma_l-1)*(h_l-q_l))
1678 _primToConsJacoMat[_nVar-1]=-rho*rho*(
1679 cp_v* concentration /(rho_v*(h_v-q_v))
1680 +cp_l*(1-concentration)/(rho_l*(h_l-q_l))
1683 /******* Masse partielle **********/
1684 _primToConsJacoMat[_nVar]=rho-concentration*rho*rho*(1/rho_v-1/rho_l);
1685 _primToConsJacoMat[_nVar+1]=concentration*rho*rho*(
1686 gamma_v* concentration /(rho_v*rho_v*(gamma_v-1)*(h_v-q_v))
1687 +gamma_l*(1-concentration)/(rho_l*rho_l*(gamma_l-1)*(h_l-q_l))
1689 _primToConsJacoMat[_nVar+_nVar-1]=-concentration*rho*rho*(
1690 cp_v* concentration /(rho_v*(h_v-q_v))
1691 +cp_l*(1-concentration)/(rho_l*(h_l-q_l))
1693 /******* Quantité de mouvement **********/
1694 for(int idim=0;idim<_Ndim;idim++)
1696 _primToConsJacoMat[2*_nVar+idim*_nVar]=-vitesse[idim]*rho*rho*(1/rho_v-1/rho_l);
1697 _primToConsJacoMat[2*_nVar+idim*_nVar+1]=vitesse[idim]*rho*rho*(
1698 gamma_v* concentration /(rho_v*rho_v*(gamma_v-1)*(h_v-q_v))
1699 +gamma_l*(1-concentration)/(rho_l*rho_l*(gamma_l-1)*(h_l-q_l))
1701 _primToConsJacoMat[2*_nVar+idim*_nVar+2+idim]=rho;
1702 _primToConsJacoMat[2*_nVar+idim*_nVar+_nVar-1]=-vitesse[idim]*rho*rho*(
1703 cp_v* concentration /(rho_v*(h_v-q_v))
1704 +cp_l*(1-concentration)/(rho_l*(h_l-q_l))
1707 /******* Energie totale **********/
1708 _primToConsJacoMat[(_nVar-1)*_nVar]=rho*(h_v-h_l)-H*rho*rho*(1/rho_v-1/rho_l);
1709 _primToConsJacoMat[(_nVar-1)*_nVar+1]=H*rho*rho*(
1710 gamma_v* concentration /(rho_v*rho_v*(gamma_v-1)*(h_v-q_v))
1711 +gamma_l*(1-concentration)/(rho_l*rho_l*(gamma_l-1)*(h_l-q_l))
1713 for(int idim=0;idim<_Ndim;idim++)
1714 _primToConsJacoMat[(_nVar-1)*_nVar+2+idim]=rho*vitesse[idim];
1715 _primToConsJacoMat[(_nVar-1)*_nVar+_nVar-1]=rho*(cp_v*concentration + cp_l*(1-concentration))
1716 -rho*rho*H*(cp_v* concentration /(rho_v*(h_v-q_v))
1717 +cp_l*(1-concentration)/(rho_l*(h_l-q_l)));
1721 *_runLogFile<<"DriftModel::primToConsJacobianMatrix: eos should be StiffenedGas or StiffenedGasDellacherie"<<endl;
1722 _runLogFile->close();
1723 throw CdmathException("DriftModel::primToConsJacobianMatrix: eos should be StiffenedGas or StiffenedGasDellacherie");
1726 if(_verbose && _nbTimeStep%_freqSave ==0)
1728 cout<<" DriftModel::primToConsJacobianMatrix" << endl;
1729 displayVector(_Vi,_nVar," _Vi " );
1730 cout<<"rho_v= "<<rho_v<<", rho_l= "<<rho_l<<endl;
1731 displayMatrix(_primToConsJacoMat,_nVar," Jacobienne primToCons: ");
1735 void DriftModel::consToPrim(const double *Wcons, double* Wprim,double porosity)
1737 double m_v=Wcons[1]/porosity;
1738 double m_l=(Wcons[0]-Wcons[1])/porosity;
1739 _minm1=min(m_v,_minm1);
1740 _minm2=min(m_l,_minm2);
1741 if(m_v<-_precision || m_l<-_precision)
1743 else if(m_v< 0 && m_v>-_precision )
1745 else if(m_l< 0 && m_l>-_precision )
1747 double concentration=m_v/(m_v+m_l);
1749 for(int k=0;k<_Ndim;k++)
1750 q_2 += Wcons[k+2]*Wcons[k+2];
1751 q_2 /= Wcons[0]; //phi rho u²
1753 double rho_m_e_m=(Wcons[_nVar-1] -0.5*q_2)/porosity;
1754 double pression,Temperature;
1756 if(concentration<_precision)
1758 pression=_fluides[1]->getPressure(rho_m_e_m,m_l);
1759 Temperature=_fluides[1]->getTemperatureFromPressure(pression,m_l);
1761 else if(concentration>1-_precision)
1763 pression=_fluides[0]->getPressure(rho_m_e_m,m_v);
1764 Temperature=_fluides[0]->getTemperatureFromPressure(pression,m_v);
1766 else//Hypothèses de saturation
1769 getMixturePressureAndTemperature(concentration, m_v+m_l, rho_m_e_m, pression, Temperature);
1770 //cout<<"Temperature= "<<Temperature<<", pression= "<<pression<<endl;
1771 //throw CdmathException("DriftModel::consToPrim: Apparition de vapeur");
1773 //Seconde approche : on impose la température directement
1774 //Temperature=_Tsat;
1775 //pression=getMixturePressure(concentration, m_v+m_l, Temperature);
1777 //Troisieme approche : on impose les enthalpies de saturation
1778 //double rho_m_h_m= m_v*_hsatv + m_l*_hsatl;
1779 //pression = rho_m_h_m - rho_m_e_m;
1780 //Temperature = getMixtureTemperature(concentration, m_v+m_l, pression);
1784 cout << "pressure = "<< pression << " < 0 " << endl;
1785 *_runLogFile << "pressure = "<< pression << " < 0 " << endl;
1786 cout<<"Vecteur conservatif"<<endl;
1787 *_runLogFile<<"Vecteur conservatif"<<endl;
1788 for(int k=0;k<_nVar;k++){
1789 cout<<Wcons[k]<<endl;
1790 *_runLogFile<<Wcons[k]<<endl;
1792 _runLogFile->close();
1793 throw CdmathException("DriftModel::consToPrim: negative pressure");
1796 Wprim[0]=concentration;
1797 Wprim[1] = pression;
1798 for(int k=0;k<_Ndim;k++)
1799 Wprim[k+2] = Wcons[k+2]/Wcons[0];
1800 Wprim[_nVar-1] = Temperature;
1801 if(_verbose && _nbTimeStep%_freqSave ==0)
1803 cout<<"ConsToPrim Vecteur conservatif"<<endl;
1804 for(int k=0;k<_nVar;k++)
1805 cout<<Wcons[k]<<endl;
1806 cout<<"ConsToPrim Vecteur primitif"<<endl;
1807 for(int k=0;k<_nVar;k++)
1808 cout<<Wprim[k]<<endl;
1812 double DriftModel::getMixturePressure(double c_v, double rhom, double temperature)
1814 double gamma_v=_fluides[0]->constante("gamma");
1815 double gamma_l=_fluides[1]->constante("gamma");
1816 double Pinf_v=_fluides[0]->constante("p0");
1817 double Pinf_l=_fluides[1]->constante("p0");
1818 double q_v=_fluides[0]->constante("q");
1819 double q_l=_fluides[1]->constante("q");
1823 if( !_useDellacherieEOS)
1825 StiffenedGas* fluide0=dynamic_cast<StiffenedGas*>(_fluides[0]);
1826 StiffenedGas* fluide1=dynamic_cast<StiffenedGas*>(_fluides[1]);
1827 double e_v=fluide0->getInternalEnergy(temperature);
1828 double e_l=fluide1->getInternalEnergy(temperature);
1829 b=gamma_v*Pinf_v+gamma_l*Pinf_l -rhom*(c_l*(gamma_l-1)*(e_l-q_l)+c_v*(gamma_v-1)*(e_v-q_v));
1830 c= gamma_v*Pinf_v*gamma_l*Pinf_l
1831 -rhom*(c_l*(gamma_l-1)*(e_l-q_l)*gamma_v*Pinf_v + c_v*(gamma_v-1)*(e_v-q_v)*gamma_l*Pinf_l);
1833 else if(_useDellacherieEOS)
1835 StiffenedGasDellacherie* fluide0=dynamic_cast<StiffenedGasDellacherie*>(_fluides[0]);
1836 StiffenedGasDellacherie* fluide1=dynamic_cast<StiffenedGasDellacherie*>(_fluides[1]);
1837 double h_v=fluide0->getEnthalpy(temperature);
1838 double h_l=fluide1->getEnthalpy(temperature);
1839 b=Pinf_v+Pinf_l -rhom*(c_l*(gamma_l-1)*(h_l-q_l)/gamma_l+c_v*(gamma_v-1)*(h_v-q_v)/gamma_v);
1840 c=Pinf_v*Pinf_l-rhom*(c_l*(gamma_l-1)*(h_l-q_l)*Pinf_v + c_v*(gamma_v-1)*(h_v-q_v)*Pinf_l);
1844 *_runLogFile<<"DriftModel::getMixturePressure: phases must have the same eos"<<endl;
1845 _runLogFile->close();
1846 throw CdmathException("DriftModel::getMixturePressure: phases must have the same eos");
1849 double delta= b*b-4*a*c;
1851 if(_verbose && _nbTimeStep%_freqSave ==0)
1852 cout<<"DriftModel::getMixturePressure: a= "<<a<<", b= "<<b<<", c= "<<c<<", delta= "<<delta<<endl;
1855 *_runLogFile<<"DriftModel::getMixturePressure: cannot compute pressure, delta<0"<<endl;
1856 _runLogFile->close();
1857 throw CdmathException("DriftModel::getMixturePressure: cannot compute pressure, delta<0");
1860 return (-b+sqrt(delta))/(2*a);
1863 void DriftModel::getMixturePressureAndTemperature(double c_v, double rhom, double rhom_em, double& pression, double& temperature)
1865 double gamma_v=_fluides[0]->constante("gamma");
1866 double gamma_l=_fluides[1]->constante("gamma");
1867 double Pinf_v=_fluides[0]->constante("p0");
1868 double Pinf_l=_fluides[1]->constante("p0");
1869 double q_v=_fluides[0]->constante("q");
1870 double q_l=_fluides[1]->constante("q");
1871 double c_l=1-c_v, m_v=c_v*rhom, m_l=rhom-m_v;
1872 double a, b, c, delta;
1874 if( !_useDellacherieEOS)
1876 StiffenedGas* fluide0=dynamic_cast<StiffenedGas*>(_fluides[0]);
1877 StiffenedGas* fluide1=dynamic_cast<StiffenedGas*>(_fluides[1]);
1879 temperature= (rhom_em-m_v*fluide0->getInternalEnergy(0)-m_l*fluide1->getInternalEnergy(0))
1880 /(m_v*fluide0->constante("cv")+m_l*fluide1->constante("cv"));
1882 double e_v=fluide0->getInternalEnergy(temperature);
1883 double e_l=fluide1->getInternalEnergy(temperature);
1885 b=gamma_v*Pinf_v+gamma_l*Pinf_l -rhom*(c_l*(gamma_l-1)*(e_l-q_l)+c_v*(gamma_v-1)*(e_v-q_v));
1886 c= gamma_v*Pinf_v*gamma_l*Pinf_l
1887 -rhom*(c_l*(gamma_l-1)*(e_l-q_l)*gamma_v*Pinf_v + c_v*(gamma_v-1)*(e_v-q_v)*gamma_l*Pinf_l);
1892 *_runLogFile<<"DriftModel::getMixturePressureAndTemperature: cannot compute pressure, delta<0"<<endl;
1893 _runLogFile->close();
1894 throw CdmathException("DriftModel::getMixturePressureAndTemperature: cannot compute pressure, delta<0");
1897 pression= (-b+sqrt(delta))/(2*a);
1900 else if(_useDellacherieEOS)
1902 StiffenedGasDellacherie* fluide0=dynamic_cast<StiffenedGasDellacherie*>(_fluides[0]);
1903 StiffenedGasDellacherie* fluide1=dynamic_cast<StiffenedGasDellacherie*>(_fluides[1]);
1905 double h0_v=fluide0->getEnthalpy(0);
1906 double h0_l=fluide1->getEnthalpy(0);
1907 double cp_v = _fluides[0]->constante("cp");
1908 double cp_l = _fluides[1]->constante("cp");
1909 double denom=m_v*cp_v+m_l*cp_l;
1910 double num_v=cp_v*rhom_em+m_l*(cp_l* h0_v-cp_v* h0_l);
1911 double num_l=cp_l*rhom_em+m_v*(cp_v* h0_l-cp_l* h0_v);
1913 a=gamma_v*gamma_l-(m_v*(gamma_v-1)*gamma_l*cp_v+m_l*(gamma_l-1)*gamma_v*cp_l)/denom;
1914 b=gamma_v*gamma_l*(Pinf_v+Pinf_l)-(m_v*(gamma_v-1)*gamma_l*cp_v*Pinf_l+m_l*(gamma_l-1)*gamma_v*cp_l*Pinf_v)/denom
1915 -m_v*(gamma_v-1)*gamma_l*(num_v/denom -q_v)-m_l*(gamma_l-1)*gamma_v*(num_l/denom -q_l);
1916 c=gamma_v*gamma_l*Pinf_v*Pinf_l
1917 -m_v*(gamma_v-1)*gamma_l*(num_v/denom -q_v)*Pinf_l-m_l*(gamma_l-1)*gamma_v*(num_l/denom -q_l)*Pinf_v;
1923 *_runLogFile<<"DriftModel::getMixturePressureAndTemperature: cannot compute pressure, delta<0"<<endl;
1924 _runLogFile->close();
1925 throw CdmathException("DriftModel::getMixturePressureAndTemperature: cannot compute pressure, delta<0");
1928 pression= (-b+sqrt(delta))/(2*a);
1930 temperature=(rhom_em+pression-m_v*h0_v-m_l*h0_l)/denom;
1934 _runLogFile->close();
1935 throw CdmathException("DriftModel::getMixturePressureAndTemperature: phases must have the same eos");
1939 if(_verbose && _nbTimeStep%_freqSave ==0){
1940 cout<<"DriftModel::getMixturePressureAndTemperature: a= "<<a<<", b= "<<b<<", c= "<<c<<", delta= "<<delta<<endl;
1941 cout<<"pressure= "<<pression<<", temperature= "<<temperature<<endl;
1945 double DriftModel::getMixtureTemperature(double c_v, double rhom, double pression)
1947 double gamma_v=_fluides[0]->constante("gamma");
1948 double gamma_l=_fluides[1]->constante("gamma");
1949 double Pinf_v = _fluides[0]->constante("p0");
1950 double Pinf_l = _fluides[1]->constante("p0");
1951 double q_v=_fluides[0]->constante("q");
1952 double q_l=_fluides[1]->constante("q");
1955 if( !_useDellacherieEOS)
1957 double cv_v = _fluides[0]->constante("cv");
1958 double cv_l = _fluides[1]->constante("cv");
1959 StiffenedGas* fluide0=dynamic_cast<StiffenedGas*>(_fluides[0]);
1960 StiffenedGas* fluide1=dynamic_cast<StiffenedGas*>(_fluides[1]);
1961 double e0_v=fluide0->getInternalEnergy(0);
1962 double e0_l=fluide1->getInternalEnergy(0);
1964 double numerator=(pression+gamma_v*Pinf_v)*(pression+gamma_l*Pinf_l)/rhom
1965 - c_l*(pression+gamma_v*Pinf_v)*(gamma_l-1)*(e0_l-q_v)
1966 - c_v*(pression+gamma_l*Pinf_l)*(gamma_v-1)*(e0_v-q_l);
1967 double denominator= c_l*(pression+gamma_v*Pinf_v)*(gamma_l-1)*cv_l
1968 +c_v*(pression+gamma_l*Pinf_l)*(gamma_v-1)*cv_v;
1969 return numerator/denominator;
1971 else if(_useDellacherieEOS)
1973 double cp_v = _fluides[0]->constante("cp");
1974 double cp_l = _fluides[1]->constante("cp");
1975 StiffenedGasDellacherie* fluide0=dynamic_cast<StiffenedGasDellacherie*>(_fluides[0]);
1976 StiffenedGasDellacherie* fluide1=dynamic_cast<StiffenedGasDellacherie*>(_fluides[1]);
1977 double h0_v=fluide0->getEnthalpy(0);
1978 double h0_l=fluide1->getEnthalpy(0);
1980 double numerator= gamma_v*(pression+Pinf_v)* gamma_l*(pression+Pinf_l)/rhom
1981 - c_l*gamma_v*(pression+Pinf_v)*(gamma_l-1)*(h0_l-q_l)
1982 - c_v*gamma_l*(pression+Pinf_l)*(gamma_v-1)*(h0_v-q_v);
1983 double denominator= c_l*gamma_v*(pression+Pinf_v)*(gamma_l-1)*cp_l
1984 +c_v*gamma_l*(pression+Pinf_l)*(gamma_v-1)*cp_v;
1985 return numerator/denominator;
1989 *_runLogFile<<"DriftModel::getMixtureTemperature: phases must have the same eos"<<endl;
1990 _runLogFile->close();
1991 throw CdmathException("DriftModel::getMixtureTemperature: phases must have the same eos");
1996 void DriftModel::getMixturePressureDerivatives(double m_v, double m_l, double pression, double Tm)
1998 double rho_v=_fluides[0]->getDensity(pression,Tm);
1999 double rho_l=_fluides[1]->getDensity(pression,Tm);
2000 double alpha_v=m_v/rho_v,alpha_l=1-alpha_v;
2001 double gamma_v=_fluides[0]->constante("gamma");
2002 double gamma_l=_fluides[1]->constante("gamma");
2003 double q_v=_fluides[0]->constante("q");
2004 double q_l=_fluides[1]->constante("q");
2005 double temp1, temp2, denom;
2007 if( !_useDellacherieEOS)
2008 {//Classical stiffened gas with linear law e(T)
2009 double cv_v = _fluides[0]->constante("cv");
2010 double cv_l = _fluides[1]->constante("cv");
2012 double e_v=_fluides[0]->getInternalEnergy(Tm, rho_v);//Actually does not depends on rho_v
2013 double e_l=_fluides[1]->getInternalEnergy(Tm, rho_l);//Actually does not depends on rho_l
2015 //On estime les dérivées discrètes de la pression (cf doc)
2016 temp1= m_v* cv_v + m_l* cv_l;
2017 denom=(alpha_v/((gamma_v-1)*rho_v*(e_v-q_v))+alpha_l/((gamma_l-1)*rho_l*(e_l-q_l)))*temp1;
2018 temp2=alpha_v*cv_v/ (e_v-q_v)+alpha_l*cv_l/ (e_l-q_l);
2020 _khi=(temp1/rho_l-e_l*temp2)/denom;
2021 _ksi=(temp1*(rho_l-rho_v)/(rho_v*rho_l)+(e_l-e_v)*temp2)/denom;
2025 else if( _useDellacherieEOS)
2026 {//S. Dellacherie stiffened gas with linear law h(T)
2027 double cp_v = _fluides[0]->constante("cp");
2028 double cp_l = _fluides[1]->constante("cp");
2030 double h_v=_fluides[0]->getEnthalpy(Tm, rho_v);//Actually does not depends on rho_v
2031 double h_l=_fluides[1]->getEnthalpy(Tm, rho_l);//Actually does not depends on rho_l
2033 //On estime les dérivées discrètes de la pression (cf doc)
2034 temp1= m_v* cp_v + m_l* cp_l;
2035 temp2= alpha_v*cp_v/(h_v-q_v)+alpha_l*cp_l/(h_l-q_l);
2036 //denom=temp1*(alpha_v/(pression+Pinf_v)+alpha_l/(pression+Pinf_l))-temp2;
2037 denom=temp1*(alpha_v*gamma_v/((gamma_v-1)*rho_v*(h_v-q_v))+alpha_l*gamma_l/((gamma_l-1)*rho_l*(h_l-q_l)))-temp2;
2039 //On estime les dérivées discrètes de la pression (cf doc)
2040 _khi=(temp1/rho_l - h_l*temp2)/denom;
2041 _ksi=(temp1*(1/rho_v-1/rho_l)+(h_l-h_v)*temp2)/denom;
2045 if(_verbose && _nbTimeStep%_freqSave ==0)
2047 cout<<"rho_l= "<<rho_l<<", temp1= "<<temp1<<", temp2= "<<temp2<<", denom= "<<denom<<endl;
2048 cout<<"khi= "<<_khi<<", ksi= "<<_ksi<<", kappa= "<<_kappa<<endl;
2052 void DriftModel::entropicShift(double* n)//TO do: make sure _Vi and _Vj are well set
2054 //_Vi is (cm, p, u, T)
2055 //_Ui is (rhom, rhom cm, rhom um, rhom Em)
2056 consToPrim(_Ui,_Vi,_porosityi);
2057 double umnl=0, ul_2=0, umnr=0, ur_2=0; //valeur de u.normale et |u|²
2058 for(int i=0;i<_Ndim;i++)
2060 ul_2 += _Vi[2+i]*_Vi[2+i];
2061 umnl += _Vi[2+i]*n[i];//vitesse normale
2062 ur_2 += _Vj[2+i]*_Vj[2+i];
2063 umnr += _Vj[2+i]*n[i];//vitesse normale
2069 double Hm=(_Ui[_nVar-1]+_Vi[1])/rhom;
2070 if(_verbose && _nbTimeStep%_freqSave ==0)
2071 cout<<"Entropic shift left state rhom="<<rhom<<" cm= "<<cm<<"Hm= "<<Hm<<endl;
2072 double Tm=_Vi[_nVar-1];
2073 double hm=Hm-0.5*ul_2;
2074 double m_v=cm*rhom, m_l=rhom-m_v;
2075 double pression=getMixturePressure( cm, rhom, Tm);
2077 getMixturePressureDerivatives( m_v, m_l, pression, Tm);
2079 if(_kappa*hm+_khi+cm*_ksi<0)
2081 *_runLogFile<<"DriftModel::entropicShift : vitesse du son cellule gauche complexe"<<endl;
2082 _runLogFile->close();
2083 throw CdmathException("Valeurs propres jacobienne: vitesse du son complexe");
2085 double aml=sqrt(_kappa*hm+_khi+cm*_ksi);//vitesse du son du melange
2086 //cout<<"_khi= "<<_khi<<", _kappa= "<< _kappa << ", _ksi= "<<_ksi <<", am= "<<am<<endl;
2089 consToPrim(_Uj,_Vj,_porosityj);
2092 Hm=(_Uj[_nVar-1]+_Vj[1])/rhom;
2093 if(_verbose && _nbTimeStep%_freqSave ==0)
2094 cout<<"Entropic shift right state rhom="<<rhom<<" cm= "<<cm<<"Hm= "<<Hm<<endl;
2099 pression=getMixturePressure( cm, rhom, Tm);
2101 getMixturePressureDerivatives( m_v, m_l, pression, Tm);
2103 if(_kappa*hm+_khi+cm*_ksi<0){
2104 *_runLogFile<<"DriftModel::entropicShift: vitesse du son cellule droite complexe"<<endl;
2105 _runLogFile->close();
2106 throw CdmathException("Valeurs propres jacobienne: vitesse du son complexe");
2109 double amr=sqrt(_kappa*hm+_khi+cm*_ksi);//vitesse du son du melange
2110 //cout<<"_khi= "<<_khi<<", _kappa= "<< _kappa << ", _ksi= "<<_ksi <<", am= "<<am<<endl;
2112 _entropicShift[0]=abs(umnl-aml - (umnr-amr));
2113 _entropicShift[1]=abs(umnl - umnr);
2114 _entropicShift[2]=abs(umnl+aml - (umnr+amr));
2116 if(_verbose && _nbTimeStep%_freqSave ==0)
2118 cout<<"Entropic shift left eigenvalues: "<<endl;
2119 cout<<"("<< umnl-aml<< ", " <<umnl<<", "<<umnl+aml << ")";
2121 cout<<"Entropic shift right eigenvalues: "<<endl;
2122 cout<<"("<< umnr-amr<< ", " <<umnr<<", "<<umnr+amr << ")";
2124 cout<<"eigenvalue jumps "<<endl;
2125 cout<< _entropicShift[0] << ", " << _entropicShift[1] << ", "<< _entropicShift[2] <<endl;
2129 Vector DriftModel::convectionFlux(Vector U,Vector V, Vector normale, double porosity){
2130 if(_verbose && _nbTimeStep%_freqSave ==0)
2132 cout<<"DriftModel::convectionFlux start"<<endl;
2133 cout<<"Ucons"<<endl;
2135 cout<<"Vprim"<<endl;
2139 double rho_m=U(0);//phi rhom
2140 double m_v=U(1), m_l=rho_m-m_v;//phi rhom cv, phi rhom cl
2142 for(int i=0;i<_Ndim;i++)
2145 double concentration_v=V(0);
2146 double concentration_l= 1 - concentration_v;
2147 double pression=V(1);
2148 Vector vitesse_melange(_Ndim);
2149 for(int i=0;i<_Ndim;i++)
2150 vitesse_melange(i)=V(2+i);
2151 double Temperature= V(2+_Ndim);//(rho_m_e_m-m_v*internal_energy(0, e0v,c0v,T0)-m_l*internal_energy(0, e0l,c0l,T0))/(m_v*c0v+m_l*c0l);
2153 Vector vitesse_relative=relative_velocity(concentration_v,vitesse_melange,rho_m);
2154 Vector vitesse_v=vitesse_melange+concentration_l*vitesse_relative;
2155 Vector vitesse_l=vitesse_melange-concentration_v*vitesse_relative;
2156 double vitesse_vn=vitesse_v*normale;
2157 double vitesse_ln=vitesse_l*normale;
2158 //double rho_m_e_m=rho_m_E_m -0.5*(q_m*vitesse_melange + rho_m*concentration_v*concentration_l*vitesse_relative*vitesse_relative)
2159 double rho_v=_fluides[0]->getDensity(pression,Temperature);
2160 double rho_l=_fluides[1]->getDensity(pression,Temperature);
2161 double e_v=_fluides[0]->getInternalEnergy(Temperature,rho_v);
2162 double e_l=_fluides[1]->getInternalEnergy(Temperature,rho_l);
2165 F(0)=m_v*vitesse_vn+m_l*vitesse_ln;
2166 F(1)=m_v*vitesse_vn;
2167 for(int i=0;i<_Ndim;i++)
2168 F(2+i)=m_v*vitesse_vn*vitesse_v(i)+m_l*vitesse_ln*vitesse_l(i)+pression*normale(i)*porosity;
2169 F(2+_Ndim)=m_v*(e_v+0.5*vitesse_v*vitesse_v+pression/rho_v)*vitesse_vn+m_l*(e_l+0.5*vitesse_l*vitesse_l+pression/rho_l)*vitesse_ln;
2171 if(_verbose && _nbTimeStep%_freqSave ==0)
2173 cout<<"DriftModel::convectionFlux end"<<endl;
2174 cout<<"Flux F(U,V)"<<endl;
2181 Vector DriftModel::staggeredVFFCFlux()
2183 if(_verbose && _nbTimeStep%_freqSave ==0)
2184 cout<<"DriftModel::staggeredVFFCFlux()start"<<endl;
2186 if(_spaceScheme!=staggered || _nonLinearFormulation!=VFFC)
2188 *_runLogFile<<"DriftModel::staggeredVFFCFlux: staggeredVFFCFlux method should be called only for VFFC formulation and staggered upwinding"<<endl;
2189 _runLogFile->close();
2190 throw CdmathException("DriftModel::staggeredVFFCFlux: staggeredVFFCFlux method should be called only for VFFC formulation and staggered upwinding");
2192 else//_spaceScheme==staggered
2196 double uijn=0, phiqn=0, uin=0, ujn=0;
2197 for(int idim=0;idim<_Ndim;idim++)
2199 uijn+=_vec_normal[idim]*_Uroe[2+idim];//URoe = rho, cm, u, H
2200 uin+=_vec_normal[idim]*_Ui[2+idim];
2201 ujn+=_vec_normal[idim]*_Uj[2+idim];
2203 if( (uin>0 && ujn >0) || (uin>=0 && ujn <=0 && uijn>0) ) // formerly (uijn>_precision)
2205 for(int idim=0;idim<_Ndim;idim++)
2206 phiqn+=_vec_normal[idim]*_Ui[2+idim];//phi rho u n
2208 Fij(1)=_Vi[0]*phiqn;
2209 for(int idim=0;idim<_Ndim;idim++)
2210 Fij(2+idim)=phiqn*_Vi[2+idim]+_Vj[1]*_vec_normal[idim]*_porosityj;
2211 Fij(_nVar-1)=phiqn/_Ui[0]*(_Ui[_nVar-1]+_Vj[1]*sqrt(_porosityj/_porosityi));
2213 else if( (uin<0 && ujn <0) || (uin>=0 && ujn <=0 && uijn<0) ) // formerly (uijn<-_precision)
2215 for(int idim=0;idim<_Ndim;idim++)
2216 phiqn+=_vec_normal[idim]*_Uj[2+idim];//phi rho u n
2218 Fij(1)=_Vj[0]*phiqn;
2219 for(int idim=0;idim<_Ndim;idim++)
2220 Fij(2+idim)=phiqn*_Vj[2+idim]+_Vi[1]*_vec_normal[idim]*_porosityi;
2221 Fij(_nVar-1)=phiqn/_Uj[0]*(_Uj[_nVar-1]+_Vi[1]*sqrt(_porosityi/_porosityj));
2223 else//case (uin<=0 && ujn >=0) or (uin>=0 && ujn <=0 && uijn==0), apply centered scheme
2225 Vector Ui(_nVar), Uj(_nVar), Vi(_nVar), Vj(_nVar), Fi(_nVar), Fj(_nVar);
2226 Vector normale(_Ndim);
2227 for(int i1=0;i1<_Ndim;i1++)
2228 normale(i1)=_vec_normal[i1];
2229 for(int i1=0;i1<_nVar;i1++)
2236 Fi=convectionFlux(Ui,Vi,normale,_porosityi);
2237 Fj=convectionFlux(Uj,Vj,normale,_porosityj);
2238 Fij=(Fi+Fj)/2+_maxvploc*(Ui-Uj)/2;
2240 //case nil velocity on the interface, apply smoothed scheme
2242 double theta=uijn/(.01);
2243 Vector F1(_nVar),F2(_nVar);
2244 for(int idim=0;idim<_Ndim;idim++)
2245 phiqn+=_vec_normal[idim]*_Ui[2+idim];//phi rho u n
2248 for(int idim=0;idim<_Ndim;idim++)
2249 F1(2+idim)=phiqn*_Vi[2+idim]+_Vj[1]*_vec_normal[idim]*_porosityj;
2250 F1(_nVar-1)=phiqn/_Ui[0]*(_Ui[_nVar-1]+_Vj[1]*sqrt(_porosityj/_porosityi));
2252 for(int idim=0;idim<_Ndim;idim++)
2253 phiqn+=_vec_normal[idim]*_Uj[2+idim];//phi rho u n
2256 for(int idim=0;idim<_Ndim;idim++)
2257 F2(2+idim)=phiqn*_Vj[2+idim]+_Vi[1]*_vec_normal[idim]*_porosityi;
2258 F2(_nVar-1)=phiqn/_Uj[0]*(_Uj[_nVar-1]+_Vi[1]*sqrt(_porosityi/_porosityj));
2260 Fij=(1+theta)/2*F1+(1-theta)/2*F2;
2263 if(_verbose && _nbTimeStep%_freqSave ==0)
2265 cout<<"DriftModel::staggeredVFFCFlux() end uijn="<<uijn<<endl;
2272 void DriftModel::staggeredVFFCMatricesConservativeVariables(double u_mn)
2274 if(_verbose && _nbTimeStep%_freqSave ==0)
2275 cout<<"DriftModel::staggeredVFFCMatricesConservativeVariables()"<<endl;
2277 if(_spaceScheme!=staggered || _nonLinearFormulation!=VFFC)
2279 *_runLogFile<<"DriftModel::staggeredVFFCMatricesConservativeVariables: staggeredVFFCMatrices method should be called only for VFFC formulation and staggered upwinding"<<endl;
2280 _runLogFile->close();
2281 throw CdmathException("DriftModel::staggeredVFFCMatricesConservativeVariables: staggeredVFFCMatrices method should be called only for VFFC formulation and staggered upwinding");
2283 else//_spaceScheme==staggered && _nonLinearFormulation==VFFC
2285 double ui_n=0, ui_2=0, uj_n=0, uj_2=0;//vitesse normale et carré du module
2286 double H;//enthalpie totale (expression particulière)
2287 consToPrim(_Ui,_Vi,_porosityi);
2288 consToPrim(_Uj,_Vj,_porosityj);
2290 for(int i=0;i<_Ndim;i++)
2292 ui_2 += _Vi[2+i]*_Vi[2+i];
2293 ui_n += _Vi[2+i]*_vec_normal[i];
2294 uj_2 += _Vj[2+i]*_Vj[2+i];
2295 uj_n += _Vj[2+i]*_vec_normal[i];
2298 double rhomi=_Ui[0]/_porosityi;
2299 double mi_v=_Ui[1]/_porosityi;
2300 double mi_l=rhomi-mi_v;
2303 double Tmi=_Vi[_Ndim+2];
2304 double Emi=_Ui[_Ndim+2]/(rhomi*_porosityi);
2305 double ecini=0.5*ui_2;
2306 double emi=Emi-0.5*ui_2;
2307 double hmi=emi+pi/rhomi;
2310 double rhomj=_Uj[0]/_porosityj;
2311 double mj_v =_Uj[1]/_porosityj;
2312 double mj_l=rhomj-mj_v;
2314 double Tmj=_Vj[_Ndim+2];
2315 double Emj=_Uj[_Ndim+2]/(rhomj*_porosityj);
2316 double ecinj=0.5*uj_2;
2317 double emj=Emj-0.5*uj_2;
2318 double hmj=emj+pj/rhomj;
2320 double Hm;//value depends on the sign of interfacial velocity
2325 if(_verbose && _nbTimeStep%_freqSave ==0)
2326 cout<<"VFFC Staggered state rhomi="<<rhomi<<" cmi= "<<cmi<<" Hm= "<<Hm<<" Tmi= "<<Tmi<<" pj= "<<pj<<endl;
2328 /***********Calcul des valeurs propres ********/
2329 vector<complex<double> > vp_dist(3,0);
2330 getMixturePressureDerivatives( mj_v, mj_l, pj, Tmj);
2331 if(_kappa*hmj+_khi+cmj*_ksi<0){
2332 *_runLogFile<<"staggeredVFFCMatricesConservativeVariables: vitesse du son complexe"<<endl;
2333 _runLogFile->close();
2334 throw CdmathException("staggeredVFFCMatricesConservativeVariables: vitesse du son complexe");
2336 double amj=sqrt(_kappa*hmj+_khi+cmj*_ksi);//vitesse du son du melange
2338 if(_verbose && _nbTimeStep%_freqSave ==0)
2339 cout<<"_khi= "<<_khi<<", _kappa= "<< _kappa << ", _ksi= "<<_ksi <<", amj= "<<amj<<endl;
2341 //On remplit les valeurs propres
2342 vp_dist[0]=ui_n+amj;
2343 vp_dist[1]=ui_n-amj;
2346 _maxvploc=fabs(ui_n)+amj;
2347 if(_maxvploc>_maxvp)
2350 /******** Construction de la matrice A^+ *********/
2351 _AroePlus[0*_nVar+0]=0;
2352 _AroePlus[0*_nVar+1]=0;
2353 for(int i=0;i<_Ndim;i++)
2354 _AroePlus[0*_nVar+2+i]=_vec_normal[i];
2355 _AroePlus[0*_nVar+2+_Ndim]=0;
2356 _AroePlus[1*_nVar+0]=-ui_n*cmi;
2357 _AroePlus[1*_nVar+1]=ui_n;
2358 for(int i=0;i<_Ndim;i++)
2359 _AroePlus[1*_nVar+2+i]=cmi*_vec_normal[i];
2360 _AroePlus[1*_nVar+2+_Ndim]=0;
2361 for(int i=0;i<_Ndim;i++)
2363 _AroePlus[(2+i)*_nVar+0]=-ui_n*_Vi[2+i];
2364 _AroePlus[(2+i)*_nVar+1]=0;
2365 for(int j=0;j<_Ndim;j++)
2366 _AroePlus[(2+i)*_nVar+2+j]=_Vi[2+i]*_vec_normal[j];
2367 _AroePlus[(2+i)*_nVar+2+i]+=ui_n;
2368 _AroePlus[(2+i)*_nVar+2+_Ndim]=0;
2370 _AroePlus[(2+_Ndim)*_nVar+0]=-Hm*ui_n;
2371 _AroePlus[(2+_Ndim)*_nVar+1]=0;
2372 for(int i=0;i<_Ndim;i++)
2373 _AroePlus[(2+_Ndim)*_nVar+2+i]=Hm*_vec_normal[i];
2374 _AroePlus[(2+_Ndim)*_nVar+2+_Ndim]=ui_n;
2376 /******** Construction de la matrice A^- *********/
2377 _AroeMinus[0*_nVar+0]=0;
2378 _AroeMinus[0*_nVar+1]=0;
2379 for(int i=0;i<_Ndim;i++)
2380 _AroeMinus[0*_nVar+2+i]=0;
2381 _AroeMinus[0*_nVar+2+_Ndim]=0;
2382 _AroeMinus[1*_nVar+0]=0;
2383 _AroeMinus[1*_nVar+1]=0;
2384 for(int i=0;i<_Ndim;i++)
2385 _AroeMinus[1*_nVar+2+i]=0;
2386 _AroeMinus[1*_nVar+2+_Ndim]=0;
2387 for(int i=0;i<_Ndim;i++)
2389 _AroeMinus[(2+i)*_nVar+0]=(_khi+_kappa*ecinj)*_vec_normal[i];
2390 _AroeMinus[(2+i)*_nVar+1]=_ksi*_vec_normal[i];
2391 for(int j=0;j<_Ndim;j++)
2392 _AroeMinus[(2+i)*_nVar+2+j]=-_kappa*_vec_normal[i]*_Vj[2+j];
2393 _AroeMinus[(2+i)*_nVar+2+i]+=0;
2394 _AroeMinus[(2+i)*_nVar+2+_Ndim]=_kappa*_vec_normal[i];
2396 _AroeMinus[(2+_Ndim)*_nVar+0]=(_khi+_kappa*ecinj)*ui_n;
2397 _AroeMinus[(2+_Ndim)*_nVar+1]=_ksi*ui_n;
2398 for(int i=0;i<_Ndim;i++)
2399 _AroeMinus[(2+_Ndim)*_nVar+2+i]=-_kappa*uj_n*_Vi[2+i];
2400 _AroeMinus[(2+_Ndim)*_nVar+2+_Ndim]=_kappa*ui_n;
2402 else if(u_mn<-_precision)
2405 if(_verbose && _nbTimeStep%_freqSave ==0)
2406 cout<<"VFFC Staggered state rhomj="<<rhomj<<" cmj= "<<cmj<<" Hm= "<<Hm<<" Tmj= "<<Tmj<<" pi= "<<pi<<endl;
2408 /***********Calcul des valeurs propres ********/
2409 vector<complex<double> > vp_dist(3,0);
2410 getMixturePressureDerivatives( mi_v, mi_l, pi, Tmi);
2411 if(_kappa*hmi+_khi+cmi*_ksi<0)
2413 *_runLogFile<<"staggeredVFFCMatricesConservativeVariables: vitesse du son complexe"<<endl;
2414 _runLogFile->close();
2415 throw CdmathException("staggeredVFFCMatricesConservativeVariables: vitesse du son complexe");
2417 double ami=sqrt(_kappa*hmi+_khi+cmi*_ksi);//vitesse du son du melange
2418 if(_verbose && _nbTimeStep%_freqSave ==0)
2419 cout<<"_khi= "<<_khi<<", _kappa= "<< _kappa << ", _ksi= "<<_ksi <<", ami= "<<ami<<endl;
2421 //On remplit les valeurs propres
2422 vp_dist[0]=uj_n+ami;
2423 vp_dist[1]=uj_n-ami;
2426 _maxvploc=fabs(uj_n)+ami;
2427 if(_maxvploc>_maxvp)
2430 /******** Construction de la matrice A^+ *********/
2431 _AroePlus[0*_nVar+0]=0;
2432 _AroePlus[0*_nVar+1]=0;
2433 for(int i=0;i<_Ndim;i++)
2434 _AroePlus[0*_nVar+2+i]=0;
2435 _AroePlus[0*_nVar+2+_Ndim]=0;
2436 _AroePlus[1*_nVar+0]=0;
2437 _AroePlus[1*_nVar+1]=0;
2438 for(int i=0;i<_Ndim;i++)
2439 _AroePlus[1*_nVar+2+i]=0;
2440 _AroePlus[1*_nVar+2+_Ndim]=0;
2441 for(int i=0;i<_Ndim;i++)
2443 _AroePlus[(2+i)*_nVar+0]=(_khi+_kappa*ecini)*_vec_normal[i];
2444 _AroePlus[(2+i)*_nVar+1]=_ksi*_vec_normal[i];
2445 for(int j=0;j<_Ndim;j++)
2446 _AroePlus[(2+i)*_nVar+2+j]=-_kappa*_vec_normal[i]*_Vi[2+j];
2447 _AroePlus[(2+i)*_nVar+2+i]+=0;
2448 _AroePlus[(2+i)*_nVar+2+_Ndim]=_kappa*_vec_normal[i];
2450 _AroePlus[(2+_Ndim)*_nVar+0]=(_khi+_kappa*ecini)*ui_n;
2451 _AroePlus[(2+_Ndim)*_nVar+1]=_ksi*ui_n;
2452 for(int i=0;i<_Ndim;i++)
2453 _AroePlus[(2+_Ndim)*_nVar+2+i]=-_kappa*uj_n*_Vj[2+i];
2454 _AroePlus[(2+_Ndim)*_nVar+2+_Ndim]=_kappa*ui_n;
2456 /******** Construction de la matrice A^- *********/
2457 _AroeMinus[0*_nVar+0]=0;
2458 _AroeMinus[0*_nVar+1]=0;
2459 for(int i=0;i<_Ndim;i++)
2460 _AroeMinus[0*_nVar+2+i]=_vec_normal[i];
2461 _AroeMinus[0*_nVar+2+_Ndim]=0;
2462 _AroeMinus[1*_nVar+0]=-uj_n*cmj;
2463 _AroeMinus[1*_nVar+1]=uj_n;
2464 for(int i=0;i<_Ndim;i++)
2465 _AroeMinus[1*_nVar+2+i]=cmj*_vec_normal[i];
2466 _AroeMinus[1*_nVar+2+_Ndim]=0;
2467 for(int i=0;i<_Ndim;i++)
2469 _AroeMinus[(2+i)*_nVar+0]=-uj_n*_Vj[2+i];
2470 _AroeMinus[(2+i)*_nVar+1]=0;
2471 for(int j=0;j<_Ndim;j++)
2472 _AroeMinus[(2+i)*_nVar+2+j]=_Vj[2+i]*_vec_normal[j];
2473 _AroeMinus[(2+i)*_nVar+2+i]+=uj_n;
2474 _AroeMinus[(2+i)*_nVar+2+_Ndim]=0;
2476 _AroeMinus[(2+_Ndim)*_nVar+0]=-Hm*uj_n;
2477 _AroeMinus[(2+_Ndim)*_nVar+1]=0;
2478 for(int i=0;i<_Ndim;i++)
2479 _AroeMinus[(2+_Ndim)*_nVar+2+i]=Hm*_vec_normal[i];
2480 _AroeMinus[(2+_Ndim)*_nVar+2+_Ndim]=uj_n;
2482 else//case nil velocity on the interface
2484 double rhom=_Uroe[0];
2486 double Hm=_Uroe[_nVar-1];
2487 double m_v=cm*rhom, m_l=rhom-m_v;
2490 Vector vitesse(_Ndim);
2491 for(int idim=0;idim<_Ndim;idim++)
2493 vitesse[idim]=_Uroe[2+idim];
2494 umn += _Uroe[2+idim]*_vec_normal[idim];//vitesse normale
2495 u_2 += _Uroe[2+idim]*_Uroe[2+idim];
2497 double hm=Hm-0.5*u_2;
2499 if(cm<_precision)//pure liquid
2500 Tm=_fluides[1]->getTemperatureFromEnthalpy(hm,rhom);
2501 else if(cm>1-_precision)
2502 Tm=_fluides[0]->getTemperatureFromEnthalpy(hm,rhom);
2503 else//Hypothèse de saturation
2506 double pression= getMixturePressure(cm, rhom, Tm);
2508 getMixturePressureDerivatives( m_v, m_l, pression, Tm);//EOS is involved to express pressure jump and sound speed
2509 if(_kappa*hm+_khi+cm*_ksi<0){
2510 *_runLogFile<<"Calcul matrice de Roe: vitesse du son complexe"<<endl;
2511 _runLogFile->close();
2512 throw CdmathException("Calcul matrice de Roe: vitesse du son complexe");
2514 double am=sqrt(_kappa*hm+_khi+cm*_ksi);//vitesse du son du melange
2515 _maxvploc=fabs(umn)+am;
2516 if(_maxvploc>_maxvp)
2520 /******* Construction de la matrice A^+ ********/
2521 //A^+=0.5*jacobienne Flux (U_i)+0.5 _maxvploc Id
2523 getMixturePressureDerivatives( mi_v, mi_l, pi, Tmi);
2524 _AroePlus[0*_nVar+0]=0;
2525 _AroePlus[0*_nVar+1]=0;
2526 for(int i=0;i<_Ndim;i++)
2527 _AroePlus[0*_nVar+2+i]=0.5*_vec_normal[i];
2528 _AroePlus[0*_nVar+2+_Ndim]=0;
2529 _AroePlus[1*_nVar+0]=-0.5*ui_n*cmi;
2530 _AroePlus[1*_nVar+1]=0.5*ui_n;
2531 for(int i=0;i<_Ndim;i++)
2532 _AroePlus[1*_nVar+2+i]=0.5*cmi*_vec_normal[i];
2533 _AroePlus[1*_nVar+2+_Ndim]=0;
2534 for(int i=0;i<_Ndim;i++)
2536 _AroePlus[(2+i)*_nVar+0]=0.5*(_khi+_kappa*ecini)*_vec_normal[i]-0.5*ui_n*_Vi[2+i];
2537 _AroePlus[(2+i)*_nVar+1]=0.5*_ksi*_vec_normal[i];
2538 for(int j=0;j<_Ndim;j++)
2539 _AroePlus[(2+i)*_nVar+2+j]=0.5*_Vi[2+i]*_vec_normal[j]-0.5*_kappa*_vec_normal[i]*_Vi[2+j];
2540 _AroePlus[(2+i)*_nVar+2+i]+=0.5*ui_n;
2541 _AroePlus[(2+i)*_nVar+2+_Ndim]=0.5*_kappa*_vec_normal[i];
2543 _AroePlus[(2+_Ndim)*_nVar+0]=-0.5*Hm*ui_n;
2544 _AroePlus[(2+_Ndim)*_nVar+1]=0;
2545 for(int i=0;i<_Ndim;i++)
2546 _AroePlus[(2+_Ndim)*_nVar+2+i]=0.5*Hm*_vec_normal[i];
2547 _AroePlus[(2+_Ndim)*_nVar+2+_Ndim]=0.5*ui_n;
2549 for(int i=0;i<_nVar;i++)
2550 _AroePlus[i*_nVar+i]+=_maxvploc/2;
2551 /****** Construction de la matrice A^- *******/
2552 //A^-=0.5*jacobienne Flux (U_j)-0.5 _maxvploc Id
2554 getMixturePressureDerivatives( mj_v, mj_l, pj, Tmj);
2555 _AroeMinus[0*_nVar+0]=0;
2556 _AroeMinus[0*_nVar+1]=0;
2557 for(int i=0;i<_Ndim;i++)
2558 _AroeMinus[0*_nVar+2+i]=0.5*_vec_normal[i];
2559 _AroeMinus[0*_nVar+2+_Ndim]=0;
2560 _AroeMinus[1*_nVar+0]=-0.5*uj_n*cmj;
2561 _AroeMinus[1*_nVar+1]=0.5*uj_n;
2562 for(int i=0;i<_Ndim;i++)
2563 _AroeMinus[1*_nVar+2+i]=0.5*cmj*_vec_normal[i];
2564 _AroeMinus[1*_nVar+2+_Ndim]=0;
2565 for(int i=0;i<_Ndim;i++)
2567 _AroeMinus[(2+i)*_nVar+0]=0.5*(_khi+_kappa*ecinj)*_vec_normal[i]-0.5*uj_n*_Vj[2+i];
2568 _AroeMinus[(2+i)*_nVar+1]=0.5*_ksi*_vec_normal[i];
2569 for(int j=0;j<_Ndim;j++)
2570 _AroeMinus[(2+i)*_nVar+2+j]=0.5*_Vj[2+i]*_vec_normal[j]-0.5*_kappa*_vec_normal[i]*_Vj[2+j];
2571 _AroeMinus[(2+i)*_nVar+2+i]+=0.5*uj_n;
2572 _AroeMinus[(2+i)*_nVar+2+_Ndim]=0.5*_kappa*_vec_normal[i];
2574 _AroeMinus[(2+_Ndim)*_nVar+0]=-0.5*Hm*uj_n;
2575 _AroeMinus[(2+_Ndim)*_nVar+1]=0;
2576 for(int i=0;i<_Ndim;i++)
2577 _AroeMinus[(2+_Ndim)*_nVar+2+i]=0.5*Hm*_vec_normal[i];
2578 _AroeMinus[(2+_Ndim)*_nVar+2+_Ndim]=0.5*uj_n;
2580 for(int i=0;i<_nVar;i++)
2581 _AroeMinus[i*_nVar+i]-=_maxvploc/2;
2584 double _AroePlusi[_nVar*_nVar], _AroePlusj[_nVar*_nVar];
2585 double _AroeMinusi[_nVar*_nVar], _AroeMinusj[_nVar*_nVar];
2587 //Matrices côté gauche
2589 ****** Construction de la matrice A^+ i (gauche) *******
2590 _AroePlusi[0*_nVar+0]=0;
2591 _AroePlusi[0*_nVar+1]=0;
2592 for(int i=0;i<_Ndim;i++)
2593 _AroePlusi[0*_nVar+2+i]=_vec_normal[i];
2594 _AroePlusi[0*_nVar+2+_Ndim]=0;
2595 _AroePlusi[1*_nVar+0]=-ui_n*cmi;
2596 _AroePlusi[1*_nVar+1]=ui_n;
2597 for(int i=0;i<_Ndim;i++)
2598 _AroePlusi[1*_nVar+2+i]=cmi*_vec_normal[i];
2599 _AroePlusi[1*_nVar+2+_Ndim]=0;
2600 for(int i=0;i<_Ndim;i++)
2602 _AroePlusi[(2+i)*_nVar+0]=-ui_n*_Vi[2+i];
2603 _AroePlusi[(2+i)*_nVar+1]=0;
2604 for(int j=0;j<_Ndim;j++)
2605 _AroePlusi[(2+i)*_nVar+2+j]=_Vi[2+i]*_vec_normal[j];
2606 _AroePlusi[(2+i)*_nVar+2+i]+=ui_n;
2607 _AroePlusi[(2+i)*_nVar+2+_Ndim]=0;
2609 _AroePlusi[(2+_Ndim)*_nVar+0]=-Hm*ui_n;
2610 _AroePlusi[(2+_Ndim)*_nVar+1]=0;
2611 for(int i=0;i<_Ndim;i++)
2612 _AroePlusi[(2+_Ndim)*_nVar+2+i]=Hm*_vec_normal[i];
2613 _AroePlusi[(2+_Ndim)*_nVar+2+_Ndim]=ui_n;
2615 ****** Construction de la matrice A^- i (coté gauche)*******
2616 _AroeMinusi[0*_nVar+0]=0;
2617 _AroeMinusi[0*_nVar+1]=0;
2618 for(int i=0;i<_Ndim;i++)
2619 _AroeMinusi[0*_nVar+2+i]=0;
2620 _AroeMinusi[0*_nVar+2+_Ndim]=0;
2621 _AroeMinusi[1*_nVar+0]=0;
2622 _AroeMinusi[1*_nVar+1]=0;
2623 for(int i=0;i<_Ndim;i++)
2624 _AroeMinusi[1*_nVar+2+i]=0;
2625 _AroeMinusi[1*_nVar+2+_Ndim]=0;
2626 for(int i=0;i<_Ndim;i++)
2628 _AroeMinusi[(2+i)*_nVar+0]=(_khi+_kappa*ecinj)*_vec_normal[i];
2629 _AroeMinusi[(2+i)*_nVar+1]=_ksi*_vec_normal[i];
2630 for(int j=0;j<_Ndim;j++)
2631 _AroeMinusi[(2+i)*_nVar+2+j]=-_kappa*_vec_normal[i]*_Vj[2+j];
2632 _AroeMinusi[(2+i)*_nVar+2+i]+=0;
2633 _AroeMinusi[(2+i)*_nVar+2+_Ndim]=_kappa*_vec_normal[i];
2635 _AroeMinusi[(2+_Ndim)*_nVar+0]=(_khi+_kappa*ecinj)*ui_n;
2636 _AroeMinusi[(2+_Ndim)*_nVar+1]=_ksi*ui_n;
2637 for(int i=0;i<_Ndim;i++)
2638 _AroeMinusi[(2+_Ndim)*_nVar+2+i]=-_kappa*uj_n*_Vi[2+i];
2639 _AroeMinusi[(2+_Ndim)*_nVar+2+_Ndim]=_kappa*ui_n;
2641 //Matrices côté droit
2644 ****** Construction de la matrice A^+ j (coté droit)*******
2645 _AroePlusj[0*_nVar+0]=0;
2646 _AroePlusj[0*_nVar+1]=0;
2647 for(int i=0;i<_Ndim;i++)
2648 _AroePlusj[0*_nVar+2+i]=0;
2649 _AroePlusj[0*_nVar+2+_Ndim]=0;
2650 _AroePlusj[1*_nVar+0]=0;
2651 _AroePlusj[1*_nVar+1]=0;
2652 for(int i=0;i<_Ndim;i++)
2653 _AroePlusj[1*_nVar+2+i]=0;
2654 _AroePlusj[1*_nVar+2+_Ndim]=0;
2655 for(int i=0;i<_Ndim;i++)
2657 _AroePlusj[(2+i)*_nVar+0]=(_khi+_kappa*ecini)*_vec_normal[i];
2658 _AroePlusj[(2+i)*_nVar+1]=_ksi*_vec_normal[i];
2659 for(int j=0;j<_Ndim;j++)
2660 _AroePlusj[(2+i)*_nVar+2+j]=-_kappa*_vec_normal[i]*_Vi[2+j];
2661 _AroePlusj[(2+i)*_nVar+2+i]+=0;
2662 _AroePlusj[(2+i)*_nVar+2+_Ndim]=_kappa*_vec_normal[i];
2664 _AroePlusj[(2+_Ndim)*_nVar+0]=(_khi+_kappa*ecini)*ui_n;
2665 _AroePlusj[(2+_Ndim)*_nVar+1]=_ksi*ui_n;
2666 for(int i=0;i<_Ndim;i++)
2667 _AroePlusj[(2+_Ndim)*_nVar+2+i]=-_kappa*uj_n*_Vj[2+i];
2668 _AroePlusj[(2+_Ndim)*_nVar+2+_Ndim]=_kappa*ui_n;
2670 ****** Construction de la matrice A^- j (coté droit)*******
2671 _AroeMinusj[0*_nVar+0]=0;
2672 _AroeMinusj[0*_nVar+1]=0;
2673 for(int i=0;i<_Ndim;i++)
2674 _AroeMinusj[0*_nVar+2+i]=_vec_normal[i];
2675 _AroeMinusj[0*_nVar+2+_Ndim]=0;
2676 _AroeMinusj[1*_nVar+0]=-uj_n*cmj;
2677 _AroeMinusj[1*_nVar+1]=uj_n;
2678 for(int i=0;i<_Ndim;i++)
2679 _AroeMinusj[1*_nVar+2+i]=cmj*_vec_normal[i];
2680 _AroeMinusj[1*_nVar+2+_Ndim]=0;
2681 for(int i=0;i<_Ndim;i++)
2683 _AroeMinusj[(2+i)*_nVar+0]=-uj_n*_Vj[2+i];
2684 _AroeMinusj[(2+i)*_nVar+1]=0;
2685 for(int j=0;j<_Ndim;j++)
2686 _AroeMinusj[(2+i)*_nVar+2+j]=_Vj[2+i]*_vec_normal[j];
2687 _AroeMinusj[(2+i)*_nVar+2+i]+=uj_n;
2688 _AroeMinusj[(2+i)*_nVar+2+_Ndim]=0;
2690 _AroeMinusj[(2+_Ndim)*_nVar+0]=-Hm*uj_n;
2691 _AroeMinusj[(2+_Ndim)*_nVar+1]=0;
2692 for(int i=0;i<_Ndim;i++)
2693 _AroeMinusj[(2+_Ndim)*_nVar+2+i]=Hm*_vec_normal[i];
2694 _AroeMinusj[(2+_Ndim)*_nVar+2+_Ndim]=uj_n;
2696 double theta=u_mn/(.01);
2697 for(int i=0; i<_nVar*_nVar;i++)
2699 _AroePlus[i] =(1+theta)/2*_AroePlusi[i] +(1-theta)/2*_AroePlusj[i];
2700 _AroeMinus[i]=(1+theta)/2*_AroeMinusi[i]+(1-theta)/2*_AroeMinusj[i];
2705 for(int i=0; i<_nVar*_nVar;i++)
2707 _Aroe[i]=_AroePlus[i]+_AroeMinus[i];
2708 _absAroe[i]=_AroePlus[i]-_AroeMinus[i];
2711 if(_timeScheme==Implicit)
2712 for(int i=0; i<_nVar*_nVar;i++)
2714 _AroeMinusImplicit[i] = _AroeMinus[i];
2715 _AroePlusImplicit[i] = _AroePlus[i];
2719 void DriftModel::staggeredVFFCMatricesPrimitiveVariables(double u_mn)
2721 if(_verbose && _nbTimeStep%_freqSave ==0)
2722 cout<<"DriftModel::staggeredVFFCMatricesPrimitiveVariables()"<<endl;
2724 if(_spaceScheme!=staggered || _nonLinearFormulation!=VFFC)
2726 *_runLogFile<< "DriftModel::staggeredVFFCMatricesPrimitiveVariables: staggeredVFFCMatricesPrimitiveVariables method should be called only for VFFC formulation and staggered upwinding" << endl;
2727 _runLogFile->close();
2728 throw CdmathException("DriftModel::staggeredVFFCMatricesPrimitiveVariables: staggeredVFFCMatricesPrimitiveVariables method should be called only for VFFC formulation and staggered upwinding");
2730 else//_spaceScheme==staggered && _nonLinearFormulation==VFFC
2732 //Calls to getDensityDerivatives needed
2733 double ui_n=0, ui_2=0, uj_n=0, uj_2=0;//vitesse normale et carré du module
2734 double H;//enthalpie totale (expression particulière)
2735 consToPrim(_Ui,_Vi,_porosityi);
2736 consToPrim(_Uj,_Vj,_porosityj);
2738 for(int i=0;i<_Ndim;i++)
2740 ui_2 += _Vi[2+i]*_Vi[2+i];
2741 ui_n += _Vi[2+i]*_vec_normal[i];
2742 uj_2 += _Vj[2+i]*_Vj[2+i];
2743 uj_n += _Vj[2+i]*_vec_normal[i];
2746 double rhomi=_Ui[0]/_porosityi;
2747 double mi_v=_Ui[1]/_porosityi;
2748 double mi_l=rhomi-mi_v;
2751 double Tmi=_Vi[_Ndim+2];
2752 double Emi=_Ui[_Ndim+2]/(rhomi*_porosityi);
2753 double ecini=0.5*ui_2;
2754 double emi=Emi-0.5*ui_2;
2755 double hmi=emi+pi/rhomi;
2756 double rho_vi=_fluides[0]->getDensity(pi,Tmi);
2757 double rho_li=_fluides[1]->getDensity(pi,Tmi);
2760 double rhomj=_Uj[0]/_porosityj;
2761 double mj_v =_Uj[1]/_porosityj;
2762 double mj_l=rhomj-mj_v;
2764 double Tmj=_Vj[_Ndim+2];
2765 double Emj=_Uj[_Ndim+2]/(rhomj*_porosityj);
2766 double ecinj=0.5*uj_2;
2767 double emj=Emj-0.5*uj_2;
2768 double hmj=emj+pj/rhomj;
2769 double rho_vj=_fluides[0]->getDensity(pj,Tmj);
2770 double rho_lj=_fluides[1]->getDensity(pj,Tmj);
2772 double gamma_v=_fluides[0]->constante("gamma");
2773 double gamma_l=_fluides[1]->constante("gamma");
2774 double q_v=_fluides[0]->constante("q");
2775 double q_l=_fluides[1]->constante("q");
2777 if(fabs(u_mn)>_precision)//non zero velocity on the interface
2779 if( !_useDellacherieEOS)
2781 StiffenedGas* fluide0=dynamic_cast<StiffenedGas*>(_fluides[0]);
2782 StiffenedGas* fluide1=dynamic_cast<StiffenedGas*>(_fluides[1]);
2783 double cv_v=fluide0->constante("cv");
2784 double cv_l=fluide1->constante("cv");
2786 double e_vi=_fluides[0]->getInternalEnergy(Tmi,rho_vi);
2787 double e_li=_fluides[1]->getInternalEnergy(Tmi,rho_li);
2788 double e_vj=_fluides[0]->getInternalEnergy(Tmj,rho_vj);
2789 double e_lj=_fluides[1]->getInternalEnergy(Tmj,rho_lj);
2793 if(_verbose && _nbTimeStep%_freqSave ==0)
2794 cout<<"VFFC Staggered state rhomi="<<rhomi<<" cmi= "<<cmi<<" Emi= "<<Emi<<" Tmi= "<<Tmi<<" pj= "<<pj<<endl;
2796 /***********Calcul des valeurs propres ********/
2797 vector<complex<double> > vp_dist(3,0);
2798 getMixturePressureDerivatives( mj_v, mj_l, pj, Tmj);
2799 if(_kappa*hmj+_khi+cmj*_ksi<0)
2801 *_runLogFile<<"staggeredVFFCMatricesPrimitiveVariables: vitesse du son complexe"<<endl;
2802 _runLogFile->close();
2803 throw CdmathException("staggeredVFFCMatricesPrimitiveVariables: vitesse du son complexe");
2805 double amj=sqrt(_kappa*hmj+_khi+cmj*_ksi);//vitesse du son du melange
2807 if(_verbose && _nbTimeStep%_freqSave ==0)
2808 cout<<"_khi= "<<_khi<<", _kappa= "<< _kappa << ", _ksi= "<<_ksi <<", amj= "<<amj<<endl;
2810 //On remplit les valeurs propres
2811 vp_dist[0]=ui_n+amj;
2812 vp_dist[1]=ui_n-amj;
2815 _maxvploc=fabs(ui_n)+amj;
2816 if(_maxvploc>_maxvp)
2819 /******** Construction de la matrice A^+ *********/
2820 _AroePlusImplicit[0*_nVar+0]=-rhomi*rhomi*ui_n*(1/rho_vi-1/rho_li);
2821 _AroePlusImplicit[0*_nVar+1]= rhomi*rhomi*ui_n*(cmi/(rho_vi*rho_vi*(gamma_v-1)*(e_vi-q_v))+(1-cmi)/(rho_li*rho_li*(gamma_l-1)*(e_li-q_l)));
2822 for(int i=0;i<_Ndim;i++)
2823 _AroePlusImplicit[0*_nVar+2+i]=rhomi*_vec_normal[i];
2824 _AroePlusImplicit[0*_nVar+2+_Ndim]=-rhomi*rhomi*ui_n*(cv_v*cmi/(rho_vi*(e_vi-q_v))+cv_l*(1-cmi)/(rho_li*(e_li-q_l)));
2825 _AroePlusImplicit[1*_nVar+0]=rhomi*ui_n-cmi*rhomi*rhomi*ui_n*(1/rho_vi-1/rho_li);
2826 _AroePlusImplicit[1*_nVar+1]=-cmi*rhomi*rhomi*ui_n*(cmi/(rho_vi*rho_vi*(gamma_v-1)*(e_vi-q_v))+(1-cmi)/(rho_li*rho_li*(gamma_l-1)*(e_li-q_l)));
2827 for(int i=0;i<_Ndim;i++)
2828 _AroePlusImplicit[1*_nVar+2+i]=cmi*rhomi*_vec_normal[i];
2829 _AroePlusImplicit[1*_nVar+2+_Ndim]=-cmi*rhomi*rhomi*ui_n*(cv_v*cmi/(rho_vi*(e_vi-q_v))+cv_l*(1-cmi)/(rho_li*(e_li-q_l)));
2830 for(int i=0;i<_Ndim;i++)
2832 _AroePlusImplicit[(2+i)*_nVar+0]=-rhomi*rhomi*ui_n*(1/rho_vi-1/rho_li)*_Vi[2+i];
2833 _AroePlusImplicit[(2+i)*_nVar+1]=rhomi*rhomi*ui_n*(cmi/(rho_vi*rho_vi*(gamma_v-1)*(e_vi-q_v))+(1-cmi)/(rho_li*rho_li*(gamma_l-1)*(e_li-q_l)))*_Vi[2+i];
2834 for(int j=0;j<_Ndim;j++)
2835 _AroePlusImplicit[(2+i)*_nVar+2+j]=rhomi*_Vi[2+i]*_vec_normal[j];
2836 _AroePlusImplicit[(2+i)*_nVar+2+i]+=rhomi*ui_n;
2837 _AroePlusImplicit[(2+i)*_nVar+2+_Ndim]=-rhomi*rhomi*ui_n*(cv_v*cmi/(rho_vi*(e_vi-q_v))+(1-cmi)/(rho_li*(e_li-q_l)))*_Vi[2+i];
2839 _AroePlusImplicit[(2+_Ndim)*_nVar+0]=ui_n*(rhomi*(e_vi-e_li)-Emi*rhomi*rhomi*(1/rho_vi-1/rho_li));
2840 _AroePlusImplicit[(2+_Ndim)*_nVar+1]=rhomi*rhomi*Emi*(cmi/(rho_vi*rho_vi*(gamma_v-1)*(e_vi-q_v))+(1-cmi)/(rho_li*rho_li*(gamma_l-1)*(e_li-q_l)));
2841 for(int i=0;i<_Ndim;i++)
2842 _AroePlusImplicit[(2+_Ndim)*_nVar+2+i]=(rhomi*Emi+pj)*_vec_normal[i] + ui_n*rhomi*_Vi[2+i];
2843 _AroePlusImplicit[(2+_Ndim)*_nVar+2+_Ndim]=ui_n*(rhomi*(cv_v*cmi+cv_l*(1-cmi))-Emi*rhomi*rhomi*(cv_v*cmi/(rho_vi*(e_vi-q_v))+cv_l*(1-cmi)/(rho_li*(e_li-q_l))));
2846 /******** Construction de la matrice A^- *********/
2847 _AroeMinusImplicit[0*_nVar+0]=0;
2848 _AroeMinusImplicit[0*_nVar+1]=0;
2849 for(int i=0;i<_Ndim;i++)
2850 _AroeMinusImplicit[0*_nVar+2+i]=0;
2851 _AroeMinusImplicit[0*_nVar+2+_Ndim]=0;
2852 _AroeMinusImplicit[1*_nVar+0]=0;
2853 _AroeMinusImplicit[1*_nVar+1]=0;
2854 for(int i=0;i<_Ndim;i++)
2855 _AroeMinusImplicit[1*_nVar+2+i]=0;
2856 _AroeMinusImplicit[1*_nVar+2+_Ndim]=0;
2857 for(int i=0;i<_Ndim;i++)
2859 _AroeMinusImplicit[(2+i)*_nVar+0]=0;
2860 _AroeMinusImplicit[(2+i)*_nVar+1]=_vec_normal[i];
2861 for(int j=0;j<_Ndim;j++)
2862 _AroeMinusImplicit[(2+i)*_nVar+2+j]=0;
2863 _AroeMinusImplicit[(2+i)*_nVar+2+i]+=0;
2864 _AroeMinusImplicit[(2+i)*_nVar+2+_Ndim]=0;
2866 _AroeMinusImplicit[(2+_Ndim)*_nVar+0]=0;
2867 _AroeMinusImplicit[(2+_Ndim)*_nVar+1]=ui_n;
2868 for(int i=0;i<_Ndim;i++)
2869 _AroeMinusImplicit[(2+_Ndim)*_nVar+2+i]=0;
2870 _AroeMinusImplicit[(2+_Ndim)*_nVar+2+_Ndim]=0;
2872 else if(u_mn<-_precision)
2874 if(_verbose && _nbTimeStep%_freqSave ==0)
2875 cout<<"VFFC Staggered state rhomj="<<rhomj<<" cmj= "<<cmj<<" Emj= "<<Emj<<" Tmj= "<<Tmj<<" pi= "<<pi<<endl;
2877 /***********Calcul des valeurs propres ********/
2878 vector<complex<double> > vp_dist(3,0);
2879 getMixturePressureDerivatives( mi_v, mi_l, pi, Tmi);
2880 if(_kappa*hmi+_khi+cmi*_ksi<0)
2882 *_runLogFile<<"staggeredVFFCMatricesPrimitiveVariables: vitesse du son complexe"<<endl;
2883 throw CdmathException("staggeredVFFCMatricesPrimitiveVariables: vitesse du son complexe");
2885 double ami=sqrt(_kappa*hmi+_khi+cmi*_ksi);//vitesse du son du melange
2886 if(_verbose && _nbTimeStep%_freqSave ==0)
2887 cout<<"_khi= "<<_khi<<", _kappa= "<< _kappa << ", _ksi= "<<_ksi <<", ami= "<<ami<<endl;
2889 //On remplit les valeurs propres
2890 vp_dist[0]=uj_n+ami;
2891 vp_dist[1]=uj_n-ami;
2894 _maxvploc=fabs(uj_n)+ami;
2895 if(_maxvploc>_maxvp)
2898 /******** Construction de la matrice A^+ *********/
2899 _AroePlusImplicit[0*_nVar+0]=0;
2900 _AroePlusImplicit[0*_nVar+1]=0;
2901 for(int i=0;i<_Ndim;i++)
2902 _AroePlusImplicit[0*_nVar+2+i]=0;
2903 _AroePlusImplicit[0*_nVar+2+_Ndim]=0;
2904 _AroePlusImplicit[1*_nVar+0]=0;
2905 _AroePlusImplicit[1*_nVar+1]=0;
2906 for(int i=0;i<_Ndim;i++)
2907 _AroePlusImplicit[1*_nVar+2+i]=0;
2908 _AroePlusImplicit[1*_nVar+2+_Ndim]=0;
2909 for(int i=0;i<_Ndim;i++)
2911 _AroePlusImplicit[(2+i)*_nVar+0]=0;
2912 _AroePlusImplicit[(2+i)*_nVar+1]=_vec_normal[i];
2913 for(int j=0;j<_Ndim;j++)
2914 _AroePlusImplicit[(2+i)*_nVar+2+j]=0;
2915 _AroePlusImplicit[(2+i)*_nVar+2+i]+=0;
2916 _AroePlusImplicit[(2+i)*_nVar+2+_Ndim]=0;
2918 _AroePlusImplicit[(2+_Ndim)*_nVar+0]=0;
2919 _AroePlusImplicit[(2+_Ndim)*_nVar+1]=uj_n;
2920 for(int i=0;i<_Ndim;i++)
2921 _AroePlusImplicit[(2+_Ndim)*_nVar+2+i]=0;
2922 _AroePlusImplicit[(2+_Ndim)*_nVar+2+_Ndim]=0;
2924 /******** Construction de la matrice A^- *********/
2925 _AroeMinusImplicit[0*_nVar+0]=-rhomj*rhomj*uj_n*(1/rho_vj-1/rho_lj);
2926 _AroeMinusImplicit[0*_nVar+1]= rhomj*rhomj*uj_n*(cmj/(rho_vj*rho_vj*(gamma_v-1)*(e_vj-q_v))+(1-cmj)/(rho_lj*rho_lj*(gamma_l-1)*(e_lj-q_l)));
2927 for(int i=0;i<_Ndim;i++)
2928 _AroeMinusImplicit[0*_nVar+2+i]=rhomj*_vec_normal[i];
2929 _AroeMinusImplicit[0*_nVar+2+_Ndim]=-rhomj*rhomj*uj_n*(cv_v*cmj/(rho_vj*(e_vj-q_v))+cv_l*(1-cmj)/(rho_lj*(e_lj-q_l)));
2930 _AroeMinusImplicit[1*_nVar+0]=rhomj*uj_n-cmj*rhomj*rhomj*uj_n*(1/rho_vj-1/rho_lj);
2931 _AroeMinusImplicit[1*_nVar+1]=-cmj*rhomj*rhomj*uj_n*(cmj/(rho_vj*rho_vj*(gamma_v-1)*(e_vj-q_v))+(1-cmj)/(rho_lj*rho_lj*(gamma_l-1)*(e_lj-q_l)));
2932 for(int i=0;i<_Ndim;i++)
2933 _AroeMinusImplicit[1*_nVar+2+i]=cmj*rhomj*_vec_normal[i];
2934 _AroeMinusImplicit[1*_nVar+2+_Ndim]=-cmj*rhomj*rhomj*uj_n*(cv_v*cmj/(rho_vj*(e_vj-q_v))+cv_l*(1-cmj)/(rho_lj*(e_lj-q_l)));
2935 for(int i=0;i<_Ndim;i++)
2937 _AroeMinusImplicit[(2+i)*_nVar+0]=-rhomj*rhomj*uj_n*(1/rho_vj-1/rho_lj)*_Vj[2+i];
2938 _AroeMinusImplicit[(2+i)*_nVar+1]=rhomj*rhomj*uj_n*(cmj/(rho_vj*rho_vj*(gamma_v-1)*(e_vj-q_v))+(1-cmj)/(rho_lj*rho_lj*(gamma_l-1)*(e_lj-q_l)))*_Vj[2+i];
2939 for(int j=0;j<_Ndim;j++)
2940 _AroeMinusImplicit[(2+i)*_nVar+2+j]=rhomj*_Vj[2+i]*_vec_normal[j];
2941 _AroeMinusImplicit[(2+i)*_nVar+2+i]+=rhomj*uj_n;
2942 _AroeMinusImplicit[(2+i)*_nVar+2+_Ndim]=-rhomj*rhomj*uj_n*(cv_v*cmj/(rho_vj*(e_vj-q_v))+cv_l*(1-cmj)/(rho_lj*(e_lj-q_l)))*_Vj[2+i];
2944 _AroeMinusImplicit[(2+_Ndim)*_nVar+0]=uj_n*(rhomj*(e_vj-e_lj)-Emj*rhomj*rhomj*(1/rho_vj-1/rho_lj));
2945 _AroeMinusImplicit[(2+_Ndim)*_nVar+1]=rhomj*rhomj*Emj*(cmj/(rho_vj*rho_vj*(gamma_v-1)*(e_vj-q_v))+(1-cmj)/(rho_lj*rho_lj*(gamma_l-1)*(e_lj-q_l)));
2946 for(int i=0;i<_Ndim;i++)
2947 _AroeMinusImplicit[(2+_Ndim)*_nVar+2+i]=(rhomj*Emj+pi)*_vec_normal[i] + uj_n*rhomj*_Vj[2+i];
2948 _AroeMinusImplicit[(2+_Ndim)*_nVar+2+_Ndim]=uj_n*(rhomj*(cv_v*cmj+cv_l*(1-cmj))-Emj*rhomj*rhomj*(cv_v*cmj/(rho_vj*(e_vj-q_v))+cv_l*(1-cmj)/(rho_lj*(e_lj-q_l))));
2952 *_runLogFile<< "DriftModel::staggeredVFFCMatricesPrimitiveVariables: velocity umn should be non zero" << endl;
2953 _runLogFile->close();
2954 throw CdmathException("DriftModel::staggeredVFFCMatricesPrimitiveVariables: velocity umn should be non zero");
2957 else if(_useDellacherieEOS)
2959 StiffenedGasDellacherie* fluide0=dynamic_cast<StiffenedGasDellacherie*>(_fluides[0]);
2960 StiffenedGasDellacherie* fluide1=dynamic_cast<StiffenedGasDellacherie*>(_fluides[1]);
2961 double cp_v=fluide0->constante("cp");
2962 double cp_l=fluide1->constante("cp");
2964 double h_vi=_fluides[0]->getEnthalpy(Tmi,rho_vi);
2965 double h_li=_fluides[1]->getEnthalpy(Tmi,rho_li);
2966 double h_vj=_fluides[0]->getEnthalpy(Tmj,rho_vj);
2967 double h_lj=_fluides[1]->getEnthalpy(Tmj,rho_lj);
2968 double Hmi=Emi+pi/rho_vi;
2969 double Hmj=Emj+pj/rho_vj;
2973 if(_verbose && _nbTimeStep%_freqSave ==0)
2974 cout<<"VFFC Staggered state rhomi="<<rhomi<<" cmi= "<<cmi<<" Hmi= "<<Hmi<<" Tmi= "<<Tmi<<" pj= "<<pj<<endl;
2976 /***********Calcul des valeurs propres ********/
2977 vector<complex<double> > vp_dist(3,0);
2978 getMixturePressureDerivatives( mj_v, mj_l, pj, Tmj);
2979 if(_kappa*hmj+_khi+cmj*_ksi<0)
2981 *_runLogFile<<"staggeredVFFCMatricesPrimitiveVariables: vitesse du son complexe"<<endl;
2982 throw CdmathException("staggeredVFFCMatricesPrimitiveVariables: vitesse du son complexe");
2984 double amj=sqrt(_kappa*hmj+_khi+cmj*_ksi);//vitesse du son du melange
2986 if(_verbose && _nbTimeStep%_freqSave ==0)
2987 cout<<"_khi= "<<_khi<<", _kappa= "<< _kappa << ", _ksi= "<<_ksi <<", amj= "<<amj<<endl;
2989 //On remplit les valeurs propres
2990 vp_dist[0]=ui_n+amj;
2991 vp_dist[1]=ui_n-amj;
2994 _maxvploc=fabs(ui_n)+amj;
2995 if(_maxvploc>_maxvp)
2998 /******** Construction de la matrice A^+ *********/
2999 _AroePlusImplicit[0*_nVar+0]=-rhomi*rhomi*ui_n*(1/rho_vi-1/rho_li);
3000 _AroePlusImplicit[0*_nVar+1]= rhomi*rhomi*ui_n*(gamma_v*cmi/(rho_vi*rho_vi*(gamma_v-1)*(h_vi-q_v))+gamma_l*(1-cmi)/(rho_li*rho_li*(gamma_l-1)*(h_li-q_l)));
3001 for(int i=0;i<_Ndim;i++)
3002 _AroePlusImplicit[0*_nVar+2+i]=rhomi*_vec_normal[i];
3003 _AroePlusImplicit[0*_nVar+2+_Ndim]=-rhomi*rhomi*ui_n*(cp_v*cmi/(rho_vi*(h_vi-q_v))+cp_l*(1-cmi)/(rho_li*(h_li-q_l)));
3004 _AroePlusImplicit[1*_nVar+0]=rhomi*ui_n-cmi*rhomi*rhomi*ui_n*(1/rho_vi-1/rho_li);
3005 _AroePlusImplicit[1*_nVar+1]=-cmi*rhomi*rhomi*ui_n*(gamma_v*cmi/(rho_vi*rho_vi*(gamma_v-1)*(h_vi-q_v))+gamma_l*(1-cmi)/(rho_li*rho_li*(gamma_l-1)*(h_li-q_l)));
3006 for(int i=0;i<_Ndim;i++)
3007 _AroePlusImplicit[1*_nVar+2+i]=cmi*rhomi*_vec_normal[i];
3008 _AroePlusImplicit[1*_nVar+2+_Ndim]=-cmi*rhomi*rhomi*ui_n*(cp_v*cmi/(rho_vi*(h_vi-q_v))+cp_l*(1-cmi)/(rho_li*(h_li-q_l)));
3009 for(int i=0;i<_Ndim;i++)
3011 _AroePlusImplicit[(2+i)*_nVar+0]=-rhomi*rhomi*ui_n*(1/rho_vi-1/rho_li)*_Vi[2+i];
3012 _AroePlusImplicit[(2+i)*_nVar+1]=rhomi*rhomi*ui_n*(gamma_v*cmi/(rho_vi*rho_vi*(gamma_v-1)*(h_vi-q_v))+gamma_l*(1-cmi)/(rho_li*rho_li*(gamma_l-1)*(h_li-q_l)))*_Vi[2+i];
3013 for(int j=0;j<_Ndim;j++)
3014 _AroePlusImplicit[(2+i)*_nVar+2+j]=rhomi*_Vi[2+i]*_vec_normal[j];
3015 _AroePlusImplicit[(2+i)*_nVar+2+i]+=rhomi*ui_n;
3016 _AroePlusImplicit[(2+i)*_nVar+2+_Ndim]=-rhomi*rhomi*ui_n*(cp_v*cmi/(rho_vi*(h_vi-q_v))+cp_l*(1-cmi)/(rho_li*(h_li-q_l)))*_Vi[2+i];
3018 _AroePlusImplicit[(2+_Ndim)*_nVar+0]=ui_n*(rhomi*(h_vi-h_li)-Hmi*rhomi*rhomi*(1/rho_vi-1/rho_li));
3019 _AroePlusImplicit[(2+_Ndim)*_nVar+1]=rhomi*rhomi*Hmi*(gamma_v*cmi/(rho_vi*rho_vi*(gamma_v-1)*(h_vi-q_v))+gamma_l*(1-cmi)/(rho_li*rho_li*(gamma_l-1)*(h_li-q_l)));
3020 for(int i=0;i<_Ndim;i++)
3021 _AroePlusImplicit[(2+_Ndim)*_nVar+2+i]=(rhomi*Emi+pj)*_vec_normal[i] + ui_n*rhomi*_Vi[2+i];
3022 _AroePlusImplicit[(2+_Ndim)*_nVar+2+_Ndim]=ui_n*(rhomi*(cp_v*cmi+cp_l*(1-cmi))-Hmi*rhomi*rhomi*(cp_v*cmi/(rho_vi*(h_vi-q_v))+cp_l*(1-cmi)/(rho_li*(h_li-q_l))));
3025 /******** Construction de la matrice A^- *********/
3026 _AroeMinusImplicit[0*_nVar+0]=0;
3027 _AroeMinusImplicit[0*_nVar+1]=0;
3028 for(int i=0;i<_Ndim;i++)
3029 _AroeMinusImplicit[0*_nVar+2+i]=0;
3030 _AroeMinusImplicit[0*_nVar+2+_Ndim]=0;
3031 _AroeMinusImplicit[1*_nVar+0]=0;
3032 _AroeMinusImplicit[1*_nVar+1]=0;
3033 for(int i=0;i<_Ndim;i++)
3034 _AroeMinusImplicit[1*_nVar+2+i]=0;
3035 _AroeMinusImplicit[1*_nVar+2+_Ndim]=0;
3036 for(int i=0;i<_Ndim;i++)
3038 _AroeMinusImplicit[(2+i)*_nVar+0]=0;
3039 _AroeMinusImplicit[(2+i)*_nVar+1]=_vec_normal[i];
3040 for(int j=0;j<_Ndim;j++)
3041 _AroeMinusImplicit[(2+i)*_nVar+2+j]=0;
3042 _AroeMinusImplicit[(2+i)*_nVar+2+i]+=0;
3043 _AroeMinusImplicit[(2+i)*_nVar+2+_Ndim]=0;
3045 _AroeMinusImplicit[(2+_Ndim)*_nVar+0]=0;
3046 _AroeMinusImplicit[(2+_Ndim)*_nVar+1]=ui_n;
3047 for(int i=0;i<_Ndim;i++)
3048 _AroeMinusImplicit[(2+_Ndim)*_nVar+2+i]=0;
3049 _AroeMinusImplicit[(2+_Ndim)*_nVar+2+_Ndim]=0;
3051 else if(u_mn<-_precision)
3053 if(_verbose && _nbTimeStep%_freqSave ==0)
3054 cout<<"VFFC Staggered state rhomj="<<rhomj<<" cmj= "<<cmj<<" Hmj= "<<Hmj<<" Tmj= "<<Tmj<<" pi= "<<pi<<endl;
3056 /***********Calcul des valeurs propres ********/
3057 vector<complex<double> > vp_dist(3,0);
3058 getMixturePressureDerivatives( mi_v, mi_l, pi, Tmi);
3059 if(_kappa*hmi+_khi+cmi*_ksi<0)
3061 *_runLogFile<<"staggeredVFFCMatricesPrimitiveVariables: vitesse du son complexe"<<endl;
3062 throw CdmathException("staggeredVFFCMatricesPrimitiveVariables: vitesse du son complexe");
3064 double ami=sqrt(_kappa*hmi+_khi+cmi*_ksi);//vitesse du son du melange
3065 if(_verbose && _nbTimeStep%_freqSave ==0)
3066 cout<<"_khi= "<<_khi<<", _kappa= "<< _kappa << ", _ksi= "<<_ksi <<", ami= "<<ami<<endl;
3068 //On remplit les valeurs propres
3069 vp_dist[0]=uj_n+ami;
3070 vp_dist[1]=uj_n-ami;
3073 _maxvploc=fabs(uj_n)+ami;
3074 if(_maxvploc>_maxvp)
3077 /******** Construction de la matrice A^+ *********/
3078 _AroePlusImplicit[0*_nVar+0]=0;
3079 _AroePlusImplicit[0*_nVar+1]=0;
3080 for(int i=0;i<_Ndim;i++)
3081 _AroePlusImplicit[0*_nVar+2+i]=0;
3082 _AroePlusImplicit[0*_nVar+2+_Ndim]=0;
3083 _AroePlusImplicit[1*_nVar+0]=0;
3084 _AroePlusImplicit[1*_nVar+1]=0;
3085 for(int i=0;i<_Ndim;i++)
3086 _AroePlusImplicit[1*_nVar+2+i]=0;
3087 _AroePlusImplicit[1*_nVar+2+_Ndim]=0;
3088 for(int i=0;i<_Ndim;i++)
3090 _AroePlusImplicit[(2+i)*_nVar+0]=0;
3091 _AroePlusImplicit[(2+i)*_nVar+1]=_vec_normal[i];
3092 for(int j=0;j<_Ndim;j++)
3093 _AroePlusImplicit[(2+i)*_nVar+2+j]=0;
3094 _AroePlusImplicit[(2+i)*_nVar+2+i]+=0;
3095 _AroePlusImplicit[(2+i)*_nVar+2+_Ndim]=0;
3097 _AroePlusImplicit[(2+_Ndim)*_nVar+0]=0;
3098 _AroePlusImplicit[(2+_Ndim)*_nVar+1]=uj_n;
3099 for(int i=0;i<_Ndim;i++)
3100 _AroePlusImplicit[(2+_Ndim)*_nVar+2+i]=0;
3101 _AroePlusImplicit[(2+_Ndim)*_nVar+2+_Ndim]=0;
3103 /******** Construction de la matrice A^- *********/
3104 _AroeMinusImplicit[0*_nVar+0]=-rhomj*rhomj*uj_n*(1/rho_vj-1/rho_lj);
3105 _AroeMinusImplicit[0*_nVar+1]= rhomj*rhomj*uj_n*(gamma_v*cmj/(rho_vj*rho_vj*(gamma_v-1)*(h_vj-q_v))+gamma_l*(1-cmj)/(rho_lj*rho_lj*(gamma_l-1)*(h_lj-q_l)));
3106 for(int i=0;i<_Ndim;i++)
3107 _AroeMinusImplicit[0*_nVar+2+i]=rhomj*_vec_normal[i];
3108 _AroeMinusImplicit[0*_nVar+2+_Ndim]=-rhomj*rhomj*uj_n*(cp_v*cmj/(rho_vj*(h_vj-q_v))+cp_l*(1-cmj)/(rho_lj*(h_lj-q_l)));
3109 _AroeMinusImplicit[1*_nVar+0]=rhomj*uj_n-cmj*rhomj*rhomj*uj_n*(1/rho_vj-1/rho_lj);
3110 _AroeMinusImplicit[1*_nVar+1]=-cmj*rhomj*rhomj*uj_n*(gamma_v*cmj/(rho_vj*rho_vj*(gamma_v-1)*(h_vj-q_v))+gamma_l*(1-cmj)/(rho_lj*rho_lj*(gamma_l-1)*(h_lj-q_l)));
3111 for(int i=0;i<_Ndim;i++)
3112 _AroeMinusImplicit[1*_nVar+2+i]=cmj*rhomj*_vec_normal[i];
3113 _AroeMinusImplicit[1*_nVar+2+_Ndim]=-cmj*rhomj*rhomj*uj_n*(cp_v*cmj/(rho_vj*(h_vj-q_v))+cp_l*(1-cmj)/(rho_lj*(h_lj-q_l)));
3114 for(int i=0;i<_Ndim;i++)
3116 _AroeMinusImplicit[(2+i)*_nVar+0]=-rhomj*rhomj*uj_n*(1/rho_vj-1/rho_lj)*_Vj[2+i];
3117 _AroeMinusImplicit[(2+i)*_nVar+1]=rhomj*rhomj*uj_n*(gamma_v*cmj/(rho_vj*rho_vj*(gamma_v-1)*(h_vj-q_v))+gamma_l*(1-cmj)/(rho_lj*rho_lj*(gamma_l-1)*(h_lj-q_l)))*_Vj[2+i];
3118 for(int j=0;j<_Ndim;j++)
3119 _AroeMinusImplicit[(2+i)*_nVar+2+j]=rhomj*_Vj[2+i]*_vec_normal[j];
3120 _AroeMinusImplicit[(2+i)*_nVar+2+i]+=rhomj*uj_n;
3121 _AroeMinusImplicit[(2+i)*_nVar+2+_Ndim]=-rhomj*rhomj*uj_n*(cp_v*cmj/(rho_vj*(h_vj-q_v))+cp_l*(1-cmj)/(rho_lj*(h_lj-q_l)))*_Vj[2+i];
3123 _AroeMinusImplicit[(2+_Ndim)*_nVar+0]=uj_n*(rhomj*(h_vj-h_lj)-Hmj*rhomj*rhomj*(1/rho_vj-1/rho_lj));
3124 _AroeMinusImplicit[(2+_Ndim)*_nVar+1]=rhomj*rhomj*Hmj*(gamma_v*cmj/(rho_vj*rho_vj*(gamma_v-1)*(h_vj-q_v))+gamma_l*(1-cmj)/(rho_lj*rho_lj*(gamma_l-1)*(h_lj-q_l)));
3125 for(int i=0;i<_Ndim;i++)
3126 _AroeMinusImplicit[(2+_Ndim)*_nVar+2+i]=(rhomj*Emj+pi)*_vec_normal[i] + uj_n*rhomj*_Vj[2+i];
3127 _AroeMinusImplicit[(2+_Ndim)*_nVar+2+_Ndim]=uj_n*(rhomj*(cp_v*cmj+cp_l*(1-cmj))-Hmj*rhomj*rhomj*(cp_v*cmj/(rho_vj*(h_vj-q_v))+cp_l*(1-cmj)/(rho_lj*(h_lj-q_l))));
3131 *_runLogFile<< "DriftModel::staggeredVFFCMatricesPrimitiveVariables: velocity umn should be non zero" << endl;
3132 _runLogFile->close();
3133 throw CdmathException("DriftModel::staggeredVFFCMatricesPrimitiveVariables: velocity umn should be non zero");
3137 throw CdmathException("DriftModel::staggeredVFFCMatricesPrimitiveVariables: eos should be StiffenedGas or StiffenedGasDellacherie");
3139 else//case nil velocity on the interface, multiply by jacobian matrix
3142 primToConsJacobianMatrix(_Vj);
3143 Poly.matrixProduct(_AroeMinus, _nVar, _nVar, _primToConsJacoMat, _nVar, _nVar, _AroeMinusImplicit);
3144 primToConsJacobianMatrix(_Vi);
3145 Poly.matrixProduct(_AroePlus, _nVar, _nVar, _primToConsJacoMat, _nVar, _nVar, _AroePlusImplicit);
3150 void DriftModel::applyVFRoeLowMachCorrections(bool isBord, string nameOfGroup)
3152 if(_nonLinearFormulation!=VFRoe)
3153 throw CdmathException("DriftModel::applyVFRoeLowMachCorrections: applyVFRoeLowMachCorrections method should be called only for VFRoe formulation");
3154 else//_nonLinearFormulation==VFRoe
3156 if(_spaceScheme==lowMach){
3158 for(int i=0;i<_Ndim;i++)
3159 u_2 += _Uroe[2+i]*_Uroe[2+i];
3161 double c = _maxvploc;//mixture sound speed
3162 double M=sqrt(u_2)/c;//Mach number
3163 _Vij[1]=M*_Vij[1]+(1-M)*(_Vi[1]+_Vj[1])/2;
3165 else if(_spaceScheme==pressureCorrection)
3168 * option 1 : no pressure correction except for inner walls where p^*=p_int
3169 * option 2 : Clerc pressure correction everywhere, except for inner walls where p^*=p_int
3170 * option 3 : Clerc pressure correction everywhere, even for inner walls
3171 * option 4 : Clerc pressure correction only inside, upwind on wall and innerwall boundaries
3172 * option 5 : Clerc pressure correction inside the domain and special pressure correction involving gravity at wall and inner wall boundaries
3174 bool isWall=false, isInnerWall=false;
3177 isWall = (_limitField[nameOfGroup].bcType==Wall);
3178 isInnerWall = (_limitField[nameOfGroup].bcType==InnerWall);
3181 if( (!isInnerWall && _pressureCorrectionOrder==2) || _pressureCorrectionOrder==3 ||
3182 (!isBord && _pressureCorrectionOrder==4) || (!isWall && !isInnerWall && _pressureCorrectionOrder==5) )//Clerc pressure correction
3184 double norm_uij=0, uij_n=0, ui_n=0, uj_n=0;
3185 for(int i=0;i<_Ndim;i++)
3187 norm_uij += _Uroe[2+i]*_Uroe[2+i];
3188 uij_n += _Uroe[2+i]*_vec_normal[i];
3189 ui_n += _Vi[2+i]*_vec_normal[i];
3190 uj_n += _Vj[2+i]*_vec_normal[i];
3192 norm_uij=sqrt(norm_uij);
3193 if(norm_uij>_precision)//avoid division by zero
3194 _Vij[1]=(_Vi[1]+_Vj[1])/2 + uij_n/norm_uij*(_Vj[1]-_Vi[1])/4 - _Uroe[0]*norm_uij*(uj_n-ui_n)/4;
3196 _Vij[1]=(_Vi[1]+_Vj[1])/2 - _Uroe[0]*norm_uij*(uj_n-ui_n)/4;
3198 else if((isWall || isInnerWall) && _pressureCorrectionOrder==5)//correction de pression gravitaire
3200 double g_n=0;//scalar product of gravity and normal vector
3201 for(int i=0;i<_Ndim;i++)
3202 g_n += _GravityField3d[i]*_vec_normal[i];
3203 _Vij[1]=_Vi[1]+ _Ui[0]*g_n/_inv_dxi/2;
3205 else if(isInnerWall && (_pressureCorrectionOrder==1 || _pressureCorrectionOrder==2) )
3208 else if(_spaceScheme==staggered)
3211 for(int i=0;i<_Ndim;i++)
3212 uij_n += _Uroe[1+i]*_vec_normal[i];
3214 if(uij_n>_precision){
3217 for(int i=0;i<_Ndim;i++)
3219 _Vij[_nVar-1]=_Vi[_nVar-1];
3221 else if(uij_n<-_precision){
3224 for(int i=0;i<_Ndim;i++)
3226 _Vij[_nVar-1]=_Vj[_nVar-1];
3229 _Vij[0]=(_Vi[0]+_Vj[0])/2;
3230 _Vij[1]=(_Vi[1]+_Vj[1])/2;
3231 for(int i=0;i<_Ndim;i++)
3232 _Vij[2+i]=(_Vi[2+i]+_Vj[2+i])/2;
3233 _Vij[_nVar-1]=(_Vi[_nVar-1]+_Vj[_nVar-1])/2;
3236 primToCons(_Vij,0,_Uij,0);
3239 void DriftModel::RoeMatrixConservativeVariables(double cm, double umn,double ecin, double Hm,Vector vitesse)
3241 //On remplit la matrice de Roe en variables conservatives : F(u_L)-F(U_R)=Aroe (U_L-U_R)
3244 for(int i=0;i<_Ndim;i++)
3245 _Aroe[0*_nVar+2+i]=_vec_normal[i];
3246 _Aroe[0*_nVar+2+_Ndim]=0;
3247 _Aroe[1*_nVar+0]=-umn*cm;
3248 _Aroe[1*_nVar+1]=umn;
3249 for(int i=0;i<_Ndim;i++)
3250 _Aroe[1*_nVar+2+i]=cm*_vec_normal[i];
3251 _Aroe[1*_nVar+2+_Ndim]=0;
3252 for(int i=0;i<_Ndim;i++)
3254 _Aroe[(2+i)*_nVar+0]=(_khi+_kappa*ecin)*_vec_normal[i]-umn*vitesse[i];
3255 _Aroe[(2+i)*_nVar+1]=_ksi*_vec_normal[i];
3256 for(int j=0;j<_Ndim;j++)
3257 _Aroe[(2+i)*_nVar+2+j]=vitesse[i]*_vec_normal[j]-_kappa*_vec_normal[i]*vitesse[j];
3258 _Aroe[(2+i)*_nVar+2+i]+=umn;
3259 _Aroe[(2+i)*_nVar+2+_Ndim]=_kappa*_vec_normal[i];
3261 _Aroe[(2+_Ndim)*_nVar+0]=(_khi+_kappa*ecin-Hm)*umn;
3262 _Aroe[(2+_Ndim)*_nVar+1]=_ksi*umn;
3263 for(int i=0;i<_Ndim;i++)
3264 _Aroe[(2+_Ndim)*_nVar+2+i]=Hm*_vec_normal[i]-_kappa*umn*vitesse[i];
3265 _Aroe[(2+_Ndim)*_nVar+2+_Ndim]=(_kappa+1)*umn;
3268 void DriftModel::staggeredRoeUpwindingMatrixConservativeVariables(double cm, double umn,double ecin, double Hm, Vector vitesse)
3270 //Calcul de décentrement de type décalé pour formulation Roe
3271 if(abs(umn)>_precision)//non zero velocity at the interface
3273 _absAroe[0*_nVar+0]=0;
3274 _absAroe[0*_nVar+1]=0;
3275 for(int i=0;i<_Ndim;i++)
3276 _absAroe[0*_nVar+2+i]=_vec_normal[i];
3277 _absAroe[0*_nVar+2+_Ndim]=0;
3278 _absAroe[1*_nVar+0]=-umn*cm;
3279 _absAroe[1*_nVar+1]=umn;
3280 for(int i=0;i<_Ndim;i++)
3281 _absAroe[1*_nVar+2+i]=cm*_vec_normal[i];
3282 _absAroe[1*_nVar+2+_Ndim]=0;
3283 for(int i=0;i<_Ndim;i++)
3285 _absAroe[(2+i)*_nVar+0]=-(_khi+_kappa*ecin)*_vec_normal[i]-umn*vitesse[i];
3286 _absAroe[(2+i)*_nVar+1]=-_ksi*_vec_normal[i];
3287 for(int j=0;j<_Ndim;j++)
3288 _absAroe[(2+i)*_nVar+2+j]=vitesse[i]*_vec_normal[j]+_kappa*_vec_normal[i]*vitesse[j];
3289 _absAroe[(2+i)*_nVar+2+i]+=umn;
3290 _absAroe[(2+i)*_nVar+2+_Ndim]=-_kappa*_vec_normal[i];
3292 _absAroe[(2+_Ndim)*_nVar+0]=(-_khi-_kappa*ecin-Hm)*umn;
3293 _absAroe[(2+_Ndim)*_nVar+1]=-_ksi*umn;
3294 for(int i=0;i<_Ndim;i++)
3295 _absAroe[(2+_Ndim)*_nVar+2+i]=Hm*_vec_normal[i]+_kappa*umn*vitesse[i];
3296 _absAroe[(2+_Ndim)*_nVar+2+_Ndim]=(-_kappa+1)*umn;
3301 else // umn<-_precision
3304 for(int i=0; i<_nVar*_nVar;i++)
3305 _absAroe[i] *= signu;
3307 else//umn=0>centered scheme
3309 for(int i=0; i<_nVar*_nVar;i++)
3311 //for(int i=0; i<_nVar;i++)
3312 // _absAroe[i+_nVar*i] =_maxvploc;
3316 void DriftModel::staggeredRoeUpwindingMatrixPrimitiveVariables(double concentration, double rhom, double umn, double Hm, Vector vitesse)
3318 //Not used. Suppress or use in alternative implicitation in primitive variable of the staggered-roe scheme
3319 //Calcul de décentrement de type décalé pour formulation Roe
3320 _absAroeImplicit[0*_nVar+0]= _drho_sur_dcv*umn;
3321 _absAroeImplicit[0*_nVar+1]= _drho_sur_dp *umn;
3322 for(int i=0;i<_Ndim;i++)
3323 _absAroeImplicit[0*_nVar+2+i]=rhom*_vec_normal[i];
3324 _absAroeImplicit[0*_nVar+2+_Ndim]=_drho_sur_dT*umn;
3325 _absAroeImplicit[1*_nVar+0]= _drhocv_sur_dcv*umn;
3326 _absAroeImplicit[1*_nVar+1]= _drhocv_sur_dp *umn;
3327 for(int i=0;i<_Ndim;i++)
3328 _absAroeImplicit[1*_nVar+2+i]=rhom*concentration*_vec_normal[i];
3329 _absAroeImplicit[1*_nVar+2+_Ndim]=_drhocv_sur_dT*umn;
3330 for(int i=0;i<_Ndim;i++)
3332 _absAroeImplicit[(2+i)*_nVar+0]= _drho_sur_dcv*umn*vitesse[i];
3333 _absAroeImplicit[(2+i)*_nVar+1]= _drho_sur_dp *umn*vitesse[i]-_vec_normal[i];
3334 for(int j=0;j<_Ndim;j++)
3335 _absAroeImplicit[(2+i)*_nVar+2+j]=rhom*vitesse[i]*_vec_normal[j];
3336 _absAroeImplicit[(2+i)*_nVar+2+i]+=rhom*umn;
3337 _absAroeImplicit[(2+i)*_nVar+2+_Ndim]=_drho_sur_dT*umn*vitesse[i];
3339 _absAroeImplicit[(2+_Ndim)*_nVar+0]= _drhoE_sur_dcv *umn;
3340 _absAroeImplicit[(2+_Ndim)*_nVar+1]=( _drhoE_sur_dp+1)*umn;
3341 for(int i=0;i<_Ndim;i++)
3342 _absAroeImplicit[(2+_Ndim)*_nVar+2+i]=rhom*(Hm*_vec_normal[i]+umn*vitesse[i]);
3343 _absAroeImplicit[(2+_Ndim)*_nVar+2+_Ndim]=_drhoE_sur_dT*umn;
3346 void DriftModel::convectionMatrixPrimitiveVariables(double concentration, double rhom, double umn, double Hm, Vector vitesse)
3348 //Not used. Suppress or use in alternative implicitation in primitive variable of the staggered-roe scheme
3349 //On remplit la matrice de Roe en variables primitives : F(V_L)-F(V_R)=Aroe (V_L-V_R)
3350 //EOS is more involved with primitive variables
3351 //first call to getDensityDerivatives(double concentration, double pression, double temperature,double v2) needed
3352 _AroeImplicit[0*_nVar+0]=_drho_sur_dcv*umn;
3353 _AroeImplicit[0*_nVar+1]=_drho_sur_dp*umn;
3354 for(int i=0;i<_Ndim;i++)
3355 _AroeImplicit[0*_nVar+2+i]=rhom*_vec_normal[i];
3356 _AroeImplicit[0*_nVar+2+_Ndim]=_drho_sur_dT*umn;
3357 _AroeImplicit[1*_nVar+0]=_drhocv_sur_dcv*umn;
3358 _AroeImplicit[1*_nVar+1]=_drhocv_sur_dp*umn;
3359 for(int i=0;i<_Ndim;i++)
3360 _AroeImplicit[1*_nVar+2+i]=rhom*concentration*_vec_normal[i];
3361 _AroeImplicit[1*_nVar+2+_Ndim]=_drhocv_sur_dT*umn;
3362 for(int i=0;i<_Ndim;i++)
3364 _AroeImplicit[(2+i)*_nVar+0]=_drho_sur_dcv*umn*vitesse[i];
3365 _AroeImplicit[(2+i)*_nVar+1]=_drho_sur_dp *umn*vitesse[i]+_vec_normal[i];
3366 for(int j=0;j<_Ndim;j++)
3367 _AroeImplicit[(2+i)*_nVar+2+j]=rhom*vitesse[i]*_vec_normal[j];
3368 _AroeImplicit[(2+i)*_nVar+2+i]+=rhom*umn;
3369 _AroeImplicit[(2+i)*_nVar+2+_Ndim]=_drho_sur_dT*umn*vitesse[i];
3371 _AroeImplicit[(2+_Ndim)*_nVar+0]= _drhoE_sur_dcv *umn;
3372 _AroeImplicit[(2+_Ndim)*_nVar+1]=(_drhoE_sur_dp+1)*umn;
3373 for(int i=0;i<_Ndim;i++)
3374 _AroeImplicit[(2+_Ndim)*_nVar+2+i]=rhom*(Hm*_vec_normal[i]+umn*vitesse[i]);
3375 _AroeImplicit[(2+_Ndim)*_nVar+2+_Ndim]=_drhoE_sur_dT*umn;
3378 void DriftModel::getDensityDerivatives(double concentration, double pression, double temperature,double v2)
3380 //EOS is more involved with primitive variables
3382 double rho_v=_fluides[0]->getDensity(pression,temperature);
3383 double rho_l=_fluides[1]->getDensity(pression,temperature);
3384 double gamma_v=_fluides[0]->constante("gamma");
3385 double gamma_l=_fluides[1]->constante("gamma");
3386 double q_v=_fluides[0]->constante("q");
3387 double q_l=_fluides[1]->constante("q");
3389 double rho=concentration*rho_v+(1-concentration)*rho_l;;
3391 if( !_useDellacherieEOS)
3393 StiffenedGas* fluide0=dynamic_cast<StiffenedGas*>(_fluides[0]);
3394 StiffenedGas* fluide1=dynamic_cast<StiffenedGas*>(_fluides[1]);
3395 double e_v = fluide0->getInternalEnergy(temperature);
3396 double e_l = fluide1->getInternalEnergy(temperature);
3397 double cv_v=fluide0->constante("cv");
3398 double cv_l=fluide1->constante("cv");
3399 double e=concentration*e_v+(1-concentration)*e_l;
3402 _drho_sur_dcv=-rho*rho*(1/rho_v-1/rho_l);
3403 _drho_sur_dp=rho*rho*(
3404 concentration /(rho_v*rho_v*(gamma_v-1)*(e_v-q_v))
3405 +(1-concentration)/(rho_l*rho_l*(gamma_l-1)*(e_l-q_l))
3407 _drho_sur_dT=-rho*rho*(
3408 cv_v* concentration /(rho_v*(e_v-q_v))
3409 +cv_l*(1-concentration)/(rho_l*(e_l-q_l))
3412 _drhocv_sur_dcv=rho-concentration*rho*rho*(1/rho_v-1/rho_l);
3413 _drhocv_sur_dp=concentration*rho*rho*(
3414 concentration /(rho_v*rho_v*(gamma_v-1)*(e_v-q_v))
3415 +(1-concentration)/(rho_l*rho_l*(gamma_l-1)*(e_l-q_l))
3417 _drhocv_sur_dT=-concentration*rho*rho*(
3418 cv_v* concentration /(rho_v*(e_v-q_v))
3419 +cv_l*(1-concentration)/(rho_l*(e_l-q_l))
3421 _drhoE_sur_dcv=rho*(e_v-e_l)-E*rho*rho*(1/rho_v-1/rho_l);
3422 _drhoE_sur_dp=E*rho*rho*(
3423 concentration /(rho_v*rho_v*(gamma_v-1)*(e_v-q_v))
3424 +(1-concentration)/(rho_l*rho_l*(gamma_l-1)*(e_l-q_l))
3426 _drhoE_sur_dT=rho*(cv_v*concentration + cv_l*(1-concentration))
3427 -rho*rho*E*( cv_v* concentration /(rho_v*(e_v-q_v))
3428 +cv_l*(1-concentration)/(rho_l*(e_l-q_l)));
3430 else if(_useDellacherieEOS)
3432 StiffenedGasDellacherie* fluide0=dynamic_cast<StiffenedGasDellacherie*>(_fluides[0]);
3433 StiffenedGasDellacherie* fluide1=dynamic_cast<StiffenedGasDellacherie*>(_fluides[1]);
3434 double h_v=fluide0->getEnthalpy(temperature);
3435 double h_l=fluide1->getEnthalpy(temperature);
3436 double h=concentration*h_v+(1-concentration)*h_l;
3438 double cp_v=fluide0->constante("cp");
3439 double cp_l=fluide1->constante("cp");
3441 _drho_sur_dcv=-rho*rho*(1/rho_v-1/rho_l);
3442 _drho_sur_dp =rho*rho*(
3443 gamma_v* concentration /(rho_v*rho_v*(gamma_v-1)*(h_v-q_v))
3444 +gamma_l*(1-concentration)/(rho_l*rho_l*(gamma_l-1)*(h_l-q_l))
3446 _drho_sur_dT=-rho*rho*(
3447 cp_v* concentration /(rho_v*(h_v-q_v))
3448 +cp_l*(1-concentration)/(rho_l*(h_l-q_l))
3451 _drhocv_sur_dcv=rho-concentration*rho*rho*(1/rho_v-1/rho_l);
3452 _drhocv_sur_dp= concentration*rho*rho*(
3453 gamma_v* concentration /(rho_v*rho_v*(gamma_v-1)*(h_v-q_v))
3454 +gamma_l*(1-concentration)/(rho_l*rho_l*(gamma_l-1)*(h_l-q_l))
3456 _drhocv_sur_dT=-concentration*rho*rho*(
3457 cp_v* concentration /(rho_v*(h_v-q_v))
3458 +cp_l*(1-concentration)/(rho_l*(h_l-q_l))
3460 _drhoE_sur_dcv=rho*(h_v-h_l)-H*rho*rho*(1/rho_v-1/rho_l);
3461 _drhoE_sur_dp=H*rho*rho*(
3462 gamma_v* concentration /(rho_v*rho_v*(gamma_v-1)*(h_v-q_v))
3463 +gamma_l*(1-concentration)/(rho_l*rho_l*(gamma_l-1)*(h_l-q_l))
3465 _drhoE_sur_dT=rho*(cp_v*concentration + cp_l*(1-concentration))
3466 -rho*rho*H*( cp_v* concentration /(rho_v*(h_v-q_v))
3467 +cp_l*(1-concentration)/(rho_l*(h_l-q_l)));
3470 throw CdmathException("DriftModel::primToConsJacobianMatrix: eos should be StiffenedGas or StiffenedGasDellacherie");
3472 if(_verbose && _nbTimeStep%_freqSave ==0)
3474 cout<<"_drho_sur_dcv= "<<_drho_sur_dcv<<", _drho_sur_dp= "<<_drho_sur_dp<<", _drho_sur_dT= "<<_drho_sur_dT<<endl;
3475 cout<<"_drhocv_sur_dcv= "<<_drhocv_sur_dcv<<", _drhocv_sur_dp= "<<_drhocv_sur_dp<<", _drhocv_sur_dT= "<<_drhocv_sur_dT<<endl;
3476 cout<<"_drhoE_sur_dcv= "<<_drhoE_sur_dcv<<", _drhoE_sur_dp= "<<_drhoE_sur_dp<<", _drhoE_sur_dT= "<<_drhoE_sur_dT<<endl;
3480 void DriftModel::save(){
3481 string prim(_path+"/DriftModelPrim_");
3482 string cons(_path+"/DriftModelCons_");
3483 string allFields(_path+"/");
3486 allFields+=_fileName;
3489 for (long i = 0; i < _Nmailles; i++){
3490 // j = 0 : concentration, j=1 : pressure; j = _nVar - 1: temperature; j = 2,..,_nVar-2: velocity
3491 for (int j = 0; j < _nVar; j++){
3493 VecGetValues(_primitiveVars,1,&Ii,&_VV(i,j));
3496 if(_saveConservativeField){
3497 for (long i = 0; i < _Nmailles; i++){
3498 // j = 0 : total density; j = 1 : vapour density; j = _nVar - 1 : energy j = 2,..,_nVar-2: momentum
3499 for (int j = 0; j < _nVar; j++){
3501 VecGetValues(_conservativeVars,1,&Ii,&_UU(i,j));
3504 _UU.setTime(_time,_nbTimeStep);
3506 _VV.setTime(_time,_nbTimeStep);
3507 // create mesh and component info
3508 if (_nbTimeStep ==0 || _restartWithNewFileName){
3509 if (_restartWithNewFileName)
3510 _restartWithNewFileName=false;
3511 string suppress_previous_runs ="rm -rf *"+_fileName+"_*";
3512 system(suppress_previous_runs.c_str());//Nettoyage des précédents calculs identiques
3514 if(_saveConservativeField){
3515 _UU.setInfoOnComponent(0,"Total_Density");// (kg/m^3)
3517 _UU.setInfoOnComponent(1,"Partial_Density");// (kg/m^3)
3518 _UU.setInfoOnComponent(2,"Momentum_x");// (kg/m^2/s)
3520 _UU.setInfoOnComponent(3,"Momentum_y");// (kg/m^2/s)
3522 _UU.setInfoOnComponent(4,"Momentum_z");// (kg/m^2/s)
3524 _UU.setInfoOnComponent(_nVar-1,"Energy_(J/m^3)");
3539 _VV.setInfoOnComponent(0,"Concentration");
3540 _VV.setInfoOnComponent(1,"Pressure_(Pa)");
3541 _VV.setInfoOnComponent(2,"Velocity_x_(m/s)");
3543 _VV.setInfoOnComponent(3,"Velocity_y_(m/s)");
3545 _VV.setInfoOnComponent(4,"Velocity_z_(m/s)");
3546 _VV.setInfoOnComponent(_nVar-1,"Temperature_(K)");
3560 // do not create mesh
3565 _VV.writeVTK(prim,false);
3568 _VV.writeMED(prim,false);
3574 if(_saveConservativeField){
3578 _UU.writeVTK(cons,false);
3581 _UU.writeMED(cons,false);
3590 for (long i = 0; i < _Nmailles; i++){
3591 // j = 0 : concentration, j=1 : pressure; j = _nVar - 1: temperature; j = 2,..,_nVar-2: velocity
3592 for (int j = 0; j < _Ndim; j++)//On récupère les composantes de vitesse
3594 int Ii = i*_nVar +2+j;
3595 VecGetValues(_primitiveVars,1,&Ii,&_Vitesse(i,j));
3597 for (int j = _Ndim; j < 3; j++)//On met à zero les composantes de vitesse si la dimension est <3
3600 _Vitesse.setTime(_time,_nbTimeStep);
3601 if (_nbTimeStep ==0 || _restartWithNewFileName){
3602 _Vitesse.setInfoOnComponent(0,"Velocity_x_(m/s)");
3603 _Vitesse.setInfoOnComponent(1,"Velocity_y_(m/s)");
3604 _Vitesse.setInfoOnComponent(2,"Velocity_z_(m/s)");
3608 _Vitesse.writeVTK(prim+"_Velocity");
3611 _Vitesse.writeMED(prim+"_Velocity");
3614 _Vitesse.writeCSV(prim+"_Velocity");
3622 _Vitesse.writeVTK(prim+"_Velocity",false);
3625 _Vitesse.writeMED(prim+"_Velocity",false);
3628 _Vitesse.writeCSV(prim+"_Velocity");
3635 double p,Tm,cv,alpha_v,rhom,rho_v,rho_l, m_v, m_l, h_v, h_l, vx,vy,vz;
3637 for (long i = 0; i < _Nmailles; i++){
3639 VecGetValues(_conservativeVars,1,&Ii,&rhom);
3641 VecGetValues(_primitiveVars,1,&Ii,&cv);
3643 VecGetValues(_primitiveVars,1,&Ii,&p);
3644 Ii = i*_nVar +_nVar-1;
3645 VecGetValues(_primitiveVars,1,&Ii,&Tm);
3647 VecGetValues(_primitiveVars,1,&Ii,&vx);
3651 VecGetValues(_primitiveVars,1,&Ii,&vy);
3654 VecGetValues(_primitiveVars,1,&Ii,&vz);
3658 rho_v=_fluides[0]->getDensity(p,Tm);
3659 rho_l=_fluides[1]->getDensity(p,Tm);
3660 alpha_v=cv*rhom/rho_v;
3663 h_v=_fluides[0]->getEnthalpy(Tm,rho_v);
3664 h_l=_fluides[1]->getEnthalpy(Tm,rho_l);
3666 _VoidFraction(i)=alpha_v;
3667 _Enthalpy(i)=(m_v*h_v+m_l*h_l)/rhom;
3668 _Concentration(i)=cv;
3669 _mixtureDensity(i)=rhom;
3672 _DensiteLiquide(i)=rho_l;
3673 _DensiteVapeur(i)=rho_v;
3674 _EnthalpieLiquide(i)=h_l;
3675 _EnthalpieVapeur(i)=h_v;
3684 _VoidFraction.setTime(_time,_nbTimeStep);
3685 _Enthalpy.setTime(_time,_nbTimeStep);
3686 _Concentration.setTime(_time,_nbTimeStep);
3687 _mixtureDensity.setTime(_time,_nbTimeStep);
3688 _Pressure.setTime(_time,_nbTimeStep);
3689 _Temperature.setTime(_time,_nbTimeStep);
3690 _DensiteLiquide.setTime(_time,_nbTimeStep);
3691 _DensiteVapeur.setTime(_time,_nbTimeStep);
3692 _EnthalpieLiquide.setTime(_time,_nbTimeStep);
3693 _EnthalpieVapeur.setTime(_time,_nbTimeStep);
3694 _VitesseX.setTime(_time,_nbTimeStep);
3697 _VitesseY.setTime(_time,_nbTimeStep);
3699 _VitesseZ.setTime(_time,_nbTimeStep);
3701 if (_nbTimeStep ==0 || _restartWithNewFileName){
3705 _VoidFraction.writeVTK(allFields+"_VoidFraction");
3706 _Enthalpy.writeVTK(allFields+"_Enthalpy");
3707 _Concentration.writeVTK(allFields+"_Concentration");
3708 _mixtureDensity.writeVTK(allFields+"_Density");
3709 _Pressure.writeVTK(allFields+"_Pressure");
3710 _Temperature.writeVTK(allFields+"_Temperature");
3711 _DensiteLiquide.writeVTK(allFields+"_LiquidDensity");
3712 _DensiteVapeur.writeVTK(allFields+"_SteamDensity");
3713 _EnthalpieLiquide.writeVTK(allFields+"_LiquidEnthalpy");
3714 _EnthalpieVapeur.writeVTK(allFields+"_SteamEnthalpy");
3715 _VitesseX.writeVTK(allFields+"_VelocityX");
3718 _VitesseY.writeVTK(allFields+"_VelocityY");
3720 _VitesseZ.writeVTK(allFields+"_VelocityZ");
3724 _VoidFraction.writeMED(allFields+"_VoidFraction");
3725 _Enthalpy.writeMED(allFields+"_Enthalpy");
3726 _Concentration.writeMED(allFields+"_Concentration");
3727 _mixtureDensity.writeMED(allFields+"_Density");
3728 _Pressure.writeMED(allFields+"_Pressure");
3729 _Temperature.writeMED(allFields+"_Temperature");
3730 _DensiteLiquide.writeMED(allFields+"_LiquidDensity");
3731 _DensiteVapeur.writeMED(allFields+"_SteamDensity");
3732 _EnthalpieLiquide.writeMED(allFields+"_LiquidEnthalpy");
3733 _EnthalpieVapeur.writeMED(allFields+"_SteamEnthalpy");
3734 _VitesseX.writeMED(allFields+"_VelocityX");
3737 _VitesseY.writeMED(allFields+"_VelocityY");
3739 _VitesseZ.writeMED(allFields+"_VelocityZ");
3743 _VoidFraction.writeCSV(allFields+"_VoidFraction");
3744 _Enthalpy.writeCSV(allFields+"_Enthalpy");
3745 _Concentration.writeCSV(allFields+"_Concentration");
3746 _mixtureDensity.writeCSV(allFields+"_Density");
3747 _Pressure.writeCSV(allFields+"_Pressure");
3748 _Temperature.writeCSV(allFields+"_Temperature");
3749 _DensiteLiquide.writeCSV(allFields+"_LiquidDensity");
3750 _DensiteVapeur.writeCSV(allFields+"_SteamDensity");
3751 _EnthalpieLiquide.writeCSV(allFields+"_LiquidEnthalpy");
3752 _EnthalpieVapeur.writeCSV(allFields+"_SteamEnthalpy");
3753 _VitesseX.writeCSV(allFields+"_VelocityX");
3756 _VitesseY.writeCSV(allFields+"_VelocityY");
3758 _VitesseZ.writeCSV(allFields+"_VelocityZ");
3767 _VoidFraction.writeVTK(allFields+"_VoidFraction",false);
3768 _Enthalpy.writeVTK(allFields+"_Enthalpy",false);
3769 _Concentration.writeVTK(allFields+"_Concentration",false);
3770 _mixtureDensity.writeVTK(allFields+"_Density",false);
3771 _Pressure.writeVTK(allFields+"_Pressure",false);
3772 _Temperature.writeVTK(allFields+"_Temperature",false);
3773 _DensiteLiquide.writeVTK(allFields+"_LiquidDensity",false);
3774 _DensiteVapeur.writeVTK(allFields+"_SteamDensity",false);
3775 _EnthalpieLiquide.writeVTK(allFields+"_LiquidEnthalpy",false);
3776 _EnthalpieVapeur.writeVTK(allFields+"_SteamEnthalpy",false);
3777 _VitesseX.writeVTK(allFields+"_VelocityX",false);
3780 _VitesseY.writeVTK(allFields+"_VelocityY",false);
3782 _VitesseZ.writeVTK(allFields+"_VelocityZ",false);
3786 _VoidFraction.writeMED(allFields+"_VoidFraction",false);
3787 _Enthalpy.writeMED(allFields+"_Enthalpy",false);
3788 _Concentration.writeMED(allFields+"_Concentration",false);
3789 _mixtureDensity.writeMED(allFields+"_Density",false);
3790 _Pressure.writeMED(allFields+"_Pressure",false);
3791 _Temperature.writeMED(allFields+"_Temperature",false);
3792 _DensiteLiquide.writeMED(allFields+"_LiquidDensity",false);
3793 _DensiteVapeur.writeMED(allFields+"_SteamDensity",false);
3794 _EnthalpieLiquide.writeMED(allFields+"_LiquidEnthalpy",false);
3795 _EnthalpieVapeur.writeMED(allFields+"_SteamEnthalpy",false);
3796 _VitesseX.writeMED(allFields+"_VelocityX",false);
3799 _VitesseY.writeMED(allFields+"_VelocityY",false);
3801 _VitesseZ.writeMED(allFields+"_VelocityZ",false);
3805 _VoidFraction.writeCSV(allFields+"_VoidFraction");
3806 _Enthalpy.writeCSV(allFields+"_Enthalpy");
3807 _Concentration.writeCSV(allFields+"_Concentration");
3808 _mixtureDensity.writeCSV(allFields+"_Density");
3809 _Pressure.writeCSV(allFields+"_Pressure");
3810 _Temperature.writeCSV(allFields+"_Temperature");
3811 _DensiteLiquide.writeCSV(allFields+"_LiquidDensity");
3812 _DensiteVapeur.writeCSV(allFields+"_SteamDensity");
3813 _EnthalpieLiquide.writeCSV(allFields+"_LiquidEnthalpy");
3814 _EnthalpieVapeur.writeCSV(allFields+"_SteamEnthalpy");
3815 _VitesseX.writeCSV(allFields+"_VelocityX");
3818 _VitesseY.writeCSV(allFields+"_VelocityY");
3820 _VitesseZ.writeCSV(allFields+"_VelocityZ");
3844 if(_saveConservativeField){
3863 _Vitesse.writeVTK(prim+"_Velocity");
3866 _Vitesse.writeMED(prim+"_Velocity");
3869 _Vitesse.writeCSV(prim+"_Velocity");
3880 _VoidFraction.writeVTK(allFields+"_VoidFraction");
3881 _Enthalpy.writeVTK(allFields+"_Enthalpy");
3882 _Concentration.writeVTK(allFields+"_Concentration");
3883 _mixtureDensity.writeVTK(allFields+"_Density");
3884 _Pressure.writeVTK(allFields+"_Pressure");
3885 _Temperature.writeVTK(allFields+"_Temperature");
3886 _DensiteLiquide.writeVTK(allFields+"_LiquidDensity");
3887 _DensiteVapeur.writeVTK(allFields+"_SteamDensity");
3888 _EnthalpieLiquide.writeVTK(allFields+"_LiquidEnthalpy");
3889 _EnthalpieVapeur.writeVTK(allFields+"_SteamEnthalpy");
3890 _VitesseX.writeVTK(allFields+"_VelocityX");
3893 _VitesseY.writeVTK(allFields+"_VelocityY");
3895 _VitesseZ.writeVTK(allFields+"_VelocityZ");
3899 _VoidFraction.writeMED(allFields+"_VoidFraction");
3900 _Enthalpy.writeMED(allFields+"_Enthalpy");
3901 _Concentration.writeMED(allFields+"_Concentration");
3902 _mixtureDensity.writeMED(allFields+"_Density");
3903 _Pressure.writeMED(allFields+"_Pressure");
3904 _Temperature.writeMED(allFields+"_Temperature");
3905 _DensiteLiquide.writeMED(allFields+"_LiquidDensity");
3906 _DensiteVapeur.writeMED(allFields+"_SteamDensity");
3907 _EnthalpieLiquide.writeMED(allFields+"_LiquidEnthalpy");
3908 _EnthalpieVapeur.writeMED(allFields+"_SteamEnthalpy");
3909 _VitesseX.writeMED(allFields+"_VelocityX");
3912 _VitesseY.writeMED(allFields+"_VelocityY");
3914 _VitesseZ.writeMED(allFields+"_VelocityZ");
3918 _VoidFraction.writeCSV(allFields+"_VoidFraction");
3919 _Enthalpy.writeCSV(allFields+"_Enthalpy");
3920 _Concentration.writeCSV(allFields+"_Concentration");
3921 _mixtureDensity.writeCSV(allFields+"_Density");
3922 _Pressure.writeCSV(allFields+"_Pressure");
3923 _Temperature.writeCSV(allFields+"_Temperature");
3924 _DensiteLiquide.writeCSV(allFields+"_LiquidDensity");
3925 _DensiteVapeur.writeCSV(allFields+"_SteamDensity");
3926 _EnthalpieLiquide.writeCSV(allFields+"_LiquidEnthalpy");
3927 _EnthalpieVapeur.writeCSV(allFields+"_SteamEnthalpy");
3928 _VitesseX.writeCSV(allFields+"_VelocityX");
3931 _VitesseY.writeCSV(allFields+"_VelocityY");
3933 _VitesseZ.writeCSV(allFields+"_VelocityZ");
3940 if (_restartWithNewFileName)
3941 _restartWithNewFileName=false;
3944 void DriftModel::testConservation()
3946 double SUM, DELTA, x;
3948 for(int i=0; i<_nVar; i++)
3952 cout << "Masse totale (kg): ";
3954 cout << "Masse partielle 1 (kg): ";
3958 cout << "Energie totale (J): ";
3960 cout << "Quantite de mouvement direction "<<i-1<<" (kg.m.s^-1): ";
3966 for(int j=0; j<_Nmailles; j++)
3968 if(!_usePrimitiveVarsInNewton)
3969 VecGetValues(_old, 1, &I, &x);//on recupere la valeur du champ
3971 VecGetValues(_primitiveVars, 1, &I, &x);//on recupere la valeur du champ
3972 SUM += x*_mesh.getCell(j).getMeasure();
3973 VecGetValues(_newtonVariation, 1, &I, &x);//on recupere la variation du champ
3974 DELTA += x*_mesh.getCell(j).getMeasure();
3977 if(fabs(SUM)>_precision)
3978 cout << SUM << ", variation relative: " << fabs(DELTA /SUM) << endl;
3980 cout << " a une somme nulle, variation absolue: " << fabs(DELTA) << endl;