Subversion Repositories sifel

Compare Revisions

Ignore whitespace Rev 776 → Rev 777

/trunk/SIFEL/METR/SRC/consol_awf2c.cpp
1249,7 → 1249,7
capwu = sw*alpha;//p. 95
 
n = get_porosity(ipp);
dsr_depsv = 0.0;//dsr_depsv = Tm->givenontransq(der_saturation_deg_depsv, ipp); //actual derivative of saturation degree with respect to volumetric strain
dsr_depsv = Tm->givenontransq(der_saturation_deg_depsv, ipp); //actual derivative of saturation degree with respect to volumetric strain
capwu = capwu + n*dsr_depsv;//influence of volumetric strain on saturation degree
break;
}
1322,7 → 1322,6
}
}
 
//capgu = 0.0;//debug!!??
return(capgu);
 
}
/trunk/SIFEL/METR/SRC/consol_hawf3c.cpp
1173,7 → 1173,6
}
 
 
 
/**
function returns conductivity coefficient of the general material
@param pw - water pressure
1261,9 → 1260,11
break;
}
case lewis_and_schrefler3_mefel_coup:{//Lewis and Schrefler's book
n = get_porosity(pw,pg,t,ipp);
dsr_depsv = Tm->givenontransq(der_saturation_deg_depsv, ipp); //actual derivative of saturation degree with respect to volumetric strain
capwu = capwu + n*(rhow - rhogw)*dsr_depsv;//influence of volumetric strain on saturation degree
//n = get_porosity(pw,pg,t,ipp);
//dsr_depsv = Tm->givenontransq(der_saturation_deg_depsv, ipp); //actual derivative of saturation degree with respect to volumetric strain
//capwu = capwu + n*(rhow - rhogw)*dsr_depsv;//influence of volumetric strain on saturation degree
 
capwu = 0.0;// this part is in TRFEL already included
break;
}
default:{
1325,9 → 1326,11
break;
}
case lewis_and_schrefler3_mefel_coup:{//Lewis and Schrefler's book
n = get_porosity(pw,pg,t,ipp);
dsr_depsv = Tm->givenontransq(der_saturation_deg_depsv, ipp); //actual derivative of saturation degree with respect to volumetric strain
capgu = capgu - n*rhoga*dsr_depsv;//influence of volumetric strain on saturation degree
//n = get_porosity(pw,pg,t,ipp);
//dsr_depsv = Tm->givenontransq(der_saturation_deg_depsv, ipp); //actual derivative of saturation degree with respect to volumetric strain
//capgu = capgu - n*rhoga*dsr_depsv;//influence of volumetric strain on saturation degree
 
capgu = 0.0;// this part is in TRFEL already included
break;
}
default:{
1335,6 → 1338,7
fprintf (stderr,"\n in function (%s, line %d).\n",__FILE__,__LINE__);
}
}
 
return(capgu);
}
 
1361,6 → 1365,7
case lewis_and_schrefler3_mefel_coup:{//Lewis and Schrefler's book
betas = get_betas(ipp);
kut = -1.0*betas/3.0;//p. 348
 
kut = 0.0;// this part is in MEFEL already included
break;
}
1436,9 → 1441,11
break;
}
case lewis_and_schrefler3_mefel_coup:{//Lewis and Schrefler's book
n = get_porosity(pw,pg,t,ipp);
dsr_depsv = Tm->givenontransq(der_saturation_deg_depsv, ipp); //actual derivative of saturation degree with respect to volumetric strain
captu = captu - dhvap*rhow*n*dsr_depsv;//influence of volumetric strain on saturation degree
//n = get_porosity(pw,pg,t,ipp);
//dsr_depsv = Tm->givenontransq(der_saturation_deg_depsv, ipp); //actual derivative of saturation degree with respect to volumetric strain
//captu = captu - dhvap*rhow*n*dsr_depsv;//influence of volumetric strain on saturation degree
 
captu = 0.0;// this part is in TRFEL already included
break;
}
default:{
1447,6 → 1454,7
}
}
 
captu = 0.0;//debug??!!
return(captu);
}
 
2236,7 → 2244,7
// Gaw-Maj-Sch "Numerical analysis of hygro thermal behaviour and damage of concrete at high T"
// dhvap = 2.672e+5 * pow((t-tcr),0.38) (eq. 49)
 
dhvap = dhvap*0.0;//debug??!!
//dhvap = dhvap*0.0;//debug??!!
 
return(dhvap);
}
/trunk/SIFEL/METR/SRC/globmatc.cpp
679,7 → 679,7
Gtm->give_code_numbers (i,rcn.a);
//elemvalues (i, r);
prescvalues (r.a,rcn.a,tndofe);
prescvalues (r.a,rcn.a,mndofe); //old wrong: prescvalues (r.a,rcn.a,tndofe);
 
mxv (km,r,f);
cmulv (-1.0,f);
688,7 → 688,7
 
// contribution from time derivative of prescribed values on boundaries
// prescribed values from the previous time step
prevprescvalues (pr.a,rcn.a,tndofe);
prevprescvalues (pr.a,rcn.a,mndofe);//old: prevprescvalues (pr.a,rcn.a,tndofe);
// the length of time step
dt=Tp->timecont.forwarddt;
// computation of the time derivative of
718,10 → 718,10
 
lower_cond_coupl_mat (i,lcid,km);
Gtm->give_code_numbers (i,rcn.a);
Gtt->give_code_numbers (i,rcn.a);//old:Gtm->give_code_numbers (i,rcn.a);
//elemvalues (i, r);
prescvalues (r.a,rcn.a,mndofe);
prescvalues (r.a,rcn.a,tndofe);//old: prescvalues (r.a,rcn.a,mndofe);
 
mxv (km,r,f);
cmulv (-1.0,f);
730,7 → 730,7
 
// contribution from time derivative of prescribed values on boundaries
// prescribed values from the previous time step
prevprescvalues (pr.a,rcn.a,mndofe);
prevprescvalues (pr.a,rcn.a,tndofe);//old: prevprescvalues (pr.a,rcn.a,mndofe);
// the length of time step
dt=Tp->timecont.forwarddt;
// computation of the time derivative of
/trunk/SIFEL/TRFEL/SRC/C60bazant.cpp
10,9 → 10,11
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "aliast.h"
#include "constrel.h"
#include "C60bazant.h"
#include "globalt.h"
#include "globmatt.h"
 
C60bazmat::C60bazmat()
{
146,6 → 148,10
fprintf (Outt,"s = %lf\n",sw);
fprintf (Outt,"mt = %lf\n",mt);
*/
 
//fprintf(Outt,"sw = %e, pc = %e, t = %e, ipp = %ld\n",sw,pc,t,0);
//fflush(Outt);
 
return(sw);
}
 
207,6 → 213,9
dsw_dpc = dsdc;
 
//fprintf(Outt,"dsw_dpc = %e, pc = %e, t = %e, ipp = %ld\n",dsw_dpc,pc,t,0);
//fflush(Outt);
 
return(dsw_dpc);
}
 
287,6 → 296,9
dsw_dt = dsdt;
 
//fprintf(Outt,"dsw_dt = %e, pc = %e, t = %e, ipp = %ld\n",dsw_dt,pc,t,0);
//fflush(Outt);
 
return(dsw_dt);
}
 
363,7 → 375,7
// set a lower limit of krg (e.g. 0.0001)to avoid numerical problems (pag.10)
if (krg<0.0001){
krg=0.0001;
print_err("relative permeability was modified",__FILE__,__LINE__,__func__);
//print_err("relative permeability was modified",__FILE__,__LINE__,__func__);
}
return(krg);
412,8 → 424,8
 
tau = C60baz_tau(pc,t);
dd = tau*mw*ma/gasr;//according to Fortran code
//dd = tau;//according to Francesco's PhD thesis?
//dd = tau*mw*ma/gasr;//according to Fortran code
dd = tau;//according to multiphase.cpp and Francesco's PhD thesis
 
return(dd);
}
803,3 → 815,27
*/
void C60bazmat::print(FILE *out)
{}
 
 
/**
The function returns ordered dof names of primary unknowns
required by the model.
 
@param dofname - array of uknown name for particular nodal dofs (output)
dofname[i] = name of i-th nodal unknown (for names see aliast.h - enum namevart)
@param ntm - number of transported media = number of nodal dof = length of array dofname
 
Created by Tomas Krejci according to Tomas Koudelka, 09/12/2022
*/
void C60bazmat::give_dof_names(namevart *dofname, long ntm)
{
if (ntm < 1)
{
print_err("the model defines %ld unknowns but number of transported media is %ld",
__FILE__, __LINE__, __func__, 1, ntm);
abort();
}
dofname[0] = trf_capillary_press;
dofname[1] = trf_gas_press;
dofname[2] = trf_temperature;
}
/trunk/SIFEL/TRFEL/SRC/C60bazant.h
43,6 → 43,8
void read(XFILE *in);
void print(FILE *out);
 
void give_dof_names(namevart *dofname, long ntm);
 
private:
double mw;
double ma;
/trunk/SIFEL/TRFEL/SRC/Makefile
31,7 → 31,7
####################################################################
# List of source files
#
SRCS = aepointst.cpp backupsolt.cpp baroghelB.cpp bazped.cpp bnodvalt.cpp boundfluxes.cpp carbmat1.cpp cemhydmat.cpp
SRCS = aepointst.cpp backupsolt.cpp baroghelB.cpp baroghel_reten.cpp bazant_reten.cpp bazped.cpp bnodvalt.cpp boundfluxes.cpp carbmat1.cpp cemhydmat.cpp
SRCS += cerny_concrete.cpp climatcond.cpp climatcond2.cpp concreteB.cpp consol_awf1.cpp consol_wf1.cpp consol_awf2.cpp consol_hawf3.cpp consol_hwf2.cpp consol_wf2.cpp constrel.cpp cpnnpsolvert.cpp
SRCS += cpnpsolvert.cpp crsection1d.cpp crsection2d.cpp crsection3d.cpp C30baroghel.cpp
SRCS += C60baroghel.cpp C30bazant.cpp C60bazant.cpp damisotrmat.cpp dampermeability.cpp devriesmat.cpp discisotrmat.cpp discmat.cpp
/trunk/SIFEL/TRFEL/SRC/aliast.h
34,10 → 34,10
 
// aliases for name/type of primary variables/unknowns (MUST be numbered from 1 and one by one)
enum namevart {trf_temperature=1,trf_rel_humidity=2,trf_water_press=3,trf_gas_press=4, trf_pore_press=5,
trf_salt_conc=6,trf_hydraulic_head=7,trf_vapor_content=8,trf_moisture=9, trf_press_water_vapor=10};
trf_salt_conc=6,trf_hydraulic_head=7,trf_vapor_content=8,trf_moisture=9, trf_press_water_vapor=10,trf_capillary_press=3};
const enumstr namevartstr[] = {{"trf_temperature",1}, {"trf_rel_humidity",2}, {"trf_water_press",3},
{"trf_gas_press",4}, {"trf_pore_press",5}, {"trf_salt_conc",6},
{"trf_hydraulic_head",7},{"trf_vapor_content",8},{"trf_moisture",9}, {"trf_press_water_vapor",10}};
{"trf_hydraulic_head",7},{"trf_vapor_content",8},{"trf_moisture",9}, {"trf_press_water_vapor",10}, {"trf_press_capillary_vapor",11}};
const kwdset namevart_kwdset(sizeof(namevartstr)/sizeof(*namevartstr),namevartstr);
 
 
170,8 → 170,8
const kwdset airwaterflowtype_kwdset(sizeof(airwaterflowtypestr)/sizeof(*airwaterflowtypestr),airwaterflowtypestr);
 
// aliases for models of heat, air and water (three phase) flow in deforming porous medium (soils) consol_hawf3 = consolidation
enum heatairwaterflowtype {lewis_and_schrefler3=1,lewis_and_schrefler3_2=2,van_genuchten3=4,lewis_and_schrefler3_mefel=7};
const enumstr heatairwaterflowtypestr[] = {{"lewis_and_schrefler3",1},{"lewis_and_schrefler3_2",2}, {"van_genuchten3",4},{"lewis_and_schrefler3_mefel",7}};
enum heatairwaterflowtype {lewis_and_schrefler3=1,lewis_and_schrefler3_2=2,artificial3=3,van_genuchten3=4,lewis_and_schrefler3_mefel=7};
const enumstr heatairwaterflowtypestr[] = {{"lewis_and_schrefler3",1},{"lewis_and_schrefler3_2",2},{"artificial3",3}, {"van_genuchten3",4},{"lewis_and_schrefler3_mefel",7}};
const kwdset heatairwaterflowtype_kwdset(sizeof(heatairwaterflowtypestr)/sizeof(*heatairwaterflowtypestr),heatairwaterflowtypestr);
 
 
188,9 → 188,9
enum sourceloc {nod=1,elem=2};
 
// aliases for models of retention curve type
enum sr_type {lewis_and_schrefler_sr=1,gardner_exponential_sr=2,potts_log_linear_sr=3,van_genuchten_sr=4,van_genuchten2_sr=5,mefel_sr=7,table_sr=8};
enum sr_type {lewis_and_schrefler_sr=1,gardner_exponential_sr=2,potts_log_linear_sr=3,van_genuchten_sr=4,van_genuchten2_sr=5,mefel_sr=7,table_sr=8,bazant_sr=9,baroghel_sr=10};
const enumstr sr_typestr[] = {{"lewis_and_schrefler_sr",1}, {"gardner_exponential_sr",2}, {"potts_log_linear_sr",3}, {"van_genuchten_sr",4},
{"mefel_sr",7},{"table_sr",8}};
{"mefel_sr",7},{"table_sr",8},{"bazant_sr",9},{"baroghel_sr",10}};
const kwdset sr_type_kwdset(sizeof(sr_typestr)/sizeof(*sr_typestr),sr_typestr);
 
 
/trunk/SIFEL/TRFEL/SRC/baroghelB.cpp
8,9 → 8,11
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "aliast.h"
#include "constrel.h"
#include "baroghelB.h"
#include "globalt.h"
#include "globmatt.h"
 
