Subversion Repositories sifel

Compare Revisions

Ignore whitespace Rev 772 → Rev 773

/trunk/SIFEL/METR/SRC/consol_wf1c.cpp
0,0 → 1,1149
/*
File: consol_wf1c.cpp
Author: Tomas Krejci, 12/09/2008, revised 23/08/2012
Purpose: material model for saturated-nonsaturated one-phase flow in deforming medium
sources: Lewis and Schrefler's book pp. 86-89, initial material parameters are set for benchmark on page n. 168;
 
unknowns: TRFEL: number of unknowns=1, pw = liqud water pressure; MEFEL: according to problem_type
*/
 
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "constrel.h"
#include "onemediumc.h"
#include "coupmatu.h"
#include "coupmatl.h"
#include "globalc.h"
#include "global.h"
#include "globalt.h"
#include "consol_wf1c.h"
#include "intpoints.h"
 
 
con_wf1matc::con_wf1matc()
{
compress = 0; //compressible grains: 0=no; 1=yes
//
alpha = 1.0;//Biot's constant [-] alpha = 1 - kt/ks
//
rhos = 2000.0;//solid grain density [kg/m^3]
//
phi0 = 0.297;//initial porosity [-]
//
rhow = 1000.0;//water density [kg/m^3]
//
}
 