baroghelmat::baroghelmat()
{
566,3 → 568,29
{
fprintf (out," %e %e %e %e %e %e %e %e %e", ab, bb, c, lambdab, cpb, rhosb, phib, emod, nu);
}
 
 
 
 
/**
The function returns ordered dof names of primary unknowns
required by the model.
 
@param dofname - array of uknown name for particular nodal dofs (output)
dofname[i] = name of i-th nodal unknown (for names see aliast.h - enum namevart)
@param ntm - number of transported media = number of nodal dof = length of array dofname
 
Created by Tomas Krejci according to Tomas Koudelka, 09/12/2022
*/
void baroghelmat::give_dof_names(namevart *dofname, long ntm)
{
if (ntm < 1)
{
print_err("the model defines %ld unknowns but number of transported media is %ld",
__FILE__, __LINE__, __func__, 1, ntm);
abort();
}
dofname[0] = trf_capillary_press;
dofname[1] = trf_gas_press;
dofname[2] = trf_temperature;
}
/trunk/SIFEL/TRFEL/SRC/baroghelB.h
39,6 → 39,8
void read(XFILE *in);
void print(FILE *out);
 
void give_dof_names(namevart *dofname, long ntm);
 
private:
double gasr;
/trunk/SIFEL/TRFEL/SRC/baroghel_reten.cpp
0,0 → 1,204
/*
File: baroghel_reten.cpp
Author: Tomas Krejci, 09/12/2022
Purpose: Calculates the retention curve assuming Baroghel formulation extended for high temperature
sources: C60sifel.f90 and SATBGMsifel.f90 from Padova
*/
 
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "aliast.h"
#include "baroghel_reten.h"
#include "globalt.h"
 
baroghel_reten::baroghel_reten()
{
t0 = 273.15; //reference temperature
tcr = 647.3; //critical point of water
 
//for concrete BO:
b = 2.2748; //retention curve coefficient b - exponent
q3 = 18.6237e6; //retention curve coefficient q3[Pa]
q2 = 7.0e6; //retention curve coefficient q2[Pa]
n = 1.2; //retention curve coefficient N - exponent
z = 0.5; //z is governing the transition through critical temperature
 
//for concrete BH:
//b = 2.0601; //retention curve coefficient b - exponent
//q3 = 46.9364e6; //retention curve coefficient q3[Pa]
//q2 = 7.0e6; //retention curve coefficient q2[Pa]
//n = 1.2; //retention curve coefficient N - exponent
//z = 0.5; //z is governing the transition through critical temperature
 
//for concrete CO:
//b = 2.1684; //retention curve coefficient b - exponent
//q3 = 37.5479e6; //retention curve coefficient q3[Pa]
//q2 = 7.0e6; //retention curve coefficient q2[Pa]
//n = 1.2; //retention curve coefficient N - exponent
//z = 0.5; //z is governing the transition through critical temperature
 
//for concrete CH:
//b = 1.9570; //retention curve coefficient b - exponent
//q3 = 96.2837e6; //retention curve coefficient q3[Pa]
//q2 = 7.0e6; //retention curve coefficient q2[Pa]
//n = 1.2; //retention curve coefficient N - exponent
//z = 0.5; //z is governing the transition through critical temperature
}
 
baroghel_reten::~baroghel_reten()
{
}
 
 
/**
function reads parameters
@param in - input file
*/
void baroghel_reten::read(XFILE *in)
{
xfscanf (in,"%lf %lf %lf %lf %lf", &b, &q2, &q3, &n, &z);
}
 
 
/**
function prints parameters
@param out - output file
*/
void baroghel_reten::print(FILE *out)
{
fprintf (out," %e %e %e %e %e\n", b, q2, q3, n, z);
}
 
 
/**
function computes degree of saturation (sorption curve)
@param pc - capillary pressure
@param t - temperature
 
@retval sw - degree of saturation
*/
double baroghel_reten::sw(double pc,double t)
{
double sw,n,a,b,tt,q0,q2,q3,z,e,e0,g;
b = 2.27;
q3 = 18.62e6;
q2 = 7.0e6;
n = 1.2;
z = 0.5; //z is governing the transition through critical temperature
if (t <= 373.15)
a = q3;
else{
tt=(t-373.15)/(647.15-373.15);
q0=(q3-q2)*((2.*(tt*tt*tt))-(3.0*(tt*tt))+1.0);
q2 = 25.0e6;
a = q0 + q2;
}
 
if (t <= tcr)
e = pow(((tcr - t0)/(tcr - t)),n);
else{
e0=pow(((tcr-293.15)/z),n);
e = n/z*e0*t + e0 - n/z*e0*(tcr - z);
}
 
g = pow((e/a*pc),(b/(b-1.0)));
sw = pow((g + 1.0),(-1.0/b));
return(sw);
}
 
/**
function computes partial derivative of degree of saturation with respect to pc
@param pc - capillary pressure
@param t - temperature
 
@retval ds_dpc - partial derivative of degree of saturation with respect to pc
*/
double baroghel_reten::dsw_dpc(double pc,double t)
{
double dsw_dpc;
double n,a,b,dg_dpc,tt,q0,q2,q3,z,e,e0,g;
b = 2.27;
q3 = 18.62e6 ;
q2 = 7.0e6;
n = 1.2;
z = 0.5; //z is governing the transition through critical temperature
if (t <= 373.15)
a = q3;
else{
tt=(t-373.15)/(647.15-373.15);
q0=(q3-q2)*((2.*(tt*tt*tt))-(3.0*(tt*tt))+1.0);
q2 = 25.0e6;
a = q0 + q2;
}
if (t <= tcr){
e = pow(((tcr - t0)/(tcr - t)),n);
}
else{
e0=pow(((tcr-293.15)/z),n);
e = n/z*e0*t + e0 - n/z*e0*(tcr - z);
}
g = pow((e/a*pc),(b/(b-1.0)));
dg_dpc = (b/(b-1.0))*e/a*pow((e/a*pc),(b/(b-1.0)-1.0));
dsw_dpc = (-1.0/b)*pow((g + 1.0),(-1.0-1.0/b))*dg_dpc;
 
return(dsw_dpc);
}
 
/**
function computes partial derivative of degree of saturation with respect to t
@param pc - capillary pressure
@param t - temperature
 
@retval dsw_dt - partial derivative of degree of saturation with respect to t
*/
double baroghel_reten::dsw_dt(double pc,double t)
{
double dsw_dt,e,b,n,de_dt,g,a,da_dt,dea_dt,dg_dt,tt,dtt_dt,q0,dq0_dt,q2,q3,e0,z;
b = 2.27;
q3 = 18.62e6 ;
q2 = 7.0e6;
n = 1.2;
z = 0.5; //z is governing the transition through critical temperature
if (t <= 373.15){
a = q3;
da_dt = 0.0;
}
else{
tt=(t-373.15)/(647.15-373.15);
q0=(q3-q2)*((2.*(tt*tt*tt))-(3.0*(tt*tt))+1.0);
dtt_dt = 1.0/(647.15-373.15);
dq0_dt = (q3-q2)*(6.0*tt*tt*dtt_dt - 6.0*tt*dtt_dt);
q2 = 25.0e6;
a = q0 + q2;
da_dt = dq0_dt;
}
if (t <= tcr){
e = pow(((tcr - t0)/(tcr - t)),n);
de_dt = n*pow(((tcr - t0)/(tcr - t)),(n-1.0))*(tcr - t0)/(tcr - t)/(tcr - t);
}
else{
e0=pow(((tcr-293.15)/z),n);
e = n/z*e0*t + e0 - n/z*e0*(tcr - z);
de_dt = n/z*e0;
}
 
g = pow((e/a*pc),(b/(b-1.0)));
dea_dt = (de_dt*a - da_dt*e)/a/a;
dg_dt = (b/(b-1.0))*pc*pow((e/a*pc),(b/(b-1.0)-1.0))*dea_dt;
dsw_dt = (-1.0/b)*pow((g + 1.0),(-1.0-1.0/b))*dg_dt;
 
return(dsw_dt);
}
/trunk/SIFEL/TRFEL/SRC/baroghel_reten.h
0,0 → 1,23
#ifndef BAROGHEL_RETEN_H
#define BAROGHEL_RETEN_H
 
#include "genfile.h"
 
class baroghel_reten
{
public:
baroghel_reten(); //constructor
~baroghel_reten(); //destructor
 
void read(XFILE *in);
void print(FILE *out);
double sw(double pc,double t);
double dsw_dpc(double pc,double t);
double dsw_dt(double pc,double t);
 
private:
double t0,tcr,b,q2,q3,n,z;
};
 
#endif
/trunk/SIFEL/TRFEL/SRC/bazant_reten.cpp
0,0 → 1,233
/*
File: bazant_reten.cpp
Author: Tomas Krejci, 09/12/2022
Purpose: Calculates the retention curve assuming Bazant formulation extended for high temperature
sources: C60sifel.f90 and SATBGMsifel.f90 from Padova
*/
 
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "aliast.h"
#include "constrel.h"
#include "bazant_reten.h"
#include "globalt.h"
 
bazant_reten::bazant_reten()
{
t0 = 273.15; //reference temperature
p0 = 101325.0; //pressure
tcr = 647.3; //critical point of water
}
bazant_reten::~bazant_reten()
{}
 
/**
function computes degree of saturation (desorption curve)
 
@param pc - capillary pressure
@param t - temperature
 
@retval sw - degree of saturation
*/
double bazant_reten::sat(double pc,double t)
{
double sw,rh,tamb,tem,tt1,mt1,fcm,mt;
state_eq tt;
 
tamb= 25.0;
rh = tt.get_rh(pc,t);
//temperature in Celsius degrees
if (t < tcr){
tem = t - 273.15;
if (tem < tamb)
tem = tamb;
}
else
tem = tcr - 273.15;
 
//function SATUR_A.f90
tt1 = ((tem+10.0)/(25.0+10.0))*((tem+10.0)/(25.0+10.0));//modified temperature from Bazant formulation
mt1 = 1.04 - (tt1/(24.0+tt1));//function m from Bazant formulation
if(t < tcr){
fcm = (tcr-t)/(tcr-298.15);
mt = mt1*fcm;//bazant parameter (below critical point)
if(mt < 1.0e-4)
mt = 1.0e-4;
}
else
mt = 1.0e-4; //Bazant parameter (above critical point)
//conversion pow(a,b) -> exp(b*log(a))
sw = pow(rh,(1.0/mt));//saturation calculation
//sw = exp(log(rh)/mt);//saturation calculation
if(sw < 1.0e-3)//limit for saturation level
{
sw=1.0e-3;
fprintf (Outt,"\n\n Uprava saturace, function sat (C60bazant.cpp)");
}
 
/*
fprintf (Outt,"\n\n");
fprintf (Outt,"pc = %lf\n",pc);
fprintf (Outt,"t = %lf\n",t);
fprintf (Outt,"rh = %lf\n",rh);
fprintf (Outt,"s = %lf\n",sw);
fprintf (Outt,"mt = %lf\n",mt);
*/
//fprintf(Outt,"sw = %e, pc = %e, t = %e, ipp = %ld\n",sw,pc,t,0);
//fflush(Outt);
 
return(sw);
}
 
/**
function computes partial derivative of degree of saturation with respect to pc
@param pc - capillary pressure
@param t - temperature
 
@retval ds_dpc - partial derivative of degree of saturation with respect to pc
*/
double bazant_reten::dsat_dpc(double pc,double t)
{
double dsw_dpc;
double sw,rh,tamb,tem,tt1,mt1,fcm,mt;
double ssmh,drhdc,dsdc;
state_eq tt;
 
tamb= 25.0;
rh = tt.get_rh(pc,t);
drhdc = tt.get_drh_dpc(pc,t);
 
//temperature in Celsius degrees
if (t < tcr){
tem = t - 273.15;
if (tem < tamb)
tem = tamb;
}
else
tem = tcr - 273.15;
 
//function SATUR_A.f90
tt1 = ((tem+10.0)/(25.0+10.0))*((tem+10.0)/(25.0+10.0));//modified temperature from Bazant formulation
mt1 = 1.04 - (tt1/(24.0+tt1));//function m from Bazant formulation
if(t < tcr){
fcm = (tcr-t)/(tcr-298.15);
mt = mt1*fcm;//bazant parameter (below critical point)
if(mt<1.0e-4) mt=1.0e-4;
}
else{
mt = 1.0e-4; //Bazant parameter (above critical point)
}
//conversion pow(a,b) -> exp(b*log(a))
sw = pow(rh,(1.0/mt));//saturation calculation
//sw = exp(1.0/mt*log(rh));//saturation calculation
//sw = sat(pc,t);
if(sw < 1.0e-3)//limit for saturation level
{
sw=1.0e-3;
fprintf (Outt,"\n\n Uprava saturace, function dsat_dpc (C60bazant.cpp)");
}
 
ssmh = sw/mt/rh;
dsdc = ssmh*drhdc;
//if(sw <= 1.0e-10)//limit for saturation level
//dsdc=0.0;
dsw_dpc = dsdc;
 
//fprintf(Outt,"dsw_dpc = %e, pc = %e, t = %e, ipp = %ld\n",dsw_dpc,pc,t,0);
//fflush(Outt);
 
return(dsw_dpc);
}
 
/**
function computes partial derivative of degree of saturation with respect to t
@param pc - capillary pressure
@param t - temperature
 
@retval dsw_dt - partial derivative of degree of saturation with respect to t
*/
double bazant_reten::dsat_dt(double pc,double t)
{
double dsw_dt;
double sw,rh,tamb,tem,tt1,mt1,fcm,mt,dtt1dt,dm1tdt,dfcmdt,dmtdt;
double lnh,xfun,drhdt,ssmh,dsdt;
state_eq tt;
 
tamb= 25.0;
rh = tt.get_rh(pc,t);
drhdt = tt.get_drh_dt(pc,t);
 
//temperature in Celsius degrees
if (t < tcr){
tem = t - 273.15;
if (tem < tamb)
tem = tamb;
}
else
tem = tcr - 273.15;
 
//function SATUR_A.f90
tt1 = ((tem+10.)/(25.0+10.))*((tem+10.)/(25.0+10.));//modified temperature from Bazant formulation
mt1 = 1.04 - (tt1/(24.+tt1));//function m from Bazant formulation
if(t < tcr){
fcm = (tcr-t)/(tcr-298.15);
mt = mt1*fcm;//bazant parameter (below critical point)
dtt1dt = 2.0*(tem+10.0)/((25.0+10.0)*(25.0+10.0));
dm1tdt = -24.0*dtt1dt/((24.0+tt1)*(24.0+tt1));
dfcmdt = -1.0/(tcr-298.15);
if(mt<1.0e-4) mt=1.0e-4;
}
else{
fcm = 0.0;
mt = 1.0e-4; //Bazant parameter (above critical point)
dtt1dt = 0.0;
dm1tdt = 0.0;
dfcmdt = 0.0;
}
dmtdt = dm1tdt*fcm + mt1*dfcmdt;
//conversion pow(a,b) -> exp(b*log(a))
sw = pow(rh,(1.0/mt));//saturation calculation
//sw = exp(1.0/mt*log(rh));//saturation calculation
lnh = log(rh);
if(sw < 1.0e-3)//limit for saturation level
{
sw=1.0e-3;
fprintf (Outt,"\n\n Uprava saturace, function dsat_dt (C60bazant.cpp)");
}
if(t < tcr){//t=temperature on gauss point ; tcr is critical temperature
xfun = drhdt - rh/mt*lnh*dmtdt;
ssmh = sw/mt/rh;
dsdt = ssmh*xfun;
}
else{
xfun = drhdt - rh/mt*lnh*dmtdt ;
ssmh = sw/mt/rh;
dsdt = 0.0;
}
//if(sw <= 1.0e-10)//limit for saturation level
//dsdt=0.0;//debug??!!
dsw_dt = dsdt;
 
//fprintf(Outt,"dsw_dt = %e, pc = %e, t = %e, ipp = %ld\n",dsw_dt,pc,t,0);
//fflush(Outt);
 
return(dsw_dt);
}
/trunk/SIFEL/TRFEL/SRC/bazant_reten.h
0,0 → 1,22
#ifndef BAZANT_RETEN_H
#define BAZANT_RETEN_H
 
#include "genfile.h"
 
class bazant_reten
{
public:
bazant_reten(); //constructor
~bazant_reten(); //destructor
 
double sat(double pc,double t);
double dsat_dpc(double pc,double t);
double dsat_dt(double pc,double t);
 
private:
double t0;
double p0;
double tcr;
};
 
#endif
/trunk/SIFEL/TRFEL/SRC/consol_awf1.cpp
786,10 → 786,7
*/
double con_awf1mat::get_kww(double pw, long ipp)
{
double s,ds;
double krw,kww;
double e,eps_v,sigma_m_eff,e0_bar;
double a,b,c,e0_hat,lam_hat;
switch (model_type){
1360,7 → 1357,6
*/
void con_awf1mat::waterpress_check(double &pw,long ipp)
{
 
}
 
/**
/trunk/SIFEL/TRFEL/SRC/consol_awf2.cpp
99,7 → 99,6
switch (model_type){
case lewis_and_schrefler2:{//Lewis and Schrefler's model approach
lewis_ret.read(in);
break;
}
case lewis_and_schrefler2_mefel:{//Lewis and Schrefler's approach coupled with mechanics
1761,18 → 1760,27
*/
void con_awf2mat::gaspress_check(double pw,double &pg,long ipp)
{
switch (model_type){
case lewis_and_schrefler2://Lewis and Schrefler's book
case van_genuchten2:
case lewis_and_schrefler2_mefel:{//Lewis and Schrefler's book, saturation degree ind its derivative are obtained from mefel;
break;
double pgw;
//general
//gas pressure check
pgw = get_pgw(pw,pg,t0);
if (pgw > pg){
pg = pgw;
fprintf(Outt,"gas pressure was modified in integration point No. %ld",ipp);
}
default:{
print_err("unknown model type is required", __FILE__, __LINE__, __func__);
abort();
 
//general
//fully saturated state
if(pgw < 100.0){
pg = 101325.0;
}
 
// M. Starnoni 24-11-2010
if (pg < 0.0){
print_err("gas pressure was modified in integration point No. %ld",__FILE__,__LINE__,__func__,ipp);
pg = 0.0;
}
}
 
 
1786,20 → 1794,6
*/
void con_awf2mat::waterpress_check(double &pw,double pg,long ipp)
{
switch (model_type){
case lewis_and_schrefler2:{//Lewis and Schrefler's book
break;
}
case van_genuchten2:
case lewis_and_schrefler2_mefel:{//Lewis and Schrefler's book, saturation degree ind its derivative are obtained from mefel;
break;
}
default:{
print_err("unknown model type is required", __FILE__, __LINE__, __func__);
abort();
}
}
 
}
 
/**
/trunk/SIFEL/TRFEL/SRC/consol_hawf3.cpp
23,15 → 23,16
#include <stdlib.h>
#include <math.h>
#include "aliast.h"
#include "constrel.h"
#include "consol_hawf3.h"
#include "globalt.h"
#include "C60bazant.h" //testing
#include "C60baroghel.h" //testing
#include "baroghel_reten.h"
#include "bazant_reten.h"
#include "globmatt.h"
 
con_hawf3mat::con_hawf3mat()
{
compress = 0; //compressible grains: 0=no; 1=yes
compress = 0; //compressible grains: 0=no; 1=yes
kintr_type = 0; //intrinsic permability calculation type
krw_type = 0; //relative permeability calculation type
krg_type = 0; //relative permability calculation type
137,6 → 138,10
 
pw_bc = 0.0; //free boundary pressure
rel_gas_press = 0;
 
//parameters for the artificial material:
kww0 = kwg0 = kwt0 = kgw0 = kgg0 = kgt0 = ktw0 = ktg0 = ktt0 = 0.0;
capww0 = capwg0 = capwt0 = capgw0 = capgg0 = capgt0 = captw0 = captg0 = captt0 =0.0;
}
con_hawf3mat::~con_hawf3mat()
{}
158,8 → 163,11
xfscanf (in,"%le %le %le %le %le %le %le %d %d %d %d %d %d %d %d", &alpha0, &ks0, &phi0, &rhos0, &pw_bc, &tau0, &kintr0, &kintr_type, &krw_type, &krg_type, &deff_type, &sr_type, &lambda_type, &cps_type, &betas_type);
 
switch (model_type){
case artificial3:{//artificial isotropic material for non-isotherma air-water flow
xfscanf (in,"%le %le %le %le %le %le %le %le %le", &kww0, &kwg0, &kwt0, &kgw0, &kgg0, &kgt0, &ktw0, &ktg0, &ktt0);
xfscanf (in,"%le %le %le %le %le %le %le %le %le", &capww0, &capwg0, &capwt0, &capgw0, &capgg0, &capgt0, &captw0, &captg0, &captt0);
}
case lewis_and_schrefler3:{//Lewis and Schrefler's model approach
lewis_ret.read(in);
break;
}
case lewis_and_schrefler3_mefel:{//Lewis and Schrefler's model approach coupled with mechanics
172,141 → 180,152
}
}
 
//intrinsic permebility calculation:
switch (kintr_type){
case 0:{//constant
break;
if(model_type != artificial3){
//intrinsic permebility calculation:
switch (kintr_type){
case 0:{//constant
break;
}
case 1:{//dependent on porosity
xfscanf (in,"%le %le", &bb1, &phi01);
break;
}
default:{
print_err ("unknown model type is required", __FILE__, __LINE__, __func__);
abort();
}
}
//relative permebility calculation:
switch (krw_type){
case 0:{//constant
break;
}
case 1:
case 2:{//dependent on saturation degree:
xfscanf (in,"%le %le", &sirr, &ssat);
break;
}
case 3:{//exponential
xfscanf (in,"%le", &lambda_krw);
break;
}
case 9:{//bazant
break;
}
default:{
print_err ("unknown model type is required", __FILE__, __LINE__, __func__);
abort();
}
}
//gas relative permebility calculation:
switch (krg_type){
case 0:{//constant
break;
}
case 1:
case 2:{//dependent on saturation degree:
//xfscanf (in,"%le %le", &sirr, &ssat); //not finished
break;
}
case 3:{//exponential
//xfscanf (in,"%le", &lambda_krw); //not finished
break;
}
default:{
print_err ("unknown model type is required", __FILE__, __LINE__, __func__);
abort();
}
}
//retention curve type:
switch (sr_type){
case baroghel_sr:{//Baroghle-Bouny approach
baroghel_ret.read(in);
break;
}
case bazant_sr:{//Bazant approach
break;
}
case lewis_and_schrefler_sr:{//Lewis and Schrefler's book
lewis_ret.read(in);
break;
}
case van_genuchten_sr:{//partially saturated medium =Van Genuchten model
van_genuchten_ret.read(in);
break;
}
case van_genuchten2_sr:{//partially saturated medium =Van Genuchten model
//van_genuchten_ret.read2(in);
break;
}
case mefel_sr:{//from MEFEL
break;
}
case table_sr:{//from table
//reading of retention curve:
data.read (in);
break;
}
default:{
print_err ("unknown model type is required", __FILE__, __LINE__, __func__);
abort();
}
}
//conductivity calculation:
switch (lambda_type){
case 0:{//constant
xfscanf (in,"%le", &lambda0);
break;
}
case 1:{//dependent on moisture
xfscanf (in,"%le %le", &lambda_dry, &lambda_wet);
break;
}
case 2:{//dependent on moisture
xfscanf (in,"%le %le %le %le", &lambda_dry, &lambda_wet, &sr_dry, &sr_wet);
break;
}
default:{
print_err ("unknown model type is required", __FILE__, __LINE__, __func__);
abort();
}
}
//specific heat calculation:
switch (cps_type){
case 0:{//constant
xfscanf (in,"%le", &cps0);
break;
}
case 1:{//dependent on moisture
xfscanf (in,"%le %le", &cps_dry, &cps_wet);
break;
}
default:{
print_err ("unknown model type is required", __FILE__, __LINE__, __func__);
abort();
}
}
//thermal expansion calculation:
switch (betas_type){
case 0:{//constant
xfscanf (in,"%le", &betas0);
break;
}
case 1:{//dependent on moisture
xfscanf (in,"%le %le", &betas_dry, &betas_wet);
break;
}
default:{
print_err ("unknown model type is required", __FILE__, __LINE__, __func__);
abort();
}
}
}
case 1:{//dependent on porosity
xfscanf (in,"%le %le", &bb1, &phi01);
break;
}
default:{
print_err ("unknown model type is required", __FILE__, __LINE__, __func__);
abort();
}
}
 
//relative permebility calculation:
switch (krw_type){
case 0:{//constant
break;
}
case 1:
case 2:{//dependent on saturation degree:
xfscanf (in,"%le %le", &sirr, &ssat);
break;
}
case 3:{//exponential
xfscanf (in,"%le", &lambda_krw);
break;
}
default:{
print_err ("unknown model type is required", __FILE__, __LINE__, __func__);
abort();
}
}
 
//gas relative permebility calculation:
switch (krg_type){
case 0:{//constant
break;
}
case 1:
case 2:{//dependent on saturation degree:
//xfscanf (in,"%le %le", &sirr, &ssat); //not finished
break;
}
case 3:{//exponential
//xfscanf (in,"%le", &lambda_krw); //not finished
break;
}
default:{
print_err ("unknown model type is required", __FILE__, __LINE__, __func__);
abort();
}
}
 
//retention curve type:
switch (sr_type){
case lewis_and_schrefler_sr:{//Lewis and Schrefler's book
lewis_ret.read(in);
break;
}
case van_genuchten_sr:{//partially saturated medium =Van Genuchten model
van_genuchten_ret.read(in);
break;
}
case van_genuchten2_sr:{//partially saturated medium =Van Genuchten model
//van_genuchten_ret.read2(in);
break;
}
case mefel_sr:{//from MEFEL
break;
}
case table_sr:{//from table
//reading of retention curve:
data.read (in);
break;
}
default:{
print_err ("unknown model type is required", __FILE__, __LINE__, __func__);
abort();
}
}
 
//conductivity calculation:
switch (lambda_type){
case 0:{//constant
xfscanf (in,"%le", &lambda0);
break;
}
case 1:{//dependent on moisture
xfscanf (in,"%le %le", &lambda_dry, &lambda_wet);
break;
}
case 2:{//dependent on moisture
xfscanf (in,"%le %le %le %le", &lambda_dry, &lambda_wet, &sr_dry, &sr_wet);
break;
}
default:{
print_err ("unknown model type is required", __FILE__, __LINE__, __func__);
abort();
}
}
 
//specific heat calculation:
switch (cps_type){
case 0:{//constant
xfscanf (in,"%le", &cps0);
break;
}
case 1:{//dependent on moisture
xfscanf (in,"%le %le", &cps_dry, &cps_wet);
break;
}
default:{
print_err ("unknown model type is required", __FILE__, __LINE__, __func__);
abort();
}
}
 
//thermal expansion calculation:
switch (betas_type){
case 0:{//constant
xfscanf (in,"%le", &betas0);
break;
}
case 1:{//dependent on moisture
xfscanf (in,"%le %le", &betas_dry, &betas_wet);
break;
}
default:{
print_err ("unknown model type is required", __FILE__, __LINE__, __func__);
abort();
}
}
}
 
 
444,11 → 463,17
*/
double con_hawf3mat::get_cdiff(double pw,double pg,double t)
{
double cdiff;
double cdiff,pgw;
if(rel_gas_press == 1)
pg = pg + p0;
 
//actualized - check if pgw>pg
pgw = 0.0;
get_pgw(pw,pg,t);
if (pgw > pg)
pg = pgw;
 
//conversion pow(a,b) -> exp(b*log(a))
//cdiff = dv0*p0/pg*pow((t/t0),bv0);
if (t < tcr)
512,10 → 537,16
*/
double con_hawf3mat::get_sw(double pw, double pg, double t, long ipp)
{
double sw;
double sw,pc;
sw = 0.0;
pc = 0.0;
switch (sr_type){
case bazant_sr:{//Bazant
pc = pg - pw;
sw = bazant_ret.sat(pc,t);
break;
}
case lewis_and_schrefler_sr:{//Lewis and Schrefler's book
sw = lewis_ret.sw(pw);
break;
571,10 → 602,15
*/
double con_hawf3mat::get_ds_dpc(double pw, double pg, double t, long ipp)
{
double dsw_dpc;
dsw_dpc = 0.0;
double dsw_dpc,pc;
dsw_dpc = pc = 0.0;
 
switch (sr_type){
case bazant_sr:{//Bazant
pc = pg - pw;
dsw_dpc = bazant_ret.dsat_dpc(pc,t);
break;
}
case lewis_and_schrefler_sr:{//Lewis and Schrefler's book
dsw_dpc = -lewis_ret.dsw_dpw(pw);
break;
620,10 → 656,15
*/
double con_hawf3mat::get_ds_dt(double pw, double pg, double t, long ipp)
{
double dsw_dt;
dsw_dt = 0.0;
double dsw_dt,pc;
dsw_dt = pc = 0.0;
 
switch (sr_type){
case bazant_sr:{//Bazant
pc = pg - pw;
dsw_dt = bazant_ret.dsat_dt(pc,t);
break;
}
case lewis_and_schrefler_sr:{//Lewis and Schrefler's book
dsw_dt = lewis_ret.dsw_dt(pw);
break;
652,6 → 693,7
abort();
}
}
 
return(dsw_dt);
}
 
766,6 → 808,11
krw = pow(sw,lambda_krw);
break;
}
case 9:{//bazant
sw = get_sw(pw,pg,t,ipp);
krw = exp(10.0*log(sw));
break;
}
 
default:{
print_err ("unknown model type is required", __FILE__, __LINE__, __func__);
789,7 → 836,8
*/
double con_hawf3mat::get_krg(double pw, double pg, double t, long ipp)
{
double sw,krg;
double krg;
//double sw;
krg=1.0;
switch (krg_type){
798,17 → 846,17
break;
}
case 1:{//dependent on saturation degree
sw = get_sw(pw,pg,t,ipp);
//sw = get_sw(pw,pg,t,ipp);
krg = 1.0; //not finished
break;
}
case 2:{//dependent on saturation degree
sw = get_sw(pw,pg,t,ipp);
//sw = get_sw(pw,pg,t,ipp);
krg = 1.0; //not finished
break;
}
case 3:{//exponential
sw = get_sw(pw,pg,t,ipp);
//sw = get_sw(pw,pg,t,ipp);
krg = 1.0; //not finished
break;
}
960,21 → 1008,21
double con_hawf3mat::get_dpgw_dpc(double pw,double pg,double t)
{
double dpgw_dpc,pgws,rhow,tt,pc;
 
//critical point of water check
if(t >= tcr)
tt = tcr;
tt = tcr;
else
tt = t;
 
tt = t;
pgws = get_pgws(tt);
rhow = get_rhow(tt);
 
if(rel_gas_press == 1)
pg = pg + p0;
 
pg = pg + p0;
pc = pg - pw;
dpgw_dpc = -pgws*exp(-pc*mw/rhow/gasr/tt)*mw/rhow/gasr/tt;
return(dpgw_dpc);
991,7 → 1039,7
*/
double con_hawf3mat::get_dpgw_dt(double pw,double pg,double t)
{
double dpgw_dt,dpgws_dt,tt,pgws,rhow,pc;
double dpgw_dt,dpgws_dt,tt,pgws,rhow,pc,drhow_dt;
 
//critical point of water check
if(t >= tcr)
1004,8 → 1052,12
pc = pg - pw;
dpgws_dt = get_dpgws_dt(tt);
drhow_dt = get_drhow_dt(tt);
 
if(t < tcr)
dpgw_dt = dpgws_dt*exp(-pc*mw/rhow/gasr/tt) + pgws*exp(-pc*mw/rhow/gasr/tt)*(pc*mw/rhow/gasr/tt/tt);
//actualized form:
//dpgw_dt = dpgws_dt*exp(-pc*mw/rhow/gasr/tt) + pgws*exp(-pc*mw/rhow/gasr/tt)*(pc*mw/rhow/gasr/tt)*(drhow_dt/rhow + 1/tt);
else
dpgw_dt = dpgws_dt*exp(-pc*mw/rhow/gasr/tt);
 
1138,12 → 1190,13
//d2rhow_dtdpc = 0.0;
}
rhow = rhow0;//testing??!!
//rhow = rhow0;//testing??!!
return(rhow);
}
 
 
 
/**
function computes compresibility coefficient of water
1257,6 → 1310,8
drhow_dt = 0.0;
}
drhow_dt = 0.0;//testing??!!
 
return(drhow_dt);
}
 
1487,7 → 1542,6
rhoga = (pg - pgw)*ma/gasr/t;
else
rhoga = 0.0;
//rhoga = (pg - pgw)*ma/gasr/t;
 
return(rhoga);
}
1614,6 → 1668,7
break;
}
case 3:{//
sw = get_sw(pw,pg,t,ipp);
lambdaeff = pow(lambda_wet,sw)*pow(lambda_dry,(1-sw)); //this is from Patek's thesis p. 18, but it seems to be strange
break;
}
1755,7 → 1810,11
k = get_ktg(pw,pg,t,ipp);// *scale_g;//scaling
if((ri == 2) && (ci == 2))
k = get_ktt1(pw,pg,t,ipp);// *scale_t;//scaling
 
//if(Tp->time >= 4.345400e+05)
//fprintf(Outt,"matcond k = %e, ri = %ld, ci = %ld, ipp = %ld\n",k,ri,ci,ipp);
//fflush(Outt);
 
fillm(0.0,d);
d[0][0] = k; d[0][1] = 0.0;
1847,7 → 1906,11
c = get_captg(pw,pg,t,ipp);// *scale_g;//scaling
if((ri == 2) && (ci == 2))
c = get_captt(pw,pg,t,ipp);// *scale_t;//scaling
 
//if(Tp->time >= 4.345400e+05)
//fprintf(Outt,"matcap c = %e, ri = %ld, ci = %ld, ipp = %ld, pw = %e, pg = %e, t = %e\n",c,ri,ci,ipp,pw,pg,t);
//fflush(Outt);
 
check_math_errel(ipp);
}
 
2201,6 → 2264,7
fwu = 0.0;
switch (model_type){
case artificial3:
case lewis_and_schrefler3:
break;
case lewis_and_schrefler3_mefel:{//Lewis and Schrefler's book
2255,6 → 2319,7
fgu = 0.0;
switch (model_type){
case artificial3:
case lewis_and_schrefler3:
break;
case lewis_and_schrefler3_mefel:{//Lewis and Schrefler's book
2306,6 → 2371,7
ftu = 0.0;
switch (model_type){
case artificial3:
case lewis_and_schrefler3:
break;
case lewis_and_schrefler3_mefel:{//Lewis and Schrefler's book
2534,44 → 2600,31
@param ipp - number of integration point
 
TKr, 18.5.2005
TKr 06/12/2022 - actualized
*/
void con_hawf3mat::gaspress_check(double pw,double &pg,double t,long ipp)
{
switch (model_type){
case lewis_and_schrefler3://Lewis and Schrefler's book
case van_genuchten3://partially saturated medium = Van Genuchten model
case lewis_and_schrefler3_mefel:{//Lewis and Schrefler's book, saturation degree ind its derivative are obtained from mefel;
break;
}
case lewis_and_schrefler3_2:{//Lewis and Schrefler's book p. 381
break;
}
default:{
print_err ("unknown model type is required", __FILE__, __LINE__, __func__);
abort();
}
}
 
/* double pgw;
double pgw;
//general
//gas pressure check
pgw = get_pgw(pw,pg,t);
if (pgw > pg)
pg = pgw;
if (pgw > pg){
pg = pgw;
//fprintf(Outt,"gas pressure was modified in integration point No. %ld\n",ipp);
}
 
//general
//fully saturated state
if(pw < 100.0)
pg = 101325.0;
if(pgw < 100.0){
pg = 101325.0;
}
 
// M. Starnoni 24-11-2010
if (pg<0.0){
pg = 0.0;
print_err("gas pressure was modified",__FILE__,__LINE__,__func__);
if (pg < 0.0){
print_err("gas pressure was modified in integration point No. %ld",__FILE__,__LINE__,__func__,ipp);
pg = 0.0;
}
*/
}
 
 
2585,33 → 2638,10
TKr, 18.5.2005
TKr 06/12/2022 - actualized
*/
void con_hawf3mat::waterpress_check(double &pw,double pg,double t,long ipp)
{
switch (model_type){
case lewis_and_schrefler3:{//Lewis and Schrefler's book
//this must be checked:
//if (pw >= 0.0)
//pw = 0.0;
//storing into gauss point
//Tm->ip[ipp].av[0] = pw;
break;
}
case van_genuchten3:{//partially saturated medium = Van Genuchten model
break;
}
case lewis_and_schrefler3_mefel:{//Lewis and Schrefler's book, saturation degree ind its derivative are obtained from mefel;
break;
}
case lewis_and_schrefler3_2:{//Lewis and Schrefler's book p. 381
break;
}
default:{
print_err ("unknown model type is required", __FILE__, __LINE__, __func__);
abort();
}
}
}
 
/**
2620,13 → 2650,13
@param nv - array with variables
*/
void con_hawf3mat::values_correction (vector &nv)
void con_hawf3mat::values_correction (vector &nv, long ipp)
{
// pore water pressure
waterpress_check(nv[0],nv[1],nv[2],0);
waterpress_check(nv[0],nv[1],nv[2],ipp);
// gas pressure
gaspress_check(nv[0],nv[1],nv[2],0);
gaspress_check(nv[0],nv[1],nv[2],ipp);
}
 
 
2644,35 → 2674,39
double capww;
double alpha,n,ks,sw,rhogw,sg,rhow,kw,dpgw_dpc,ds_dpc;
 
//gas pressure check
gaspress_check(pw,pg,t,ipp);
//water pressure check
waterpress_check(pw,pg,t,ipp);
 
alpha = get_alpha(pw,pg,t,ipp);
n = get_porosity(pw,pg,t,ipp);
ks = get_ks(pw,pg,t,ipp);
sw = get_sw(pw,pg,t,ipp);
rhogw = get_rhogw(pw,pg,t);
sg = 1.0 - sw;
rhow = get_rhow(t);
kw = get_kw(pw,pg,t,ipp);
dpgw_dpc = get_dpgw_dpc(pw,pg,t);
ds_dpc = get_ds_dpc(pw,pg,t,ipp);//this is equal to minus derivative with respect to suction
if(rel_gas_press == 1)
pg = pg + p0;
 
if(compress == 1){
//compressible grains:
capww = (alpha - n)/ks*sw*(rhogw*sg + rhow*sw) + rhow*sw*n/kw;
capww = capww - sg*n*mw/t/gasr*dpgw_dpc;
capww = capww - ((alpha - n)/ks*(rhogw*sg*(pw-pg) + rhow*sw*(pw-pg)) + n*(rhow - rhogw))*ds_dpc;
}
if(model_type == artificial3)
capww = capww0;
else{
//incompressible grains:
capww = -sg*n*mw/t/gasr*dpgw_dpc;
capww = capww - n*(rhow - rhogw)*ds_dpc;
//gas pressure check
gaspress_check(pw,pg,t,ipp);
//water pressure check
waterpress_check(pw,pg,t,ipp);
alpha = get_alpha(pw,pg,t,ipp);
n = get_porosity(pw,pg,t,ipp);
ks = get_ks(pw,pg,t,ipp);
sw = get_sw(pw,pg,t,ipp);
rhogw = get_rhogw(pw,pg,t);
sg = 1.0 - sw;
rhow = get_rhow(t);
kw = get_kw(pw,pg,t,ipp);
dpgw_dpc = get_dpgw_dpc(pw,pg,t);
ds_dpc = get_ds_dpc(pw,pg,t,ipp);//this is equal to minus derivative with respect to suction
if(rel_gas_press == 1)
pg = pg + p0;
if(compress == 1){
//compressible grains:
capww = (alpha - n)/ks*sw*(rhogw*sg + rhow*sw) + rhow*sw*n/kw;
capww = capww - sg*n*mw/t/gasr*dpgw_dpc;
capww = capww - ((alpha - n)/ks*(rhogw*sg*(pw-pg) + rhow*sw*(pw-pg)) + n*(rhow - rhogw))*ds_dpc;
}
else{
//incompressible grains:
capww = -sg*n*mw/t/gasr*dpgw_dpc;
capww = capww - n*(rhow - rhogw)*ds_dpc;
}
}
return(capww);
}
2693,35 → 2727,39
double capwg;
double alpha,n,ks,sw,rhogw,sg,rhow,kw,dpgw_dpc,ds_dpc;
 
//gas pressure check
gaspress_check(pw,pg,t,ipp);
//water pressure check
waterpress_check(pw,pg,t,ipp);
 
alpha = get_alpha(pw,pg,t,ipp);
n = get_porosity(pw,pg,t,ipp);
ks = get_ks(pw,pg,t,ipp);
sw = get_sw(pw,pg,t,ipp);
rhogw = get_rhogw(pw,pg,t);
sg = 1.0 - sw;
rhow = get_rhow(t);
kw = get_kw(pw,pg,t,ipp);
dpgw_dpc = get_dpgw_dpc(pw,pg,t);
ds_dpc = get_ds_dpc(pw,pg,t,ipp);//this is equal to minus derivative with respect to suction
if(rel_gas_press == 1)
pg = pg + p0;
 
if(compress == 1){
//compressible grains:
capwg = (alpha - n)/ks*sg*(rhogw*sg + rhow*sw);
capwg = capwg + sg*n*mw/t/gasr*dpgw_dpc;
capwg = capwg + ((alpha - n)/ks*(rhogw*sg*(pw-pg) + rhow*sw*(pw-pg)) + n*(rhow - rhogw))*ds_dpc;
}
if(model_type == artificial3)
capwg = capwg0;
else{
//incompressible grains:
capwg = sg*n*mw/t/gasr*dpgw_dpc;
capwg = capwg + n*(rhow - rhogw)*ds_dpc;
//gas pressure check
gaspress_check(pw,pg,t,ipp);
//water pressure check
waterpress_check(pw,pg,t,ipp);
alpha = get_alpha(pw,pg,t,ipp);
n = get_porosity(pw,pg,t,ipp);
ks = get_ks(pw,pg,t,ipp);
sw = get_sw(pw,pg,t,ipp);
rhogw = get_rhogw(pw,pg,t);
sg = 1.0 - sw;
rhow = get_rhow(t);
kw = get_kw(pw,pg,t,ipp);
dpgw_dpc = get_dpgw_dpc(pw,pg,t);
ds_dpc = get_ds_dpc(pw,pg,t,ipp);//this is equal to minus derivative with respect to suction
if(rel_gas_press == 1)
pg = pg + p0;
if(compress == 1){
//compressible grains:
capwg = (alpha - n)/ks*sg*(rhogw*sg + rhow*sw);
capwg = capwg + sg*n*mw/t/gasr*dpgw_dpc;
capwg = capwg + ((alpha - n)/ks*(rhogw*sg*(pw-pg) + rhow*sw*(pw-pg)) + n*(rhow - rhogw))*ds_dpc;
}
else{
//incompressible grains:
capwg = sg*n*mw/t/gasr*dpgw_dpc;
capwg = capwg + n*(rhow - rhogw)*ds_dpc;
}
}
return(capwg);
}
2741,39 → 2779,42
double capwt;
double betasgw,sw,sg,dpgw_dt,pgw,alpha,n,ks,rhogw,rhow,ds_dt;
 
//gas pressure check
gaspress_check(pw,pg,t,ipp);
//water pressure check
waterpress_check(pw,pg,t,ipp);
 
betasgw = get_betasgw(pw,pg,t,ipp);
sw = get_sw(pw,pg,t,ipp);
sg = 1.0 - sw;
dpgw_dt = get_dpgw_dt(pw,pg,t);
pgw = get_pgw(pw,pg,t);
alpha = get_alpha(pw,pg,t,ipp);
n = get_porosity(pw,pg,t,ipp);
ks = get_ks(pw,pg,t,ipp);
rhogw = get_rhogw(pw,pg,t);
rhow = get_rhow(t);
ds_dt = get_ds_dt(pw,pg,t,ipp);
 
if(rel_gas_press == 1)
pg = pg + p0;
 
if(compress == 1){
//compressible grains:
capwt = -betasgw;
capwt = capwt + sg*n*mw/t/gasr*(dpgw_dt-pgw/t);
capwt = capwt + ((alpha - n)/ks*(rhogw*sg*(pw-pg) + rhow*sw*(pw-pg)) + n*(rhow - rhogw))*ds_dt;
}
if(model_type == artificial3)
capwt = capwt0;
else{
//incompressible grains:
capwt = -betasgw;
capwt = capwt + sg*n*mw/t/gasr*(dpgw_dt-pgw/t);
capwt = capwt + n*(rhow - rhogw)*ds_dt;
//gas pressure check
gaspress_check(pw,pg,t,ipp);
//water pressure check
waterpress_check(pw,pg,t,ipp);
betasgw = get_betasgw(pw,pg,t,ipp);
sw = get_sw(pw,pg,t,ipp);
sg = 1.0 - sw;
dpgw_dt = get_dpgw_dt(pw,pg,t);
pgw = get_pgw(pw,pg,t);
alpha = get_alpha(pw,pg,t,ipp);
n = get_porosity(pw,pg,t,ipp);
ks = get_ks(pw,pg,t,ipp);
rhogw = get_rhogw(pw,pg,t);
rhow = get_rhow(t);
ds_dt = get_ds_dt(pw,pg,t,ipp);
if(rel_gas_press == 1)
pg = pg + p0;
if(compress == 1){
//compressible grains:
capwt = -betasgw;
capwt = capwt + sg*n*mw/t/gasr*(dpgw_dt-pgw/t);
capwt = capwt + ((alpha - n)/ks*(rhogw*sg*(pw-pg) + rhow*sw*(pw-pg)) + n*(rhow - rhogw))*ds_dt;
}
else{
//incompressible grains:
capwt = -betasgw;
capwt = capwt + sg*n*mw/t/gasr*(dpgw_dt-pgw/t);
capwt = capwt + n*(rhow - rhogw)*ds_dt;
}
}
 
return(capwt);
}
 
2791,26 → 2832,29
double kww;
double rhow,krw,muw,rhog,mg,dg,dpgw_dpc,kintr;
 
//gas pressure check
gaspress_check(pw,pg,t,ipp);
//water pressure check
waterpress_check(pw,pg,t,ipp);
 
rhow = get_rhow(t);
krw = get_krw(pw,pg,t,ipp);
muw = get_muw(t);
rhog =get_rhog(pw,pg,t);
mg = get_mg(pw,pg,t);
dg = get_dg(pw,pg,t,ipp);
dpgw_dpc = get_dpgw_dpc(pw,pg,t);
kintr = get_kintr(pw,pg,t,ipp);
 
if(rel_gas_press == 1)
pg = pg + p0;
 
kww = rhow*kintr*krw/muw;//water
kww = kww - rhog*ma*mw/mg/mg*dg/pg*dpgw_dpc;//water vapour
 
if(model_type == artificial3)
kww = kww0;
else{
//gas pressure check
gaspress_check(pw,pg,t,ipp);
//water pressure check
waterpress_check(pw,pg,t,ipp);
rhow = get_rhow(t);
krw = get_krw(pw,pg,t,ipp);
muw = get_muw(t);
rhog =get_rhog(pw,pg,t);
mg = get_mg(pw,pg,t);
dg = get_dg(pw,pg,t,ipp);
dpgw_dpc = get_dpgw_dpc(pw,pg,t);
kintr = get_kintr(pw,pg,t,ipp);
if(rel_gas_press == 1)
pg = pg + p0;
kww = rhow*kintr*krw/muw;//water
kww = kww - rhog*ma*mw/mg/mg*dg/pg*dpgw_dpc;//water vapour
}
return(kww);
}
 
2829,30 → 2873,33
double kwg;
double rhogw,krg,mug,dg,pgw,rhow,krw,muw,rhog,mg,dpgw_dpc,kintr;
//gas pressure check
gaspress_check(pw,pg,t,ipp);
//water pressure check
waterpress_check(pw,pg,t,ipp);
 
rhogw = get_rhogw(pw,pg,t);
krg = get_krg(pw,pg,t,ipp);
mug = get_mug(pw,pg,t);
dg = get_dg(pw,pg,t,ipp);
pgw = get_pgw(pw,pg,t);
rhow = get_rhow(t);
krw = get_krw(pw,pg,t,ipp);
muw = get_muw(t);
rhog = get_rhog(pw,pg,t);
mg = get_mg(pw,pg,t);
dpgw_dpc = get_dpgw_dpc(pw,pg,t);
kintr = get_kintr(pw,pg,t,ipp);
 
if(rel_gas_press == 1)
pg = pg + p0;
 
kwg = rhogw*kintr*krg/mug;//air
kwg = kwg + rhog*ma*mw/mg/mg*dg*(dpgw_dpc/pg - pgw/pg/pg);//water vapour
 
if(model_type == artificial3)
kwg = kwg0;
else{
//gas pressure check
gaspress_check(pw,pg,t,ipp);
//water pressure check
waterpress_check(pw,pg,t,ipp);
rhogw = get_rhogw(pw,pg,t);
krg = get_krg(pw,pg,t,ipp);
mug = get_mug(pw,pg,t);
dg = get_dg(pw,pg,t,ipp);
pgw = get_pgw(pw,pg,t);
rhow = get_rhow(t);
krw = get_krw(pw,pg,t,ipp);
muw = get_muw(t);
rhog = get_rhog(pw,pg,t);
mg = get_mg(pw,pg,t);
dpgw_dpc = get_dpgw_dpc(pw,pg,t);
kintr = get_kintr(pw,pg,t,ipp);
if(rel_gas_press == 1)
pg = pg + p0;
kwg = rhogw*kintr*krg/mug;//air
kwg = kwg + rhog*ma*mw/mg/mg*dg*(dpgw_dpc/pg - pgw/pg/pg);//water vapour
}
return(kwg);
}
 
2870,34 → 2917,38
{
double capgw,alpha,n,ks,sw,sg,pc,rhoga,ds_dpc,dpgw_dpc;
//gas pressure check
gaspress_check(pw,pg,t,ipp);
//water pressure check
waterpress_check(pw,pg,t,ipp);
 
alpha = get_alpha(pw,pg,t,ipp);
n = get_porosity(pw,pg,t,ipp);
ks = get_ks(pw,pg,t,ipp);
sw = get_sw(pw,pg,t,ipp);
sg = 1.0 - sw;
pc = pg - pw;
rhoga = get_rhoga(pw,pg,t);
ds_dpc = get_ds_dpc(pw,pg,t,ipp);
dpgw_dpc = get_dpgw_dpc(pw,pg,t);//this is equal to minus derivative with respect to suction
if(rel_gas_press == 1)
pg = pg + p0;
 
if(compress == 1){
//compressible grains:
capgw = (alpha - n)/ks*sw*sg*rhoga;
capgw = capgw + ((alpha - n)/ks*sg*pc + n)*rhoga*ds_dpc;
capgw = capgw + sg*n*mw/t/gasr*dpgw_dpc;
}
if(model_type == artificial3)
capgw = capgw0;
else{
//incompressible grains:
capgw = n*rhoga*ds_dpc;
capgw = capgw + sg*n*mw/t/gasr*dpgw_dpc;
//gas pressure check
gaspress_check(pw,pg,t,ipp);
//water pressure check
waterpress_check(pw,pg,t,ipp);
alpha = get_alpha(pw,pg,t,ipp);
n = get_porosity(pw,pg,t,ipp);
ks = get_ks(pw,pg,t,ipp);
sw = get_sw(pw,pg,t,ipp);
sg = 1.0 - sw;
pc = pg - pw;
rhoga = get_rhoga(pw,pg,t);
ds_dpc = get_ds_dpc(pw,pg,t,ipp);
dpgw_dpc = get_dpgw_dpc(pw,pg,t);//this is equal to minus derivative with respect to suction
if(rel_gas_press == 1)
pg = pg + p0;
if(compress == 1){
//compressible grains:
capgw = (alpha - n)/ks*sw*sg*rhoga;
capgw = capgw + ((alpha - n)/ks*sg*pc + n)*rhoga*ds_dpc;
capgw = capgw + sg*n*mw/t/gasr*dpgw_dpc;
}
else{
//incompressible grains:
capgw = n*rhoga*ds_dpc;
capgw = capgw + sg*n*mw/t/gasr*dpgw_dpc;
}
}
return(capgw);
}
2917,36 → 2968,40
{
double capgg,alpha,n,ks,sw,sg,pc,rhoga,ds_dpc,dpgw_dpc;
//gas pressure check
gaspress_check(pw,pg,t,ipp);
//water pressure check
waterpress_check(pw,pg,t,ipp);
 
alpha = get_alpha(pw,pg,t,ipp);
n = get_porosity(pw,pg,t,ipp);
ks = get_ks(pw,pg,t,ipp);
sw = get_sw(pw,pg,t,ipp);
sg = 1.0 - sw;
rhoga = get_rhoga(pw,pg,t);
pc = pg - pw;
ds_dpc = get_ds_dpc(pw,pg,t,ipp);//this is equal to minus derivative with respect to suction
dpgw_dpc = get_dpgw_dpc(pw,pg,t);
if(rel_gas_press == 1)
pg = pg + p0;
 
if(compress == 1){
//compressible grains:
capgg = (alpha - n)/ks*sg*sg*rhoga;
capgg = capgg - ((alpha - n)/ks*sg*pc + n)*rhoga*ds_dpc;
capgg = capgg + sg*n*ma/t/gasr;
capgg = capgg - sg*n*mw/t/gasr*dpgw_dpc;
}
if(model_type == artificial3)
capgg = capgg0;
else{
//incompressible grains:
capgg = -n*rhoga*ds_dpc;
capgg = capgg + sg*n*ma/t/gasr;
capgg = capgg - sg*n*mw/t/gasr*dpgw_dpc;
//gas pressure check
gaspress_check(pw,pg,t,ipp);
//water pressure check
waterpress_check(pw,pg,t,ipp);
alpha = get_alpha(pw,pg,t,ipp);
n = get_porosity(pw,pg,t,ipp);
ks = get_ks(pw,pg,t,ipp);
sw = get_sw(pw,pg,t,ipp);
sg = 1.0 - sw;
rhoga = get_rhoga(pw,pg,t);
pc = pg - pw;
ds_dpc = get_ds_dpc(pw,pg,t,ipp);//this is equal to minus derivative with respect to suction
dpgw_dpc = get_dpgw_dpc(pw,pg,t);
if(rel_gas_press == 1)
pg = pg + p0;
if(compress == 1){
//compressible grains:
capgg = (alpha - n)/ks*sg*sg*rhoga;
capgg = capgg - ((alpha - n)/ks*sg*pc + n)*rhoga*ds_dpc;
capgg = capgg + sg*n*ma/t/gasr;
capgg = capgg - sg*n*mw/t/gasr*dpgw_dpc;
}
else{
//incompressible grains:
capgg = -n*rhoga*ds_dpc;
capgg = capgg + sg*n*ma/t/gasr;
capgg = capgg - sg*n*mw/t/gasr*dpgw_dpc;
}
}
return(capgg);
}
2965,39 → 3020,42
{
double capgt,rhoga,betasg,alpha,n,ks,sw,sg,pc,ds_dt,dpgw_dt,pgw;
//gas pressure check
gaspress_check(pw,pg,t,ipp);
//water pressure check
waterpress_check(pw,pg,t,ipp);
 
rhoga = get_rhoga(pw,pg,t);
betasg = get_betasg(pw,pg,t,ipp);
alpha = get_alpha(pw,pg,t,ipp);
n = get_porosity(pw,pg,t,ipp);
ks = get_ks(pw,pg,t,ipp);
sw = get_sw(pw,pg,t,ipp);
sg = 1.0 - sw;
pc = pg - pw;
ds_dt = get_ds_dt(pw,pg,t,ipp);
dpgw_dt = get_dpgw_dt(pw,pg,t);
pgw = get_pgw(pw,pg,t);
 
if(rel_gas_press == 1)
pg = pg + p0;
 
if(compress == 1){
//compressible grains:
capgt = -1.0*rhoga*betasg;
capgt = capgt - ((alpha - n)/ks*sg*pc + n)*rhoga*ds_dt;
capgt = capgt - (sg*n*ma/t/t/gasr*pg + sg*n*ma/t/gasr*(dpgw_dt-pgw/t));
}
if(model_type == artificial3)
capgt = capgt0;
else{
//incompressible grains:
capgt = -1.0*rhoga*betasg;
capgt = capgt - n*rhoga*ds_dt;
capgt = capgt - (sg*n*ma/t/t/gasr*pg + sg*n*ma/t/gasr*(dpgw_dt-pgw/t));
//gas pressure check
gaspress_check(pw,pg,t,ipp);
//water pressure check
waterpress_check(pw,pg,t,ipp);
rhoga = get_rhoga(pw,pg,t);
betasg = get_betasg(pw,pg,t,ipp);
alpha = get_alpha(pw,pg,t,ipp);
n = get_porosity(pw,pg,t,ipp);
ks = get_ks(pw,pg,t,ipp);
sw = get_sw(pw,pg,t,ipp);
sg = 1.0 - sw;
pc = pg - pw;
ds_dt = get_ds_dt(pw,pg,t,ipp);
dpgw_dt = get_dpgw_dt(pw,pg,t);
pgw = get_pgw(pw,pg,t);
if(rel_gas_press == 1)
pg = pg + p0;
if(compress == 1){
//compressible grains:
capgt = -1.0*rhoga*betasg;
capgt = capgt - ((alpha - n)/ks*sg*pc + n)*rhoga*ds_dt;
capgt = capgt - (sg*n*ma/t/t/gasr*pg + sg*n*ma/t/gasr*(dpgw_dt-pgw/t));
}
else{
//incompressible grains:
capgt = -1.0*rhoga*betasg;
capgt = capgt - n*rhoga*ds_dt;
capgt = capgt - (sg*n*ma/t/t/gasr*pg + sg*n*ma/t/gasr*(dpgw_dt-pgw/t));
}
}
 
return(capgt);
}
 
3016,21 → 3074,24
{
double kgw,rhog,dg,mg,dpgw_dpc;
//gas pressure check
gaspress_check(pw,pg,t,ipp);
//water pressure check
waterpress_check(pw,pg,t,ipp);
rhog = get_rhog(pw,pg,t);
dg = get_dg(pw,pg,t,ipp);
dpgw_dpc = get_dpgw_dpc(pw,pg,t);
mg = get_mg(pw,pg,t);
 
if(rel_gas_press == 1)
pg = pg + p0;
kgw = rhog*ma*mw/mg/mg*dg/pg*dpgw_dpc;//water vapour
 
if(model_type == artificial3)
kgw = kgw0;
else{
//gas pressure check
gaspress_check(pw,pg,t,ipp);
//water pressure check
waterpress_check(pw,pg,t,ipp);
rhog = get_rhog(pw,pg,t);
dg = get_dg(pw,pg,t,ipp);
dpgw_dpc = get_dpgw_dpc(pw,pg,t);
mg = get_mg(pw,pg,t);
if(rel_gas_press == 1)
pg = pg + p0;
kgw = rhog*ma*mw/mg/mg*dg/pg*dpgw_dpc;//water vapour
}
return(kgw);
}
 
3049,27 → 3110,30
{
double kgg,rhoga,krg,mug,rhog,dg,mg,dpgw_dpc,pgw,kintr;
//gas pressure check
gaspress_check(pw,pg,t,ipp);
//water pressure check
waterpress_check(pw,pg,t,ipp);
 
rhoga = get_rhoga(pw,pg,t);
krg = get_krg(pw,pg,t,ipp);
mug = get_mug(pw,pg,t);
rhog = get_rhog(pw,pg,t);
dg = get_dg(pw,pg,t,ipp);
mg = get_mg(pw,pg,t);
dpgw_dpc = get_dpgw_dpc(pw,pg,t);
pgw = get_pgw(pw,pg,t);
kintr = get_kintr(pw,pg,t,ipp);
 
if(rel_gas_press == 1)
pg = pg + p0;
 
kgg = rhoga*kintr*krg/mug;//air
kgg = kgg - rhog*ma*mw/mg/mg*dg*(dpgw_dpc/pg-pgw/pg/pg);//water vapour
 
if(model_type == artificial3)
kgg = kgg0;
else{
//gas pressure check
gaspress_check(pw,pg,t,ipp);
//water pressure check
waterpress_check(pw,pg,t,ipp);
rhoga = get_rhoga(pw,pg,t);
krg = get_krg(pw,pg,t,ipp);
mug = get_mug(pw,pg,t);
rhog = get_rhog(pw,pg,t);
dg = get_dg(pw,pg,t,ipp);
mg = get_mg(pw,pg,t);
dpgw_dpc = get_dpgw_dpc(pw,pg,t);
pgw = get_pgw(pw,pg,t);
kintr = get_kintr(pw,pg,t,ipp);
if(rel_gas_press == 1)
pg = pg + p0;
kgg = rhoga*kintr*krg/mug;//air
kgg = kgg - rhog*ma*mw/mg/mg*dg*(dpgw_dpc/pg-pgw/pg/pg);//water vapour
}
return(kgg);
}
 
3108,33 → 3172,38
double con_hawf3mat::get_captw(double pw,double pg,double t,long ipp)
{
double captw,dhvap,rhow,alpha,n,ks,sw,kw,dsw_dpc;
//gas pressure check
gaspress_check(pw,pg,t,ipp);
//water pressure check
waterpress_check(pw,pg,t,ipp);
 
dhvap = get_dhvap(t);
rhow = get_rhow(t);
alpha = get_alpha(pw,pg,t,ipp);
n = get_porosity(pw,pg,t,ipp);
ks = get_ks(pw,pg,t,ipp);
sw = get_sw(pw,pg,t,ipp);
kw = get_kw(pw,pg,t,ipp);
dsw_dpc = get_ds_dpc(pw,pg,t,ipp);//this is equal to minus derivative with respect to suction
 
if(rel_gas_press == 1)
pg = pg + p0;
 
if(compress == 1){
//compressible grains:
captw = -1.0*dhvap*rhow*((alpha - n)/ks*sw*sw + sw*n/kw);
captw = captw + dhvap*rhow*((alpha - n)/ks*sw*pw - (alpha - n)/ks*sw*pg + n)*dsw_dpc;
}
if(model_type == artificial3)
captw = captw0;
else{
//incompressible grains:
captw = dhvap*rhow*n*dsw_dpc;
//gas pressure check
gaspress_check(pw,pg,t,ipp);
//water pressure check
waterpress_check(pw,pg,t,ipp);
dhvap = get_dhvap(t);
rhow = get_rhow(t);
alpha = get_alpha(pw,pg,t,ipp);
n = get_porosity(pw,pg,t,ipp);
ks = get_ks(pw,pg,t,ipp);
sw = get_sw(pw,pg,t,ipp);
kw = get_kw(pw,pg,t,ipp);
dsw_dpc = get_ds_dpc(pw,pg,t,ipp);//this is equal to minus derivative with respect to suction
if(rel_gas_press == 1)
pg = pg + p0;
if(compress == 1){
//compressible grains:
captw = -1.0*dhvap*rhow*((alpha - n)/ks*sw*sw + sw*n/kw);
captw = captw + dhvap*rhow*((alpha - n)/ks*sw*pw - (alpha - n)/ks*sw*pg + n)*dsw_dpc;
}
else{
//incompressible grains:
captw = dhvap*rhow*n*dsw_dpc;
}
}
 
return(captw);
}
 
3190,34 → 3259,44
{
double captg,dhvap,rhow,alpha,n,ks,sw,sg,kw,dsw_dpc;
//gas pressure check
gaspress_check(pw,pg,t,ipp);
//water pressure check
waterpress_check(pw,pg,t,ipp);
if(model_type == artificial3)
captg = captg0;
else{
//gas pressure check
gaspress_check(pw,pg,t,ipp);
//water pressure check
waterpress_check(pw,pg,t,ipp);
dhvap = get_dhvap(t);
rhow = get_rhow(t);
alpha = get_alpha(pw,pg,t,ipp);
n = get_porosity(pw,pg,t,ipp);
ks = get_ks(pw,pg,t,ipp);
sw = get_sw(pw,pg,t,ipp);
sg = 1.0 - sw;
kw = get_kw(pw,pg,t,ipp);
dsw_dpc = get_ds_dpc(pw,pg,t,ipp);//this is equal to minus derivative with respect to suction
if(rel_gas_press == 1)
pg = pg + p0;
//if(Tp->time >= 4.345400e+05 && ipp == 25188){
//fprintf(Outt,"dhvap = %e, rhow = %e, alpha = %e, n = %e, ks = %e, sw = %e, sg = %e, kw = %e, dsw_dpc = %e, pw = %e, pg = %e, t = %e\n",
// dhvap,rhow,alpha,n,ks,sw,sg,kw,dsw_dpc,pw,pg,t);
//fflush(Outt);
//}
 
dhvap = get_dhvap(t);
rhow = get_rhow(t);
alpha = get_alpha(pw,pg,t,ipp);
n = get_porosity(pw,pg,t,ipp);
ks = get_ks(pw,pg,t,ipp);
sw = get_sw(pw,pg,t,ipp);
sg = 1.0 - sw;
kw = get_kw(pw,pg,t,ipp);
dsw_dpc = get_ds_dpc(pw,pg,t,ipp);//this is equal to minus derivative with respect to suction
 
if(rel_gas_press == 1)
pg = pg + p0;
 
if(compress == 1){
//compressible grains:
captg = -dhvap*rhow*(alpha - n)/ks*sg*sw;//water
captg = captg - dhvap*rhow*((alpha - n)/ks*sw*pw - (alpha - n)/ks*sw*pg + n)*dsw_dpc;
if(compress == 1){
//compressible grains:
captg = -dhvap*rhow*(alpha - n)/ks*sg*sw;//water
captg = captg - dhvap*rhow*((alpha - n)/ks*sw*pw - (alpha - n)/ks*sw*pg + n)*dsw_dpc;
}
else{
//incompressible grains:
captg = -dhvap*rhow*n*dsw_dpc;
}
}
else{
//incompressible grains:
captg = -dhvap*rhow*n*dsw_dpc;
}
 
return(captg);
}
 
3235,35 → 3314,39
{
double captt,rhocp,dhvap,betasw,rhow,alpha,n,ks,sw,dsw_dt;
//gas pressure check
gaspress_check(pw,pg,t,ipp);
//water pressure check
waterpress_check(pw,pg,t,ipp);
 
rhocp = get_rhocp(pw,pg,t,ipp);
dhvap = get_dhvap(t);
betasw = get_betasw(pw,pg,t,ipp);
rhow = get_rhow(t);
alpha = get_alpha(pw,pg,t,ipp);
n = get_porosity(pw,pg,t,ipp);
ks = get_ks(pw,pg,t,ipp);
sw = get_sw(pw,pg,t,ipp);
dsw_dt = get_ds_dt(pw,pg,t,ipp);
 
if(rel_gas_press == 1)
pg = pg + p0;
 
if(compress == 1){
//compressible grains:
captt = rhocp;
captt = captt + dhvap*betasw*rhow;
captt = captt - dhvap*(rhow*((alpha - n)/ks*sw*pw - (alpha - n)/ks*sw*pg + n)*dsw_dt);
}
if(model_type == artificial3)
captt = captt0;
else{
//incompressible grains:
captt = rhocp;
captt = captt + dhvap*betasw*rhow;
captt = captt - dhvap*n*dsw_dt;
//gas pressure check
gaspress_check(pw,pg,t,ipp);
//water pressure check
waterpress_check(pw,pg,t,ipp);
rhocp = get_rhocp(pw,pg,t,ipp);
dhvap = get_dhvap(t);
betasw = get_betasw(pw,pg,t,ipp);
rhow = get_rhow(t);
alpha = get_alpha(pw,pg,t,ipp);
n = get_porosity(pw,pg,t,ipp);
ks = get_ks(pw,pg,t,ipp);
sw = get_sw(pw,pg,t,ipp);
dsw_dt = get_ds_dt(pw,pg,t,ipp);
if(rel_gas_press == 1)
pg = pg + p0;
if(compress == 1){
//compressible grains:
captt = rhocp;
captt = captt + dhvap*betasw*rhow;
captt = captt - dhvap*(rhow*((alpha - n)/ks*sw*pw - (alpha - n)/ks*sw*pg + n)*dsw_dt);
}
else{
//incompressible grains:
captt = rhocp;
captt = captt + dhvap*betasw*rhow;
captt = captt - dhvap*n*dsw_dt;
}
}
return(captt);
}
3283,19 → 3366,24
{
double dhvap,rhow,krw,muw,ktw,kintr;
//gas pressure check
gaspress_check(pw,pg,t,ipp);
//water pressure check
waterpress_check(pw,pg,t,ipp);
dhvap = get_dhvap(t);
rhow = get_rhow(t);
krw = get_krw(pw,pg,t,ipp);
muw = get_muw(t);
kintr = get_kintr(pw,pg,t,ipp);
if(model_type == artificial3)
ktw = ktw0;
else{
//gas pressure check
gaspress_check(pw,pg,t,ipp);
//water pressure check
waterpress_check(pw,pg,t,ipp);
dhvap = get_dhvap(t);
rhow = get_rhow(t);
krw = get_krw(pw,pg,t,ipp);
muw = get_muw(t);
kintr = get_kintr(pw,pg,t,ipp);
//correction 02/12/2022
ktw = -dhvap*(rhow*kintr*krw/muw);//water
}
 
ktw = dhvap*(rhow*kintr*krw/muw);//water
 
return(ktw);
}
 
3315,15 → 3403,19
{
double lambdaeff,ktt1;
//gas pressure check
gaspress_check(pw,pg,t,ipp);
//water pressure check
waterpress_check(pw,pg,t,ipp);
if(model_type == artificial3)
ktt1 = ktt0;
else{
//gas pressure check
gaspress_check(pw,pg,t,ipp);
//water pressure check
waterpress_check(pw,pg,t,ipp);
lambdaeff = get_lambdaeff(pw,pg,t,ipp);
ktt1 = lambdaeff;
}
 
lambdaeff = get_lambdaeff(pw,pg,t,ipp);
 
ktt1 = lambdaeff;
 
return(ktt1);
}
 
3341,21 → 3433,24
{
double n,sw,cpw,rhow,krw,muw,ktt2a,kintr;
//gas pressure check
gaspress_check(pw,pg,t,ipp);
//water pressure check
waterpress_check(pw,pg,t,ipp);
 
n = get_porosity(pw,pg,t,ipp);
sw = get_sw(pw,pg,t,ipp);
rhow = get_rhow(t);
cpw = get_cpw();
krw = get_krw(pw,pg,t,ipp);
muw = get_muw(t);
kintr = get_kintr(pw,pg,t,ipp);
 
ktt2a = n*sw*rhow*cpw*kintr*krw/muw;
if(model_type == artificial3)
ktt2a = 0.0;
else{
//gas pressure check
gaspress_check(pw,pg,t,ipp);
//water pressure check
waterpress_check(pw,pg,t,ipp);
n = get_porosity(pw,pg,t,ipp);
sw = get_sw(pw,pg,t,ipp);
rhow = get_rhow(t);
cpw = get_cpw();
krw = get_krw(pw,pg,t,ipp);
muw = get_muw(t);
kintr = get_kintr(pw,pg,t,ipp);
ktt2a = n*sw*rhow*cpw*kintr*krw/muw;
}
return(ktt2a);
}
 
3373,22 → 3468,26
double n,sw,sg,rhog,cpg,krg,mug,kintr;
double ktt2b;
//gas pressure check
gaspress_check(pw,pg,t,ipp);
//water pressure check
waterpress_check(pw,pg,t,ipp);
if(model_type == artificial3)
ktt2b = 0.0;
else{
//gas pressure check
gaspress_check(pw,pg,t,ipp);
//water pressure check
waterpress_check(pw,pg,t,ipp);
n = get_porosity(pw,pg,t,ipp);
sw = get_sw(pw,pg,t,ipp);
sg = 1.0 - sw;
rhog = get_rhog(pw,pg,t);
cpg = get_cpg(pw,pg,t);
krg = get_krg(pw,pg,t,ipp);
mug = get_mug(pw,pg,t);
kintr = get_kintr(pw,pg,t,ipp);
ktt2b = n*sg*rhog*cpg*kintr*krg/mug;
}
 
n = get_porosity(pw,pg,t,ipp);
sw = get_sw(pw,pg,t,ipp);
sg = 1.0 - sw;
rhog = get_rhog(pw,pg,t);
cpg = get_cpg(pw,pg,t);
krg = get_krg(pw,pg,t,ipp);
mug = get_mug(pw,pg,t);
kintr = get_kintr(pw,pg,t,ipp);
 
ktt2b = n*sg*rhog*cpg*kintr*krg/mug;
 
return(ktt2b);
}
 
3406,14 → 3505,18
{
double rhow,ktt2c;
//gas pressure check
gaspress_check(pw,pg,t,ipp);
//water pressure check
waterpress_check(pw,pg,t,ipp);
rhow = get_rhow(t);
ktt2c = rhow;
if(model_type == artificial3)
ktt2c = 0.0;
else{
//gas pressure check
gaspress_check(pw,pg,t,ipp);
//water pressure check
waterpress_check(pw,pg,t,ipp);
rhow = get_rhow(t);
ktt2c = rhow;
}
 
return(ktt2c);
}
3432,15 → 3535,19
{
double rhog,ktt2d;
//gas pressure check
gaspress_check(pw,pg,t,ipp);
//water pressure check
waterpress_check(pw,pg,t,ipp);
if(model_type == artificial3)
ktt2d = 0.0;
else{
//gas pressure check
gaspress_check(pw,pg,t,ipp);
//water pressure check
waterpress_check(pw,pg,t,ipp);
rhog = get_rhog(pw,pg,t);
ktt2d = rhog;
}
 
rhog = get_rhog(pw,pg,t);
ktt2d = rhog;
 
return(ktt2d);
}
 
3459,22 → 3566,25
double fc1;
double rhow,rhogw,rhog,krw,muw,krg,mug,kintr;
//gas pressure check
gaspress_check(pw,pg,t,ipp);
//water pressure check
waterpress_check(pw,pg,t,ipp);
 
rhow = get_rhow(t);
rhogw = get_rhogw(pw,pg,t);
rhog = get_rhog(pw,pg,t);
krw = get_krw(pw,pg,t,ipp);
muw = get_muw(t);
krg = get_krg(pw,pg,t,ipp);
mug = get_mug(pw,pg,t);
kintr = get_kintr(pw,pg,t,ipp);
 
fc1 = rhogw*kintr*krg*rhog/mug + rhow*kintr*krw*rhow/muw;
if(model_type == artificial3)
fc1 = 0.0;
else{
//gas pressure check
gaspress_check(pw,pg,t,ipp);
//water pressure check
waterpress_check(pw,pg,t,ipp);
rhow = get_rhow(t);
rhogw = get_rhogw(pw,pg,t);
rhog = get_rhog(pw,pg,t);
krw = get_krw(pw,pg,t,ipp);
muw = get_muw(t);
krg = get_krg(pw,pg,t,ipp);
mug = get_mug(pw,pg,t);
kintr = get_kintr(pw,pg,t,ipp);
fc1 = rhogw*kintr*krg*rhog/mug + rhow*kintr*krw*rhow/muw;
}
return(fc1);
}
 
3492,19 → 3602,24
double fg;
double rhoga,rhog,krg,mug,kintr;
//gas pressure check
gaspress_check(pw,pg,t,ipp);
//water pressure check
waterpress_check(pw,pg,t,ipp);
if(model_type == artificial3)
fg = 0.0;
else{
//gas pressure check
gaspress_check(pw,pg,t,ipp);
//water pressure check
waterpress_check(pw,pg,t,ipp);
rhoga = get_rhoga(pw,pg,t);
rhog = get_rhog(pw,pg,t);
krg = get_krg(pw,pg,t,ipp);
mug = get_mug(pw,pg,t);
kintr = get_kintr(pw,pg,t,ipp);
fg = rhoga*kintr*krg*rhog/mug;
fg = 0.0;//no gravity force is included for the moist air
}
 
rhoga = get_rhoga(pw,pg,t);
rhog = get_rhog(pw,pg,t);
krg = get_krg(pw,pg,t,ipp);
mug = get_mug(pw,pg,t);
kintr = get_kintr(pw,pg,t,ipp);
fg = rhoga*kintr*krg*rhog/mug;
fg = 0.0;//no gravity force is included for the moist air
return(fg);
}
 
3522,19 → 3637,23
double ft1;
double dhvap,rhow,krw,muw,kintr;
 
//gas pressure check
gaspress_check(pw,pg,t,ipp);
//water pressure check
waterpress_check(pw,pg,t,ipp);
if(model_type == artificial3)
ft1 = 0.0;
else{
//gas pressure check
gaspress_check(pw,pg,t,ipp);
//water pressure check
waterpress_check(pw,pg,t,ipp);
dhvap = get_dhvap(t);
rhow = get_rhow(t);
krw = get_krw(pw,pg,t,ipp);
muw = get_muw(t);
kintr = get_kintr(pw,pg,t,ipp);
ft1 = -dhvap*rhow*kintr*krw*rhow/muw;
}
 
dhvap = get_dhvap(t);
rhow = get_rhow(t);
krw = get_krw(pw,pg,t,ipp);
muw = get_muw(t);
kintr = get_kintr(pw,pg,t,ipp);
 
ft1 = -dhvap*rhow*kintr*krw*rhow/muw;
 
return(ft1);
}
 
3567,7 → 3686,7
pgws = get_pgws(t);
rhow = get_rhow(t);
trc=-1.0*pgws*mw*mw/gasr/gasr/t/t/rhow;//1st member of T. serie
trc=1.0*pgws*mw*mw/gasr/gasr/t/t/rhow;//1st member of T. serie
break;
}
3651,7 → 3770,7
//pgws = get_pgws(t);
//rhow = get_rhow(t);
 
trc = 1.0;//get_pwrh(bv,t);//tady??!!
trc = -1.0;//get_pwrh(bv,t);//tady??!!
//trc=-1.0*pgws*mw*mw/gasr/gasr/t/t/rhow;//1st member of T. serie
 
break;
3728,7 → 3847,8
double con_hawf3mat::get_transmission_nodval_ww(double bv,double pw,double pg,double t,long bc,long ipp)
{
double new_nodval;
state_eq tt;
 
//gas pressure check
gaspress_check(pw,pg,t,ipp);
//water pressure check
3736,7 → 3856,8
 
switch (bc){//type of prescribed variable
case 30:{
new_nodval = 0.0;
new_nodval = tt.get_pcrh(bv,t);
new_nodval = pg - new_nodval;
break;
}
case 40:{//simulation of free soil surface
3765,7 → 3886,8
*/
double con_hawf3mat::get_transmission_flux_ww(double bv,double pw,double pg,double t,long bc,long ipp)
{
double flux,trc;
double flux,trc,pc;
state_eq tt;
//gas pressure check
gaspress_check(pw,pg,t,ipp);
3774,7 → 3896,11
 
switch (bc){//type of prescribed variable
case 30:{//for testing - boundary flux
flux = 0.0;
pc = pg - pw;
 
flux = tt.get_pcrh(bv,t);
flux = -1.0*(flux - pc);
 
break;
}
case 40:{//simulation of free soil surface
3938,8 → 4064,12
other = r[0];
break;
}
case 6:{//moisture content
other = get_w(r[0],r[1],r[2],ipp);
//case 6:{//moisture content
//other = get_w(r[0],r[1],r[2],ipp);
//break;
//}
case 6:{//relative humidity
other = get_pgw(r[0],r[1],r[2])/get_pgws(r[2]);
break;
}
case 7:{//dsw_dpw
3991,8 → 4121,12
fprintf (out,"Pore water pressure (Pa) ");
break;
}
case 6:{//moisture content
fprintf (out,"Moisture content (kg/kg) ");
//case 6:{//moisture content
//fprintf (out,"Moisture content (kg/kg) ");
//break;
//}
case 6:{//relative humidity
fprintf (out,"Relative humidity () ");
break;
}
case 7:{//dS_dpw
4137,25 → 4271,6
 
pg = Tm->ip[ipp].av[1];
 
switch (model_type){
case lewis_and_schrefler3:{//Lewis and Schrefler's book
//this must be checked:
if (pg <= 0.0)
pg = 0.0;
//storing into gauss point
Tm->ip[ipp].av[1] = pg;
break;
}
case van_genuchten3://partially saturated medium = Van Genuchten model
case lewis_and_schrefler3_mefel:{//Lewis and Schrefler's book, saturation degree ind its derivative are obtained from mefel
break;
}
default:{
print_err ("unknown model type is required", __FILE__, __LINE__, __func__);
abort();
}
}
return(pg);
}
 
4175,6 → 4290,9
pw = Tm->ip[ipp].av[0];
pg = Tm->ip[ipp].av[1];
if(rel_gas_press == 1)
pg = pg + p0;
 
switch (model_type){
case lewis_and_schrefler3:{//Lewis and Schrefler's book
suction = -1.0*(pg - pw);//this is correct, because pore water pressure is negative
4309,6 → 4427,7
void con_hawf3mat::give_reqntq(long *antq)
{
switch (model_type){
case artificial3:
case lewis_and_schrefler3:{//Lewis and Schrefler's book
break;
}
/trunk/SIFEL/TRFEL/SRC/consol_hawf3.h
3,6 → 3,8
 
#include "aliast.h"
#include "genfile.h"
#include "baroghel_reten.h"
#include "bazant_reten.h"
#include "lewis_schrefler.h"
#include "van_genuchten.h"
 
93,7 → 95,7
void gaspress_check(double pw,double &pg,double t,long ipp);
void waterpress_check(double &pw,double pg,double t,long ipp);
 
void values_correction (vector &nv);
void values_correction (vector &nv, long ipp);
 
double get_kww(double pw,double pg,double t,long ipp);
double get_capww(double pw,double pg,double t,long ipp);
154,6 → 156,10
/// marks required non-transport quantities
void give_reqntq(long *antq);
 
/// Baroghel retention curve:
baroghel_reten baroghel_ret;
/// Bazant retention curve:
bazant_reten bazant_ret;
/// Lewis and Schrefler's retention curve:
lewis_reten lewis_ret;
/// van Genuchten's retention curve:
231,7 → 237,9
double sirr,ssat,lambda_krw;
double bb1,phi01; //permeability parameters
 
 
//parameters for the arificial material:
double kww0,kwg0,kwt0,kgw0,kgg0,kgt0,ktw0,ktg0,ktt0;
double capww0,capwg0,capwt0,capgw0,capgg0,capgt0,captw0,captg0,captt0;
};
 
#endif
/trunk/SIFEL/TRFEL/SRC/consol_hwf2.cpp
1432,21 → 1432,6
*/
void con_hwf2mat::waterpress_check(double &pw,double t,long ipp)
{
switch (model_type){
case lewis_and_schrefler2hw:{//Lewis and Schrefler's book
break;
}
case lewis_and_schrefler2hw_mefel:{//Lewis and Schrefler's book, saturation degree ind its derivative are obtained from mefel;
break;
}
case van_genuchten2hw:{//partially saturated medium = Van Genuchten model
break;
}
default:{
print_err ("unknown model type is required", __FILE__, __LINE__, __func__);
abort();
}
}
}
 
/**
1455,7 → 1440,7
@param nv - array with variables
*/
void con_hwf2mat::values_correction (vector &nv)
void con_hwf2mat::values_correction (vector &nv,long ipp)
{
// pore water pressure
waterpress_check(nv[0],nv[1],0);
/trunk/SIFEL/TRFEL/SRC/consol_hwf2.h
64,7 → 64,7
 
void waterpress_check(double &pw,double t,long ipp);
 
void values_correction (vector &nv);
void values_correction (vector &nv,long ipp);
 
double get_kww(double pw,double t,long ipp);
double get_capww(double pw,double t,long ipp);
/trunk/SIFEL/TRFEL/SRC/consol_wf2.cpp
34,6 → 34,8
muw0 = 1.0e-3;//water viscosity [Pa.s]
kw = 2.0e9;//bulk modulus of water [Pa]
mug0 = 1.8e-5;//air viscosity [Pa.s]
 
//PHYSICAL PROPERTIES OF SOIL set to zero
alpha = 1.0;//Biot's constant [-] alpha = 1 - kt/ks
ks = 2.167e9;//bulk modulus of solid phase (grains) [Pa]
72,7 → 74,6
switch (model_type){
case lewis_and_schrefler2:{//Lewis and Schrefler's model approach
lewis_ret.read(in);
break;
}
case lewis_and_schrefler2_mefel:{//Lewis and Schrefler's approach coupled with mechanics
1027,7 → 1028,7
rhog = 1.25;//temporarilly
kintr = get_kintr(pw,pg,ipp);
fg1 = krg*kintr/mug0*rhog;//p. 97
//fg1 = 0.0;//no gravity force is included for the air
fg1 = 0.0;//no gravity force is included for the air
return(fg1);
 
}
/trunk/SIFEL/TRFEL/SRC/constrel.cpp
130,6 → 130,9
cpw = 4181.0;
lambdaw = 0.6;
kw0 = 0.43e10;//compresibility coefficient of water
 
//minimal pc value
pcmin = 100.0;
}
state_eq::~state_eq()
{}
/trunk/SIFEL/TRFEL/SRC/constrel.h
97,6 → 97,8
double get_hydren(double pc,double pg,double t,long ipp);
double get_fste(double pc,double pg,double t,long ipp);
 
double pcmin; //minimal value of capillary pressure
 
private:
 
double gasr;//universal gas constant
/trunk/SIFEL/TRFEL/SRC/multiphase.cpp
215,7 → 215,7
d[0][0] = k; d[0][1] = 0.0;
d[1][0] = 0.0; d[1][1] = k;
 
fprintf(Outt,"\n i=%ld j=%ld kij=%e",ri,ci,k);
//fprintf(Outt,"\n i=%ld j=%ld kij=%e",ri,ci,k);//debug information
 
}
 
839,6 → 839,8
*/
void multiph::cappress_check(double &pc,double pg,double t,long ipp)
{
state_eq tt;
 
/* double pgw,pcmax,rhow;
state_eq tt;
867,9 → 869,10
*/
 
// M. Starnoni 24-11-2010
if (pc<0.0){
pc=0.0;
print_err("capillary pressure was modified",__FILE__,__LINE__,__func__);
//corrected by TKr 08/12/2022
if (pc <= tt.pcmin){
pc = tt.pcmin;
//print_err("capillary pressure was modified",__FILE__,__LINE__,__func__);
}
}
936,11 → 939,12
ds_dpc = tt.get_ds_dpc(pc,t,ipp);
ddbw = tt.get_ddbw(pc,pg,t,ipp);
 
kcc = rhow*kintr*krw/muw + ddbw*rhow;//water
kcc = kcc - deff*dpgw_dpc/t/mg;//water vapour
//kcc = rhow*kintr*krw/muw + ddbw*rhow;//water
//kcc = kcc - deff*dpgw_dpc/t/mg;//water vapour
//kcc = -1.0*kcc;//sign correction
//correction by TKr 09/12/2022
kcc = -rhow*kintr*krw/muw - ddbw*rhow + rhog*ma*mw/mg/mg*deff/pg*dpgw_dpc;//water + bound water + water vapour
 
kcc = -1.0*kcc;//sign correction
 
return(kcc);
}
 
989,6 → 993,7
double multiph::get_kcg(double pc,double pg,double t,long ipp)
{
double rhogw,kintr,krg,mug,deff,pgw,rhow,krw,muw,kcg,mg;
double rhog;
state_eq tt;
//gas pressure check
1006,11 → 1011,13
krw = tt.get_krw(pc,t,ipp);
muw = tt.get_muw(t);
mg = tt.get_mg(pc,pg,t);
rhog = tt.get_rhog(pc,pg,t);
 
kcg = -1.0*rhow*kintr*krw/muw;//water
kcg = -1.0*rhogw*kintr*krg/mug + deff*pgw/pg/t/mg;//water vapour
kcg = -1.0*kcg;//sign correction
//kcg = -1.0*rhow*kintr*krw/muw;//water
//kcg = -1.0*rhogw*kintr*krg/mug + deff*pgw/pg/t/mg;//water vapour
//kcg = -1.0*kcg;//sign correction
//correction by TKr 09/12/2022
kcg = rhow*kintr*krw/muw + rhogw*kintr*krg/mug - rhog*ma*mw/mg/mg*deff*pgw/pg/pg;//water + water vapour
 
return(kcg);
}
1067,10 → 1074,11
dpgw_dt = tt.get_dpgw_dt(pc,t);
mg = tt.get_mg(pc,pg,t);
 
kct = 0.0;//water
kct = kct - deff*dpgw_dt/t/mg;//water vapour
kct = -1.0*kct;//sign correction
//kct = 0.0;//water
//kct = kct - deff*dpgw_dt/t/mg;//water vapour
// kct = -1.0*kct;//sign correction
//correction by TKr 09/12/2022
kct = rhog*ma*mw/mg/mg*deff/pg*dpgw_dt;
 
return(kct);
}
1124,6 → 1132,7
double multiph::get_kgg(double pc,double pg,double t,long ipp)
{
double rhoga,kintr,krg,mug,rhog,deff,kgg,mg,pgw;
double dpgw_dpc;
state_eq tt;
//gas pressure check
1140,9 → 1149,12
mg = tt.get_mg(pc,pg,t);
pgw = tt.get_pgw(pc,t);
 
kgg = -1.0*rhoga*kintr*krg/mug - deff*pgw/pg/t/mg;//dry air + water vapour
//kgg = -1.0*rhoga*kintr*krg/mug - deff*pgw/pg/t/mg;//dry air + water vapour
//kgg = -1.0*kgg;//sign correction
 
kgg = -1.0*kgg;//sign correction
// correction 09/12/2022:
dpgw_dpc = tt.get_dpgw_dpc(pc,t);
kgg = rhoga*kintr*krg/mug - rhog*ma*mw/mg/mg*deff*(dpgw_dpc/pg-pgw/pg/pg);
 
return(kgg);
}
1199,10 → 1211,11
dpgw_dpc = tt.get_dpgw_dpc(pc,t);
mg = tt.get_mg(pc,pg,t);
kgc = deff*dpgw_dpc/t/mg;//water vapour
//kgc = deff*dpgw_dpc/t/mg;//water vapour
//kgc = -1.0*kgc;//sign correction
//correction 09/12/2022 by TKr
kgc = -rhog*ma*mw/mg/mg*deff/pg*dpgw_dpc;
 
kgc = -1.0*kgc;//sign correction
 
return(kgc);
}
 
1261,10 → 1274,12
dpgw_dt = tt.get_dpgw_dt(pc,t);
mg = tt.get_mg(pc,pg,t);
 
kgt = deff*dpgw_dt/t/mg;//water vapour
//kgt = deff*dpgw_dt/t/mg;//water vapour
//kgt = -1.0*kgt;//sign correction
//correction by TKr 09/12/2022
kgt = -rhog*ma*mw/mg/mg*deff/pg*dpgw_dt;
 
kgt = -1.0*kgt;//sign correction
 
return(kgt);
}
 
1321,10 → 1336,8
 
lambdaeff = tt.get_lambdaeff(pc,pg,t,ipp);
 
ktt1 = -1.0*lambdaeff;
ktt1 = lambdaeff;
 
ktt1 = -1.0*ktt1;//sign correction
 
return(ktt1);
}
 
1340,7 → 1353,7
*/
double multiph::get_ktt2a(double pc,double pg,double t,long ipp)
{
double cpw,rhow,kintr,krw,muw,ktt2a;
double cpw,rhow,kintr,krw,muw,ktt2a,phi,s;
state_eq tt;
//gas pressure check
1353,8 → 1366,11
kintr = tt.get_kintr(pc,pg,t,ipp);
krw = tt.get_krw(pc,t,ipp);
muw = tt.get_muw(t);
phi = tt.get_phi(t,ipp);
s = tt.get_s(pc,t,ipp);
 
ktt2a = cpw*rhow*kintr*krw/muw;
//correction 09/12/2022
ktt2a = phi*s*cpw*rhow*kintr*krw/muw;
return(ktt2a);
}
1396,7 → 1412,7
*/
double multiph::get_ktt2c(double pc,double pg,double t,long ipp)
{
double rhocpg,kintr,krg,mug,ktt2c;
double rhocpg,kintr,krg,mug,ktt2c,phi,s;
state_eq tt;
//gas pressure check
1408,8 → 1424,10
kintr = tt.get_kintr(pc,pg,t,ipp);
krg = tt.get_krg(pc,t,ipp);
mug = tt.get_mug(pc,pg,t);
phi = tt.get_phi(t,ipp);
s = tt.get_s(pc,t,ipp);
 
ktt2c = rhocpg*kintr*krg/mug;
ktt2c = phi*s*rhocpg*kintr*krg/mug;
 
return(ktt2c);
}
/trunk/SIFEL/TRFEL/SRC/npsolvert.cpp
527,11 → 527,12
if (Mesprt==1) fprintf (stdout,"\n\n --------------------------------------------------------------");
if (Mesprt==1) fprintf (stdout,"\n TRFEL Time step = %ld, Time %e, Time increment = %e",istep,Tp->time,dt);
if (Mesprt==1) fprintf (stdout,"\n --------------------------------------------------------------\n");
 
//fprintf (Outt,"\n\n --------------------------------------------------------------");
//fprintf (Outt,"\n TRFEL Time step = %ld, Time %e, Time increment = %e",istep,Tp->time,dt);
//fprintf (Outt,"\n --------------------------------------------------------------\n");
//fflush(Outt);
 
// number of transport degrees of freedom
n=Ndoft;
/trunk/SIFEL/TRFEL/SRC/threemedia.cpp
274,6 → 274,17
i = Tm->ip[ipp].idm;
switch (Tm->ip[ipp].tm){//material type
case concreteB:
case baroghelB:
case C60baroghelB:
case C30baroghelB:
case o30bazantB:
case C60bazantB:{
//multiph mtph;
 
//mtph.rhs_volume2(d,ri,ci,ipp);//must be completed
break;
}
case consolhawf3:{
Tm->consol_hawf3[i].rhs_volume2(c,ri,ci,ipp);
break;
/trunk/SIFEL/TRFEL/SRC/transmat.cpp
680,6 → 680,17
moisth[i].give_dof_names(dofname, ntm);
break;
}
case baroghelB:{
baroghel[i].give_dof_names(dofname, ntm);
break;
}
case C60bazantB:{
C60bazant[i].give_dof_names(dofname, ntm);
break;
}
 
/*
case kunzel2:
//kun2[i].give_dof_names(dofname, ntm);
4220,8 → 4231,8
break;
}
case moistheat:{
moisth[i].values_correction (nv);
break;
moisth[i].values_correction (nv);
break;
}
case glasgow:{
break;
4266,11 → 4277,11
break;
}
case consolhawf3:{
//consol_hawf3[i].values_correction (nv, ipp);//doplnit??!!
consol_hawf3[i].values_correction (nv, ipp);
break;
}
case consolhwf2:{
//consol_hwf2[i].values_correction (nv, ipp);//doplnit??!!
consol_hwf2[i].values_correction (nv, ipp);
break;
}
default:{
/trunk/SIFEL/TRFEL/SRC/van_genuchten.cpp
89,7 → 89,7
*/
double van_genuchten_reten::dsw_dpw(double pw)
{
double dsw_dpw,hp,expm;
double dsw_dpw,hp;
pw = pw/1000.0;//in kPa