con_wf1matc::~con_wf1matc()
{}
 
 
/**
function reads parameters
@param in - input file
 
12/9/2008, TKr
*/
void con_wf1matc::read(XFILE *in)
{
xfscanf (in,"%k%m","waterflowmechtype",&waterflowmechtype_kwdset, &model_type); //water flow model type
xfscanf (in,"%d", &compress); //compressibility of grains
switch (model_type){
case lewis_and_schrefler_coup:
case lewis_and_schrefler_mefel_coup: {//Lewis and Schrefler's model
xfscanf (in,"%lf %lf %lf %d", &alpha, &phi0, &rhos, &sr_type);
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 gardner_exponential_sr:{//dependent on saturation degree:
gardner_ret.read(in);
break;
}
case potts_log_linear_sr:{//exponential
potts_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();
}
}
// incompressible grains:
//if(compress == 0)
//alpha = 1.0;
}
 
 
/**
function prints parameters
@param out - output file
 
12/9/2008, TKr
*/
void con_wf1matc::print(FILE *out)
{
/* fprintf(out," %d",model_type);
switch (model_type){
case lewis_and_schrefler_coup:{//Lewis and Schrefler's model
fprintf (out," %le %le %le %le %le %le %le %le\n", alpha, ks, phi0, kw, rhow, muw0, kintr, rhos);
break;
}
case lewis_and_schrefler_mefel_coup:{//Lewis and Schrefler's model and moisture storage functions is included in MEFEL
fprintf (out," %le %le %le %le %le %le %le %le %le\n", mefel_units, alpha, ks, phi0, kw, rhow, muw0, kintr, rhos);
break;
}
case van_genuchten_coup:{//partially saturated medium =Van Genuchten model
fprintf (out,"\n %le %le %le %le %le %le %le %le %le \n", mefel_units, alpha, ks, phi0, kw, rhow, muw0, kintr, rhos);
van_genuchten_ret.print(out);
break;
}
case consol_camclay_coup://fully saturated isotropic consolidation with changing porosity from cam-clay model
case consol_camclay_mech_coup:{
fprintf (out," %le %le %le %le %le %le %le %le %le\n", mefel_units, alpha, ks, phi0, kw, rhow, muw0, kintr, rhos);
break;
}
default:{
fprintf (stderr,"\n unknown model type is required");
fprintf (stderr,"\n in function (%s, line %d).\n",__FILE__,__LINE__);
}
}
*/
}
 
 
/**
function computes conductivity matrix of the material
in the required integration point for upper block
@param d - conductivity %matrix of material
@param ri - row index
@param ci - column index
@param ipp - number of integration point
 
01/05/2018, TKr
*/
void con_wf1matc::matcond_u (matrix &d,long ri,long ci,long ipp)
{
long m;
m = d.m;
 
switch (m){
case 1:{
matcond1d_u (d,ri,ci,ipp);//1D
break;
}
case 3:{
matcond2d_u (d,ri,ci,ipp);//2D
break;
}
case 4:{
matcond2d_ax_u (d,ri,ci,ipp);//2D - axisymmetric
break;
}
case 6:{
matcond3d_u (d,ri,ci,ipp);//3D
break;
}
default:{
fprintf (stderr,"\n unknown number of components of stress tensor is required");
fprintf (stderr,"\n in function (%s, line %d).\n",__FILE__,__LINE__);
}
}
}
 
 
/**
function creates conductivity matrix of the material for 1D problems
for upper block
 
@param d - conductivity %matrix of the material
@param ri - row index
@param ci - column index
@param ipp - number of integration point
 
01/05/2018, TKr
*/
void con_wf1matc::matcond1d_u (matrix &d,long ri,long ci,long ipp)
{
double k;
double pw;
k = 0.0;
pw = Cmu->ip[ipp].av[0];
if((ri == 0) && (ci == 0))
k = get_kuw(pw,ipp);
d[0][0] = k;
}
 
/**
function creates conductivity matrix of the material for 2D problems
for upper block
 
@param d - conductivity %matrix of the material
@param ri - row index
@param ci - column index
@param ipp - number of integration point
 
01/05/2018, TKr
*/
void con_wf1matc::matcond2d_u (matrix &d,long ri,long ci,long ipp)
{
double k;
double pw;
k = 0.0;
pw = Cmu->ip[ipp].av[0];
if((ri == 0) && (ci == 0))
k = get_kuw(pw,ipp);
fillm(0.0,d);
 
d[0][0] = k;
d[1][0] = k;
d[2][0] = 0.0;
}
 
 
/**
function creates conductivity matrix of the material for 2D axisymmetric problems
for upper block
 
@param d - conductivity %matrix of the material
@param ri - row index
@param ci - column index
@param ipp - number of integration point
 
10/05/2018, TKr
*/
void con_wf1matc::matcond2d_ax_u (matrix &d,long ri,long ci,long ipp)
{
double k;
double pw;
k = 0.0;
pw = Cmu->ip[ipp].av[0];
if((ri == 0) && (ci == 0))
k = get_kuw(pw,ipp);
fillm(0.0,d);
 
d[0][0] = k;
d[1][0] = k;
d[2][0] = k;
d[3][0] = 0.0;
}
 
/**
function creates conductivity matrix of the material for 3D problems
for upper block
 
@param d - conductivity %matrix of the material
@param ri - row index
@param ci - column index
@param ipp - number of integration point
 
01/05/2018, TKr
*/
void con_wf1matc::matcond3d_u (matrix &d,long ri,long ci,long ipp)
{
double k;
double pw;
k = 0.0;
pw = Cmu->ip[ipp].av[0];
if((ri == 0) && (ci == 0))
k = get_kuw(pw,ipp);
fillm(0.0,d);
d[0][0]=k;
d[1][0]=k;
d[2][0]=k;
d[3][0]=0.0;
d[4][0]=0.0;
d[5][0]=0.0;
}
 
 
/**
function computes capacity matrix of the material
in the required integration point for upper block
@param c - capacity matrix of the material
@param ri - row index
@param ci - column index
@param ipp - number of integration point
01/05/2018, TKr
*/
void con_wf1matc::matcap_u (matrix &d,long ri,long ci,long ipp)
{
long m;
 
m = d.m;
 
switch (m){
case 1:{
matcap1d_u (d,ri,ci,ipp);//1D
break;
}
case 3:{
matcap2d_u (d,ri,ci,ipp);//2D
break;
}
case 4:{
matcap2d_ax_u (d,ri,ci,ipp);//2D - axisymmetric
break;
}
case 6:{
matcap3d_u (d,ri,ci,ipp);//3D
break;
}
default:{
fprintf (stderr,"\n unknown number of components of stress tensor is required");
fprintf (stderr,"\n in function (%s, line %d).\n",__FILE__,__LINE__);
}
}
}
 
 
/**
function creates capacity matrix of the material for 1D problems
for upper block
 
@param c - capacity %matrix of the material
@param ri - row index
@param ci - column index
@param ipp - number of integration point
 
01/05/2018, TKr
*/
void con_wf1matc::matcap1d_u (matrix &d,long ri,long ci,long ipp)
{
double pw,c;
pw = Cmu->ip[ipp].av[0];
if((ri == 0) && (ci == 0))
c = get_capuw(pw,ipp);
d[0][0] = c;
}
 
 
/**
function creates capacity matrix of the material for 2D problems
for upper block
 
@param c - capacity %matrix of the material
@param ri - row index
@param ci - column index
@param ipp - number of integration point
 
01/05/2018, TKr
*/
void con_wf1matc::matcap2d_u (matrix &d,long ri,long ci,long ipp)
{
double pw,c;
pw = Cmu->ip[ipp].av[0];
if((ri == 0) && (ci == 0))
c = get_capuw(pw,ipp);
fillm(0.0,d);
 
d[0][0] = c;
d[1][0] = c;
d[2][0] = 0.0;
}
 
 
/**
function creates capacity matrix of the material for 2D axisymmetric problems
for upper block
 
@param c - capacity %matrix of the material
@param ri - row index
@param ci - column index
@param ipp - number of integration point
 
10/05/2018, TKr
*/
void con_wf1matc::matcap2d_ax_u (matrix &d,long ri,long ci,long ipp)
{
double pw,c;
c = 0.0;
pw = Cmu->ip[ipp].av[0];
if((ri == 0) && (ci == 0))
c = get_capuw(pw,ipp);
 
fillm(0.0,d);
 
d[0][0] = c;
d[1][0] = c;
d[2][0] = c;
d[3][0] = 0.0;
}
 
 
/**
function creates capacity matrix of the material for 3D problems
for upper block
 
@param c - capacity %matrix of the material
@param ri - row index
@param ci - column index
@param ipp - number of integration point
 
01/05/2018, TKr
*/
void con_wf1matc::matcap3d_u (matrix &d,long ri,long ci,long ipp)
{
double pw,c;
pw = Cmu->ip[ipp].av[0];
if((ri == 0) && (ci == 0))
c = get_capuw(pw,ipp);
fillm(0.0,d);
 
d[0][0] = c;
d[1][0] = c;
d[2][0] = c;
d[3][0] = 0.0;
d[4][0] = 0.0;
d[5][0] = 0.0;
}
 
 
/**
function computes conductivity matrix of the material
in the required integration point for lower block
@param d - conductivity %matrix of material
@param ri - row index
@param ci - column index
@param ipp - number of integration point
 
01/05/2018, TKr
*/
void con_wf1matc::matcond_l (matrix &d,long ri,long ci,long ipp)
{
long m;
m = d.n;
 
switch (m){
case 1:{
matcond1d_l (d,ri,ci,ipp);//1D
break;
}
case 3:{
matcond2d_l (d,ri,ci,ipp);//2D
break;
}
case 4:{
matcond2d_ax_l (d,ri,ci,ipp);//2D - axisymmetric
break;
}
case 6:{
matcond3d_l (d,ri,ci,ipp);//3D
break;
}
default:{
fprintf (stderr,"\n unknown number of components of stress tensor is required");
fprintf (stderr,"\n in function (%s, line %d).\n",__FILE__,__LINE__);
}
}
}
 
 
/**
function creates conductivity matrix of the material for 1D problems
for lower block
 
@param d - conductivity %matrix of the material
@param ri - row index
@param ci - column index
@param ipp - number of integration point
 
01/05/2018, TKr
*/
void con_wf1matc::matcond1d_l (matrix &d,long ri,long ci,long ipp)
{
double k;
double pw;
k = 0.0;
pw = Cml->ip[ipp].av[0];
if((ri == 0) && (ci == 0))
k = get_kwu(pw,ipp);
d[0][0] = k;
}
 
/**
function creates conductivity matrix of the material for 2D problems
for lower block
 
@param d - conductivity %matrix of the material
@param ri - row index
@param ci - column index
@param ipp - number of integration point
 
01/05/2018, TKr
*/
void con_wf1matc::matcond2d_l (matrix &d,long ri,long ci,long ipp)
{
double k;
double pw;
k = 0.0;
pw = Cml->ip[ipp].av[0];
 
if((ri == 0) && (ci == 0))
k = get_kwu(pw,ipp);
fillm(0.0,d);
d[0][0] = k;
d[0][1] = k;
d[0][2] = 0.0;
}
 
/**
function creates conductivity matrix of the material for 2D axisymmetric problems
for lower block
 
@param d - conductivity %matrix of the material
@param ri - row index
@param ci - column index
@param ipp - number of integration point
 
10/05/2018, TKr
*/
void con_wf1matc::matcond2d_ax_l (matrix &d,long ri,long ci,long ipp)
{
double k;
double pw;
k = 0.0;
pw = Cml->ip[ipp].av[0];
 
if((ri == 0) && (ci == 0))
k = get_kwu(pw,ipp);
fillm(0.0,d);
d[0][0] = k;
d[0][1] = k;
d[0][2] = k;
d[0][3] = 0.0;
}
 
/**
function creates conductivity matrix of the material for 3D problems
for lower block
 
@param d - conductivity %matrix of the material
@param ri - row index
@param ci - column index
@param ipp - number of integration point
 
01/05/2018, TKr
*/
void con_wf1matc::matcond3d_l (matrix &d,long ri,long ci,long ipp)
{
double k;
double pw;
k = 0.0;
pw = Cml->ip[ipp].av[0];
if((ri == 0) && (ci == 0))
k = get_kwu(pw,ipp);
fillm(0.0,d);
d[0][0]=k;
d[0][1]=k;
d[0][2]=k;
d[0][3]=0.0;
d[0][4]=0.0;
d[0][5]=0.0;
}
 
 
/**
function computes capacity matrix of the material
in the required integration point for lower block
@param c - capacity matrix of the material
@param ri - row index
@param ci - column index
@param ipp - number of integration point
01/05/2018, TKr
*/
void con_wf1matc::matcap_l (matrix &d,long ri,long ci,long ipp)
{
long m;
 
m = d.n;
 
switch (m){
case 1:{
matcap1d_l (d,ri,ci,ipp);//1D
break;
}
case 3:{
matcap2d_l (d,ri,ci,ipp);//2D
break;
}
case 4:{
matcap2d_ax_l (d,ri,ci,ipp);//2D - axisymmetric
break;
}
case 6:{
matcap3d_l (d,ri,ci,ipp);//3D
break;
}
default:{
fprintf (stderr,"\n unknown number of components of stress tensor is required");
fprintf (stderr,"\n in function (%s, line %d).\n",__FILE__,__LINE__);
}
}
}
 
 
/**
function creates capacity matrix of the material for 1D problems
for lower block
 
@param c - capacity %matrix of the material
@param ri - row index
@param ci - column index
@param ipp - number of integration point
 
01/05/2018, TKr
*/
void con_wf1matc::matcap1d_l (matrix &d,long ri,long ci,long ipp)
{
double pw,c;
c = 0.0;
pw = Cml->ip[ipp].av[0];
if((ri == 0) && (ci == 0))
c = get_capwu(pw,ipp);
d[0][0] = c;
}
 
 
/**
function creates capacity matrix of the material for 2D problems
for lower block
 
@param c - capacity %matrix of the material
@param ri - row index
@param ci - column index
@param ipp - number of integration point
 
01/05/2018, TKr
*/
void con_wf1matc::matcap2d_l (matrix &d,long ri,long ci,long ipp)
{
double pw,c;
c = 0.0;
pw = Cml->ip[ipp].av[0];
if((ri == 0) && (ci == 0))
c = get_capwu(pw,ipp);
fillm(0.0,d);
 
d[0][0] = c;
d[0][1] = c;
d[0][2] = 0.0;
}
 
 
/**
function creates capacity matrix of the material for 2D axisymmetric problems
for lower block
 
@param c - capacity %matrix of the material
@param ri - row index
@param ci - column index
@param ipp - number of integration point
 
10/05/2018, TKr
*/
void con_wf1matc::matcap2d_ax_l (matrix &d,long ri,long ci,long ipp)
{
double pw,c;
c = 0.0;
pw = Cml->ip[ipp].av[0];
if((ri == 0) && (ci == 0))
c = get_capwu(pw,ipp);
fillm(0.0,d);
 
d[0][0] = c;
d[0][1] = c;
d[0][2] = c;
d[0][3] = 0.0;
}
 
 
/**
function creates capacity matrix of the material for 3D problems
for lower block
 
@param c - capacity %matrix of the material
@param ri - row index
@param ci - column index
@param ipp - number of integration point
 
01/05/2018, TKr
*/
void con_wf1matc::matcap3d_l (matrix &d,long ri,long ci,long ipp)
{
double pw,c;
c = 0.0;
pw = Cml->ip[ipp].av[0];
if((ri == 0) && (ci == 0))
c = get_capwu(pw,ipp);
fillm(0.0,d);
d[0][0]=c;
d[0][1]=c;
d[0][2]=c;
d[0][3]=0.0;
d[0][4]=0.0;
d[0][5]=0.0;
}
 
 
 
/**
function computes volume part of right-hand side matrix
in the required integration point
@param d - right-hand side %matrix of material
@param ri - row index
@param ci - column index
@param ipp - number of integration point
*/
void con_wf1matc::rhs_u1 (matrix &d,long ri,long ci,long ipp)
{
long m;
m = d.m;
switch (m){
case 1:{
rhs1d1 (d,ri,ci,ipp);//1D
break;
}
case 3:{
rhs2d1 (d,ri,ci,ipp);//2D
break;
}
case 4:{
rhs2d_ax_1 (d,ri,ci,ipp);//2D - axisymmetric
break;
}
case 6:{
rhs3d1 (d,ri,ci,ipp);//3D
break;
}
default:{
fprintf (stderr,"\n unknown number of components of stress tensor is required");
fprintf (stderr,"\n in function (%s, line %d).\n",__FILE__,__LINE__);
}
}
}
 
 
 
/**
function creates volume right-hand side matrix of the material for 1D problems
@param d - right-hand %matrix of the material
@param ri - row index
@param ci - column index
@param ipp - number of integration point
*/
void con_wf1matc::rhs1d1 (matrix &d,long ri,long ci,long ipp)
{
double f,pw;
pw = Cmu->ip[ipp].av[0];
f = get_fu1(pw,ipp);
fillm(0.0,d);
d[0][0] = f*Cp->gr[0];
}
 
/**
function creates volume right-hand side matrix of the material for 2D problems
@param d - right-hand %matrix of the material
@param ri - row index
@param ci - column index
@param ipp - number of integration point
*/
void con_wf1matc::rhs2d1 (matrix &d,long ri,long ci,long ipp)
{
double f,pw;
pw = Cmu->ip[ipp].av[0];
f = get_fu1(pw,ipp);
fillm(0.0,d);
d[0][0] = f*Cp->gr[0];
d[1][0] = f*Cp->gr[1];
}
 
/**
function creates volume right-hand side matrix of the material for 2D axisymmetric problems
@param d - right-hand %matrix of the material
@param ri - row index
@param ci - column index
@param ipp - number of integration point
*/
void con_wf1matc::rhs2d_ax_1 (matrix &d,long ri,long ci,long ipp)
{
double f,pw;
pw = Cmu->ip[ipp].av[0];
f = get_fu1(pw,ipp);
fillm(0.0,d);
d[0][0] = f*Cp->gr[0];
d[1][0] = f*Cp->gr[1];
}
 
/**
function creates volume right-hand side matrix of the material for 3D problems
@param d - right-hand %matrix of the material
@param ri - row index
@param ci - column index
@param ipp - number of integration point
*/
void con_wf1matc::rhs3d1 (matrix &d,long ri,long ci,long ipp)
{
double f,pw;
pw = Cmu->ip[ipp].av[0];
 
f = get_fu1(pw,ipp);
fillm(0.0,d);
d[0][0] = f*Cp->gr[0];
d[1][0] = f*Cp->gr[1];
d[2][0] = f*Cp->gr[2];
}
 
 
 
/**
function computes degree of saturation(sorption curve)
@param pw - water pressure
 
@retval sw - degree of saturation
 
12/9/2008, TKr
*/
double con_wf1matc::get_sw(double pw, long ipp)
{
double sw;
switch (sr_type){
case lewis_and_schrefler_sr:{//Lewis and Schrefler's book
sw = lewis_ret.sw(pw);
break;
}
case mefel_sr:{//saturation degree ind its derivative are obtained from mefel;
sw = Tm->givenontransq(saturation_deg, ipp); //actual saturation degree
break;
}
case table_sr:{//saturation degree and its derivative are obtained from table;
//actual saturation degree
if (data.tfunc == tab)
{
if ((pw < data.tabf->x[0]) || (pw > data.tabf->x[data.tabf->asize-1]))
{
print_err("required value %le is out of table range <%le;%le> on ip=%ld\n",
__FILE__, __LINE__, __func__, pw, data.tabf->x[0], data.tabf->x[data.tabf->asize-1], ipp);
abort();
}
}
sw = data.getval (pw);
break;
}
 
case gardner_exponential_sr:{//partially saturated medium = Exponential model, Gardner
sw = gardner_ret.sw(pw);
break;
}
case potts_log_linear_sr:{//partially saturated medium = Log linear model, Potts
sw = potts_ret.sw(pw);
break;
}
case van_genuchten_sr:{//partially saturated medium = Van Genuchten model
sw = van_genuchten_ret.sw(pw);
break;
}
case van_genuchten2_sr:{//partially saturated medium = Van Genuchten model
//sw = van_genuchten_ret.sw2(pw);
break;
}
default:{
print_err ("unknown model type is required", __FILE__, __LINE__, __func__);
abort();
}
}
return(sw);
}
 
 
/**
function returns porosity
 
@retval phi - porosity
 
12/9/2008, TKr
*/
double con_wf1matc::get_phi(double pw,long ipp)
{
double phi;
 
switch (model_type){
case lewis_and_schrefler_coup:{//Lewis and Schrefler's model
phi = phi0;
break;
}
case lewis_and_schrefler_mefel_coup:{//Lewis and Schrefler's model and moisture storage functions is included in MEFEL
phi = Tm->givenontransq(porosity, ipp); //actual porosity from models included in TRFEL
break;
}
default:{
fprintf (stderr,"\n unknown model type is required");
fprintf (stderr,"\n in function (%s, line %d).\n",__FILE__,__LINE__);
}
}
return(phi);
}
 
 
 
/**
function returns volume density of soil skeleton
 
@retval rhos - volume density of concrete skeleton
 
12/9/2008, TKr
*/
double con_wf1matc::get_rhos(double pw,long ipp)
{
return(rhos);
}
 
 
/**
function creates conductivity coefficient of the general material
@param pw - water pressure
 
@retval kuw - conductivity coefficient
 
12/9/2008, TKr
*/
double con_wf1matc::get_kuw(double pw,long ipp)
{
double s,kuw;
switch (model_type){
case lewis_and_schrefler_coup:{//Lewis and Schrefler's model
s = get_sw(pw,ipp);
kuw = -alpha*s;//p. 88
break;
}
case lewis_and_schrefler_mefel_coup:{//Lewis and Schrefler's model and moisture storage functions is included in MEFEL
kuw = 0.0;// this part is in MEFEL already included
break;
}
default:{
fprintf (stderr,"\n unknown model type is required");
fprintf (stderr,"\n in function (%s, line %d).\n",__FILE__,__LINE__);
}
}
 
return(kuw);
 
}
 
 
/**
function creates conductivity coefficient of the general material
@param pw - water pressure
 
@retval kwu - conductivity coefficient
 
12/9/2008, TKr
*/
double con_wf1matc::get_kwu(double pw,long ipp)
{
double kwu;
 
kwu = 0.0;
 
return(kwu);
 
}
 
 
/**
function creates capacity coefficient of the general material
@param pw - water pressure
 
@retval capuw - capacity coefficient
 
12/9/2008, TKr
*/
double con_wf1matc::get_capuw(double pw,long ipp)
{
double capuw;
capuw = 0.0;
 
return(capuw);
 
}
 
 
/**
function creates capacity coefficient of the general material
@param pw - water pressure
 
@retval capwu - capacity coefficient
 
12/9/2008, TKr
*/
double con_wf1matc::get_capwu(double pw,long ipp)
{
double sw,capwu,n,dsr_depsv;
capwu = 0.0;
 
switch (model_type){
case lewis_and_schrefler_coup:{//Lewis and Schrefler's model
sw = get_sw(pw,ipp);
capwu = alpha*sw;//p. 88
break;
}
case lewis_and_schrefler_mefel_coup:{//Lewis and Schrefler's model and moisture storage functions is included in MEFEL
capwu = 0.0;// this part is in TRFEL already included
break;
}
default:{
fprintf (stderr,"\n unknown model type is required");
fprintf (stderr,"\n in function (%s, line %d).\n",__FILE__,__LINE__);
}
}
return(capwu);
 
}
 
 
/**
function returns coefficient for righ-hand side of the general material
@param pw - water pressure
@param pg - air pressure
 
@retval fu1 - first part for right-hand side for balance equation
 
08/09/2021, TKr
*/
double con_wf1matc::get_fu1(double pw,long ipp)
{
double sw,n,fu1;
switch (model_type){
case lewis_and_schrefler_coup:{//Lewis and Schrefler's model
sw = get_sw(pw,ipp);
n = get_phi(pw,ipp);
fu1 = rhos*(1-n) + sw*n*rhow;//p. 88
break;
}
case lewis_and_schrefler_mefel_coup:{//Lewis and Schrefler's model and moisture storage functions is included in MEFEL
sw = get_sw(pw,ipp);
n = get_phi(pw,ipp);
fu1 = rhos*(1-n) + sw*n*rhow;//p. 88
fu1 = 0.0;// temporarily no gravity acceleration effect is assumed
break;
}
default:{
fprintf (stderr,"\n unknown model type is required");
fprintf (stderr,"\n in function (%s, line %d).\n",__FILE__,__LINE__);
}
}
return(fu1);
}
/trunk/SIFEL/METR/SRC/consol_wf1c.h
0,0 → 1,84
#ifndef CON_WF1MATC_H
#define CON_WF1MATC_H
 
#include "genfile.h"
#include "lewis_schrefler.h"
#include "gardner.h"
#include "potts.h"
#include "van_genuchten.h"
 
/**
This class defines model for water flow in soils (deformable porous medium)
*/
 
class con_wf1matc
{
public:
con_wf1matc(); //constructor
~con_wf1matc(); //destructor
 
void read(XFILE *in);
void print(FILE *out);
 
void matcond_u (matrix &d,long ri,long ci,long ipp);
void matcap_u (matrix &d,long ri,long ci,long ipp);
void matcond1d_u (matrix &d,long ri,long ci,long ipp);
void matcond2d_u (matrix &d,long ri,long ci,long ipp);
void matcond2d_ax_u (matrix &d,long ri,long ci,long ipp);
void matcond3d_u (matrix &d,long ri,long ci,long ipp);
void matcap1d_u (matrix &d,long ri,long ci,long ipp);
void matcap2d_u (matrix &d,long ri,long ci,long ipp);
void matcap2d_ax_u (matrix &d,long ri,long ci,long ipp);
void matcap3d_u (matrix &d,long ri,long ci,long ipp);
 
void matcond_l (matrix &d,long ri,long ci,long ipp);
void matcap_l (matrix &d,long ri,long ci,long ipp);
void matcond1d_l (matrix &d,long ri,long ci,long ipp);
void matcond2d_l (matrix &d,long ri,long ci,long ipp);
void matcond2d_ax_l (matrix &d,long ri,long ci,long ipp);
void matcond3d_l (matrix &d,long ri,long ci,long ipp);
void matcap1d_l (matrix &d,long ri,long ci,long ipp);
void matcap2d_l (matrix &d,long ri,long ci,long ipp);
void matcap2d_ax_l (matrix &d,long ri,long ci,long ipp);
void matcap3d_l (matrix &d,long ri,long ci,long ipp);
void rhs_u1 (matrix &d,long ri,long ci,long ipp);
void rhs1d1 (matrix &d,long ri,long ci,long ipp);
void rhs2d1 (matrix &d,long ri,long ci,long ipp);
void rhs2d_ax_1 (matrix &d,long ri,long ci,long ipp);
void rhs3d1 (matrix &d,long ri,long ci,long ipp);
 
double get_sw(double pw,long ipp);
double get_phi(double pw,long ipp);
double get_rhos(double pw,long ipp);
 
double get_kuw(double pw,long ipp);
double get_kwu(double pw,long ipp);
 
double get_capuw(double pw,long ipp);
double get_capwu(double pw,long ipp);
 
double get_fu1(double pw,long ipp);
 
double nu;
 
/// Lewis and Schrefler's retention curve:
lewis_reten lewis_ret;
/// Gardner's retention curve:
gardner_reten gardner_ret;
/// Potts' retention curve:
potts_reten potts_ret;
/// van Genuchten's retention curve:
van_genuchten_reten van_genuchten_ret;
/// general function for retention curve given by set of data
gfunct data;
private:
waterflowmechtype model_type;
int compress,sr_type;
double alpha,phi0,rhos,rhow;
};
 
#endif