| /trunk/SIFEL/GEFEL/galias.h |
|---|
| 240,13 → 240,13 |
| enum nontransquant {precons_press=1, mean_stress_eff=2, virgin_porosity=3, init_porosity=4, |
| scal_iso_damage=5,proc_zone_length=6,crack_width=7,saturation_deg=8,der_saturation_deg=9, |
| porosity=10,void_ratio=11,advect_vel_x=12,advect_vel_y=13,advect_vel_z=14,der_saturation_deg_dtemp=15, |
| strain_vol_rate=16,der_saturation_deg_depsv=17, bulk_modulus=18}; |
| strain_vol_rate=16,der_saturation_deg_depsv=17, bulk_modulus=18, mmean_stress=19}; |
| const enumstr nontransquantstr[] = {{"precons_press",1}, {"mean_stress_eff",2}, {"virgin_porosity",3}, |
| {"init_porosity",4}, {"scal_iso_damage",5}, {"proc_zone_length",6}, |
| {"crack_width",7}, {"saturation_deg",8}, {"der_saturation_deg",9}, |
| {"porosity",10}, {"void_ratio",11}, {"advect_vel_x",12}, {"advect_vel_y",13}, |
| {"advect_vel_z",14}, {"der_saturation_deg_dtemp",15}, {"strain_vol_rate",16}, |
| {"der_saturation_deg_depsv",17}, {"bulk_modulus",18}}; |
| {"der_saturation_deg_depsv",17}, {"bulk_modulus",18}, {"mmean_stress",19}}; |
| const kwdset nontransquant_kwdset(sizeof(nontransquantstr)/sizeof(*nontransquantstr), nontransquantstr); |
| /trunk/SIFEL/GEFEL/timecontr.cpp |
|---|
| 350,26 → 350,24 |
| if (tct == adaptivemax) // minimum time increment dtmin is given by timefun |
| dtmin = forwarddtu; |
| /* |
| if (adaptive_minmax) |
| { |
| dtmin = dtminfun->getval(time); |
| dtmax = dtmaxfun->getval(time); |
| } |
| */ |
| /* |
| if(adaptive || adaptive_minmax) |
| { |
| if ((bfdtu != forwarddtu) && (dt >= bsdt)) |
| if (tct == adaptive_minmax) |
| { |
| // user has prescribed change in time increment just for the given time step |
| // and solver did not require shorter time increment in the performed time step |
| forwarddt = forwarddtu; |
| bfdtu = forwarddtu; |
| dtmin = dtminfun->getval(time); |
| dtmax = dtmaxfun->getval(time); |
| } |
| } |
| */ |
| //if(adaptive || adaptive_minmax) |
| if(tct == adaptive_minmax) |
| { |
| if ((bfdtu != forwarddtu) && (dt >= bsdt)) |
| { |
| // user has prescribed change in time increment just for the given time step |
| // and solver did not require shorter time increment in the performed time step |
| forwarddt = forwarddtu; |
| bfdtu = forwarddtu; |
| } |
| } |
| } |
| if (forwarddt>dtmax)//maximum time increment |
| /trunk/SIFEL/MEFEL/PREP/hngen.cpp |
|---|
| File deleted |
| /trunk/SIFEL/MEFEL/PREP/hngen.h |
|---|
| File deleted |
| /trunk/SIFEL/MEFEL/SRC/bbm.cpp |
|---|
| 94,7 → 94,7 |
| double v_ini, v_pc0, v_lambda1; |
| double i1s, j2s, s; |
| double v_kappa1, p1, bar_pc_0, pc_0; |
| long i, ncompstr= Mm->ip[ipp].ncompstr;; |
| long i, ncompstr= Mm->ip[ipp].ncompstr; |
| vector sig(ASTCKVEC(ncompstr)),sigt(ASTCKVEC(6)), q(ASTCKVEC(2)); |
| double err=sra.give_err (); |
| // variables used for the calculation of initial value of preconsolidation pressure for the fully saturated state |
| /trunk/SIFEL/MEFEL/SRC/cpsolver.cpp |
|---|
| 182,9 → 182,9 |
| // update status indicators of prescribed displacements on elements |
| Mt->elemprescdisp(); |
| // update actual material index for all elements |
| Mm->update_actual_mat_id(lcid, 0); |
| Mm->update_actual_mat_id(lcid, 0, true); |
| // initiation of mechanical material models |
| Mm->initmaterialmodels(lcid); |
| Mm->initmaterialmodels(lcid, true); |
| rest_calc = 1; |
| print_init(-1, "at"); |
| 233,7 → 233,7 |
| // update actual material index for all old elements and call initval for those whose indeces have changed |
| if (rest_calc == 0) |
| { |
| um = Mm->update_actual_mat_id(lcid, 1); |
| um = Mm->update_actual_mat_id(lcid, 1, false); |
| if (um) fusm = yes; |
| else fusm = no; |
| } |
| 286,10 → 286,10 |
| // update actual material index for all new elements |
| // but do not initialize new active material models |
| Mm->update_actual_mat_id(lcid, 0); |
| Mm->update_actual_mat_id(lcid, 0, false); |
| // initiation of mechanical material models |
| Mm->initmaterialmodels(lcid); |
| Mm->initmaterialmodels(lcid, false); |
| // |
| /trunk/SIFEL/MEFEL/SRC/creepdam.cpp |
|---|
| 246,14 → 246,15 |
| @param ipp - integration point pointer |
| @param im - material index |
| @param ido - index of internal variables for given material in the ipp other array |
| @param[in] rinit - flag for initialization after restorage from hdbackup |
| 7.6.2005, TKo |
| */ |
| void creepdam::initvalues (long lcid, long ipp, long im, long ido) |
| void creepdam::initvalues (long lcid, long ipp, long im, long ido, bool rinit) |
| { |
| long ncompo = Mm->givencompeqother(ipp, im+1); |
| Mm->initvalues(lcid,ipp,im+1,ido); |
| Mm->initvalues(lcid,ipp,im+2,ido+ncompo); |
| Mm->initvalues(lcid,ipp,im+1,ido, rinit); |
| Mm->initvalues(lcid,ipp,im+2,ido+ncompo, rinit); |
| } |
| /trunk/SIFEL/MEFEL/SRC/creepdam.h |
|---|
| 31,7 → 31,7 |
| void nlstresses (long ipp, long im, long ido); |
| void nonloc_nlstresses (long ipp, long im, long ido); |
| void updateval (long ipp, long im, long ido); |
| void initvalues (long lcid, long ipp, long im, long ido); |
| void initvalues (long lcid, long ipp, long im, long ido, bool rinit); |
| double give_actual_ft (long ipp, long im, long ido); |
| double give_actual_fc (long ipp, long im, long ido); |
| }; |
| /trunk/SIFEL/MEFEL/SRC/effstress.cpp |
|---|
| 214,16 → 214,17 |
| Additionally, it sets flag probdesc::eigstrcomp for computation of nodal force |
| contribution due to pore pressure. |
| @param lcid - load case id |
| @param ipp - integration point pointer |
| @param im - index of the material in the ipp tm and idm arrays. The standard value is zero. |
| @param ido - index of internal variables in the ip's ipp eqother array. The standard value is zero. |
| @param[in] lcid - load case id |
| @param[in] ipp - integration point pointer |
| @param[in] im - index of the material in the ipp tm and idm arrays. The standard value is zero. |
| @param[in] ido - index of internal variables in the ip's ipp eqother array. The standard value is zero. |
| @param[in] rinit - flag for initialization after restorage from hdbackup |
| @return The function does not return anything but it changes Mm->nonmechq array. |
| Created by Tomas Koudelka, 14.10.2013 |
| */ |
| void effstress::initvalues(long lcid, long ipp, long im, long ido) |
| void effstress::initvalues(long lcid, long ipp, long im, long ido, bool rinit) |
| { |
| long i,ncomp = Mm->ip[ipp].ncompstr; |
| vector uw(ASTCKVEC(ncomp)); |
| 266,7 → 267,7 |
| // initialization of other material models from the model sequence -> stress array on integration points |
| // may be changed |
| Mm->initvalues(lcid, ipp, im+1, ido+1);//eqother values for other materials starts at 2 index |
| Mm->initvalues(lcid, ipp, im+1, ido+1, rinit);//eqother values for other materials starts at 2 index |
| //initial value of stresses |
| for (i=0; i<ncomp; i++) |
| /trunk/SIFEL/MEFEL/SRC/effstress.h |
|---|
| 34,7 → 34,7 |
| void nlstresses (long ipp, long im, long ido); |
| /// initializes mechmat::nonmechq array by constant value of pore pressure |
| void initvalues(long lcid, long ipp, long im, long ido); |
| void initvalues(long lcid, long ipp, long im, long ido, bool rinit); |
| /// marks required non-mechanical quantities |
| void give_reqnmq(long *anmq); |
| /trunk/SIFEL/MEFEL/SRC/elemswitch.cpp |
|---|
| 5802,9 → 5802,9 |
| eid = Mm->elip[app]; |
| if (Gtm->leso[eid] == 1){ |
| Mm->computenlstresses(app,Mm->ip[app]); |
| err = check_math_errel(Mm->elip[app]); |
| if (err) |
| abort(); |
| //err = check_math_errel(Mm->elip[app]); |
| //if (err) |
| //abort();/debug??!! |
| } |
| } |
| /trunk/SIFEL/MEFEL/SRC/fdsolver.cpp |
|---|
| 26,7 → 26,8 |
| for (i=0;i<Lsrs->nlc;i++){ |
| switch (Mp->tforvib){ |
| case newmark:{ |
| newmark_method (i); |
| //newmark_method (i); |
| optim_newmark_method (i); |
| break; |
| } |
| case findiff:{ |
| 196,6 → 197,167 |
| /** |
| Newmark method for solution of forced dynamic problems |
| the method is implicit |
| @param lcid - load case id |
| JK, 11. 4. 2024 |
| */ |
| void optim_newmark_method (long lcid) |
| { |
| long i,j,n; |
| double alpha,delta,time,dt,pdt,end_time,zero=1.0e-10; |
| double *a,*v,*p,*dd,*vv,*lhs,*rhs,*fs; |
| // number of degrees of freedom |
| n = Ndofm; |
| // coefficient of the method |
| alpha=Mp->alphafvn; |
| // coefficient of the method |
| delta=Mp->deltafvn; |
| // initial time |
| time=Mp->timecon.starttime (); |
| // end time |
| end_time = Mp->timecon.endtime (); |
| // first time increment |
| dt = Mp->timecon.initialtimeincr (); |
| // nodal values - displacements |
| // (assigment of initial values) |
| lhs = Lsrs->give_lhs (lcid); |
| // time derivatives of nodal values - velocities |
| // (assigment of initial values) |
| v = Lsrs->give_tdlhs (lcid); |
| // vector of dynamic load |
| rhs = Lsrs->give_rhs (2*lcid); |
| // vector of static load |
| fs = Lsrs->give_rhs (2*lcid+1); |
| // vector of acceleration |
| a = Lsrs->give_stdlhs (lcid); |
| nullv (a,n); |
| // predictor of displacements |
| dd = new double [n]; |
| // predictor of velocities |
| vv = new double [n]; |
| // auxiliary vector |
| p = new double [n]; |
| // assembling of static loads |
| nullv (fs+lcid*n, n); |
| Mb->lc[lcid].assemble (lcid,fs,NULL,1.0); |
| i=0; |
| print_init(-1, "wt"); |
| print_step(lcid, i, Mp->time, rhs); |
| print_flush(); |
| // time initialization |
| Mp->time = Mp->timecon.newtime (); |
| // time increment |
| dt = Mp->timecon.actualbacktimeincr (); |
| pdt = dt; |
| // stiffness matrix |
| stiffness_matrix (lcid); |
| // mass matrix |
| mass_matrix (lcid); |
| // damping matrix |
| damping_matrix (lcid); |
| // loop over time interval |
| do{ |
| if (Mespr==1) fprintf (stdout,"\n time %f",Mp->time); |
| if (fabs(pdt-dt)>zero){ |
| // the time step was changed |
| // therefore the matrices have to be recalculated |
| pdt = dt; |
| // stiffness matrix |
| stiffness_matrix (lcid); |
| // mass matrix |
| mass_matrix (lcid); |
| // damping matrix |
| damping_matrix (lcid); |
| // dynamic stiffness matrix |
| // (stored in array for stiffness matrix) |
| // M + delta.dt.C + alpha.dt.dt.K |
| Smat->scalgm (dt*dt*alpha); |
| Smat->addgm (1.0,*Mmat); |
| Smat->addgm (dt*delta,*Dmat); |
| } |
| // assembling of dynamic loads |
| nullv (rhs+lcid*n, n); |
| Mb->dlc[lcid].assemble (lcid,rhs,NULL,Ndofm,Mp->time); |
| // total load |
| addv (rhs,fs,n); |
| // modification of vector of prescribed nodal forces |
| cmulv (dt*dt*alpha,rhs,n); |
| // predictor vectors |
| for (j=0;j<n;j++){ |
| dd[j] = lhs[j] + dt*v[j] + dt*dt*(0.5-alpha)*a[j]; |
| vv[j] = v[j] + dt*(1.0-delta)*a[j]; |
| } |
| // M.dd |
| Mmat->gmxv (dd,p); |
| // contribution to the right hand side |
| addv (rhs,p,n); |
| // C.dd |
| Dmat->gmxv (dd,p); |
| cmulv (dt*delta,p,n); |
| addv (rhs,p,n); |
| // C.vv |
| Dmat->gmxv (vv,p); |
| cmulv (-1.0*dt*dt*alpha,p,n); |
| addv (rhs,p,n); |
| // solution of system of algebraic equations |
| //Smat->solve_system (Gtm,lhs,rhs); |
| Mp->ssle->solve_system (Gtm,Smat,lhs,rhs,Out); |
| for (j=0;j<n;j++){ |
| v[j] = vv[j] + delta/dt/alpha*(lhs[j]-dd[j]); |
| a[j] = 1.0/dt/dt/alpha*(lhs[j]-dd[j]); |
| } |
| compute_req_val (lcid); |
| print_step(lcid, i, Mp->time, rhs); |
| print_flush(); |
| // time increment |
| Mp->time = Mp->timecon.newtime (); |
| dt = Mp->timecon.actualbacktimeincr (); |
| i++; |
| }while(Mp->time<end_time); |
| print_close(); |
| delete [] p; |
| delete [] dd; |
| delete [] vv; |
| } |
| /** |
| difference method for solution of forced dynamic problems |
| the method is explicit |
| /trunk/SIFEL/MEFEL/SRC/fdsolver.h |
|---|
| 7,6 → 7,7 |
| void newmark_method (long lcid); |
| void optim_newmark_method (long lcid); |
| void difference_method (long lcid); |
| /trunk/SIFEL/MEFEL/SRC/generalmod.cpp |
|---|
| 22,7 → 22,9 |
| /****************** double-porosity structure, thermal effects ***********************************/ |
| #include <string.h> |
| #include "global.h" |
| #include "generalmod.h" |
| #include "probdesc.h" |
| using namespace std; |
| 96,6 → 98,8 |
| if (parms) |
| initialise_parameters(parms); |
| if (flag < 0) |
| return 0; |
| //Initialises state variables and parameters (flag 0) |
| if (flag == 0) { |
| 123,6 → 127,10 |
| //Calculate stiffness matrix by perturbation of current strain increment |
| if(flag == 1) { |
| direct_stiffness_matrix(strain_gen_init, stress_gen_init, qstatev_init, dstrain_gen, dtime, DDtan_gen, parms, kinc, flag); |
| //debug??!! |
| //fprintf(Out,"ipp = %ld\n",ipp); |
| //fflush(Out); |
| } |
| } |
| else { |
| 292,6 → 300,12 |
| double corr_T=+rlambda_for_se*gamma*suction*bTparam*dTemper/(aTparam+bTparam*Temper); |
| double chi=SrM; |
| //debug??!! |
| //if(rlambda != 1.0 && suction > -100.0){ |
| //fprintf(Out," rlambda = %e, chi = %e, psi = %e, s = %e\n",rlambda,chi,chi*(1-gamma*rlambda),suction); |
| //fflush(Out); |
| //} |
| for(int i=0; i<3; i++) { |
| for(int j=0; j<3; j++) { |
| dsignet_iter[3*i+j]=dsigMef[3*i+j]-kron_delta[i][j]*chi*(corr_s+corr_eM+corr_T); |
| 1081,7 → 1095,7 |
| double signet[9]; |
| double Sr=0; |
| double othervar[5]={0,0,0,0,0}; |
| convert_from_general(stress_gen, signet, othervar, ngstress-6, 1); |
| Sr=othervar[0]; |
| /trunk/SIFEL/MEFEL/SRC/hypoplunsatexptherm2.cpp |
|---|
| 219,8 → 219,9 |
| The function initializes material model. Retrieves initial values of the following |
| quantities |
| @param ipp - integration point number |
| @param ido - index of internal variables for given material in the ipp other array |
| @param[in] ipp - integration point number |
| @param[in] ido - index of internal variables for given material in the ipp other array |
| @param[in] rinit - flag for initialization after restorage from hdbackup |
| @return The function does not return anything but updates initial value of Sr and stress vector |
| in eqother array. |
| 227,7 → 228,7 |
| Created by Tomas Koudelka, 1.10.2015 |
| */ |
| void hypoplast_unsat_exp_therm::initval (long ipp,long ido) |
| void hypoplast_unsat_exp_therm::initval (long ipp,long ido, bool rinit) |
| { |
| // data handled by MEFEL |
| long ncomp = Mm->ip[ipp].ncompstr; // number of reduced stress/strain components |
| 276,7 → 277,10 |
| eps(6) = statev[1]; // set suction to the corresponding component of generalized strain vector |
| eps(7) = statev[3]; // set temperature to the corresponding component of generalized strain vector |
| // initialize model |
| herr = huet.soil_model(eps.a, sig.a, statev, deps.a, 1.0, NULL, params, NULL, rkfparams, 0, Mp->jstep, ipp); |
| int flag = 0; |
| if (rinit == true) |
| flag = -1; |
| herr = huet.soil_model(eps.a, sig.a, statev, deps.a, 1.0, NULL, params, NULL, rkfparams, flag, Mp->jstep, ipp); |
| if (herr) |
| { |
| print_err("hypoplasticity model cannot be initialized on element %ld (ip=%ld)", __FILE__, __LINE__, __func__, Mm->elip[ipp]+1, ipp+1); |
| 628,6 → 632,9 |
| rkf_statev[1] = Mm->ip[ipp].eqother[ido+2*ncomp+3+nstatev+6]; |
| rkf_statev[2] = Mm->ip[ipp].eqother[ido+2*ncomp+3+nstatev+7]; |
| rkf_statev[3] = Mm->ip[ipp].eqother[ido+2*ncomp+3+nstatev]; |
| //if(Mp->time >= 90899 && ipp == 0) |
| //printf("Time = %lf\n",Mp->time); |
| rkf_stat = huet.soil_model(eps.a, sig.a, statev.a, deps.a, dtime, dd_gen.a, NULL, rkf_statev.a, rkfparams, 1, Mp->jstep, ipp); |
| // actualize RKF state variables |
| /trunk/SIFEL/MEFEL/SRC/hypoplunsatexptherm2.h |
|---|
| 73,7 → 73,7 |
| ~hypoplast_unsat_exp_therm(); |
| void read(XFILE *in); |
| void print(FILE *out); |
| void initval (long ipp,long ido); |
| void initval (long ipp,long ido, bool rinit); |
| void matstiff (matrix &d,long ipp,long ido); |
| void nlstressesincr (long ipp, long im, long ido); |
| void givestressincr (long lcid, long ipp, long im, long ido, long fi, vector &sig); |
| /trunk/SIFEL/MEFEL/SRC/lintet.cpp |
|---|
| 2103,7 → 2103,52 |
| /** |
| function defines meaning of DOFs at nodes |
| @param eid - element id |
| TKr 11.4.2024 according to JK |
| */ |
| void lintet::define_meaning (long eid) |
| { |
| ivector cn(ASTCKIVEC(ndofe)),nod(ASTCKIVEC(nne)); |
| Mt->give_elemnodes (eid,nod); |
| Mt->give_code_numbers (eid,cn.a); |
| // displacement in x direction |
| if (cn[0]>0) Mt->nodes[nod[0]].meaning[0]=1; |
| // displacement in y direction |
| if (cn[1]>0) Mt->nodes[nod[0]].meaning[1]=2; |
| // displacement in z direction |
| if (cn[2]>0) Mt->nodes[nod[0]].meaning[2]=3; |
| // displacement in x direction |
| if (cn[3]>0) Mt->nodes[nod[1]].meaning[0]=1; |
| // displacement in y direction |
| if (cn[4]>0) Mt->nodes[nod[1]].meaning[1]=2; |
| // displacement in z direction |
| if (cn[5]>0) Mt->nodes[nod[1]].meaning[2]=3; |
| // displacement in x direction |
| if (cn[6]>0) Mt->nodes[nod[2]].meaning[0]=1; |
| // displacement in y direction |
| if (cn[7]>0) Mt->nodes[nod[2]].meaning[1]=2; |
| // displacement in z direction |
| if (cn[8]>0) Mt->nodes[nod[2]].meaning[2]=3; |
| // displacement in x direction |
| if (cn[9]>0) Mt->nodes[nod[3]].meaning[0]=1; |
| // displacement in y direction |
| if (cn[10]>0) Mt->nodes[nod[3]].meaning[1]=2; |
| // displacement in z direction |
| if (cn[11]>0) Mt->nodes[nod[3]].meaning[2]=3; |
| } |
| //////////////////// /* termitovo */ //////////////////////////////////// |
| /** |
| /trunk/SIFEL/MEFEL/SRC/lintet.h |
|---|
| 77,6 → 77,7 |
| void intpointval (long eid,vector &nodval,vector &ipval); |
| void res_eigstrain_forces (long eid,vector &nfor); |
| void eigstrain_forces (long eid,vector &nfor); |
| void define_meaning (long eid); |
| /* termitovo */ |
| void ntdbr_vector (long eid,vector &ntdbr); |
| /trunk/SIFEL/MEFEL/SRC/mechmat.cpp |
|---|
| 235,7 → 235,10 |
| delete [] mcoul; delete [] mcpar; |
| delete [] boerm; delete [] drprm; delete [] ddpm; delete [] drprm2; delete [] cclay; |
| delete [] cclayc; delete [] bbm; delete [] shpl; delete [] hisspl; |
| delete [] bbm; |
| delete [] cclayc; |
| delete [] shpl; |
| delete [] hisspl; |
| delete [] glasgmat; delete [] glasgdam; delete [] crdam; delete [] tswmat; delete [] effstr; |
| 2552,12 → 2555,13 |
| Function initializes material models with initial values. |
| @param[in] lcid - load case id |
| @param[in] rinit - flag for initialization after restorage from hdbackup |
| @return The function does not return anything. |
| Created by JK, 21.6.2004 |
| */ |
| void mechmat::initmaterialmodels (long lcid) |
| void mechmat::initmaterialmodels (long lcid, bool rinit) |
| { |
| long i,j,nip,ipp; |
| 2567,7 → 2571,7 |
| // number of the first integration point |
| ipp=Mt->elements[i].ipp[0][0]; |
| for(j=0; j<nip; j++){ |
| initvalues(lcid, ipp+j, 0, 0); |
| initvalues(lcid, ipp+j, 0, 0, rinit); |
| } |
| } |
| } |
| 2580,12 → 2584,7 |
| The function is intended for the transfer of values among meshes in coupled problems especially. |
| @param[in] lcid - load case id |
| @param[in] n - the number of required auxiliary points in the mapping array ipm |
| @param[in] ipm - integration point mapping array, |
| ipm[i].ipp < 0 => auxiliary integration point must be used |
| ipm[i].ipp >= 0 => direct mapping to the regulal integration point, |
| no computation of strains is performed, the strains are assumed |
| to be computed at the main solution procedure of the problem |
| @param[in] rinit - flag for initialization after restorage from hdbackup |
| @return The function does not return anything but it chages state of auxiliary integration points |
| at the array aip. |
| 2592,7 → 2591,7 |
| Created by Tomas Koudelka, 24.11.2017 |
| */ |
| void mechmat::aip_initmaterialmodels (long lcid) |
| void mechmat::aip_initmaterialmodels (long lcid, bool rinit) |
| { |
| long app; |
| intpoints *tmp_ip; |
| 2627,7 → 2626,7 |
| for (app=0; app<tnaip; app++) |
| { |
| if (Gtm->leso[elip[app]] == 1) |
| initvalues(lcid, app, 0, 0); |
| initvalues(lcid, app, 0, 0, rinit); |
| } |
| // switch regular integration point auxiliary integration point arrays |
| 2651,12 → 2650,13 |
| @param[in] ipp - integration point pointer |
| @param[in] im - index of material type for given ip |
| @param[in] ido - index of internal variables for given material in the ipp eqother array |
| @param[in] rinit - flag for initialization after restorage from hdbackup |
| @return The function deos not return anything. |
| Created by Tomas Koudelka |
| */ |
| void mechmat::initvalues (long lcid, long ipp,long im,long ido) |
| void mechmat::initvalues (long lcid, long ipp,long im,long ido, bool rinit) |
| { |
| long ncompo; |
| switch (ip[ipp].tm[im]){ |
| 2685,7 → 2685,7 |
| case druckerprager2: |
| ncompo = givencompeqother (ipp, im); |
| ncompo -= givencompeqother (ipp, im+1); |
| initvalues(lcid, ipp, im+1, ido+ncompo); |
| initvalues(lcid, ipp, im+1, ido+ncompo, rinit); |
| break; |
| case modcamclaymat: |
| if (Mp->strainstate == 0){ |
| 2727,7 → 2727,7 |
| compute_ipstrains(lcid); |
| Mp->strainstate = 1; |
| } |
| hypoplustherm[ip[ipp].idm[im]].initval(ipp, ido); |
| hypoplustherm[ip[ipp].idm[im]].initval(ipp, ido, rinit); |
| break; |
| case shefpl: |
| case chenplast: |
| 2735,16 → 2735,16 |
| case glasgowdamage: |
| break; |
| case creep_damage: |
| crdam[ip[ipp].idm[im]].initvalues(lcid, ipp, im, ido); |
| crdam[ip[ipp].idm[im]].initvalues(lcid, ipp, im, ido, rinit); |
| break; |
| case shrinkagemat: |
| shmat[ip[ipp].idm[im]].initvalues(lcid, ipp, im, ido); |
| shmat[ip[ipp].idm[im]].initvalues(lcid, ipp, im, ido, rinit); |
| break; |
| case time_switchmat: |
| tswmat[ip[ipp].idm[im]].initvalues(lcid, ipp, im, ido); |
| tswmat[ip[ipp].idm[im]].initvalues(lcid, ipp, im, ido, rinit); |
| break; |
| case effective_stress: |
| effstr[ip[ipp].idm[im]].initvalues(lcid, ipp, im, ido); |
| effstr[ip[ipp].idm[im]].initvalues(lcid, ipp, im, ido, rinit); |
| break; |
| case simvisc: |
| case isovisc: |
| 5536,6 → 5536,15 |
| // bulk modulus |
| ret = give_bulk_modulus(ipp); |
| break; |
| case mmean_stress:{ |
| vector sig(ip[ipp].ncompstr); |
| vector sigt(ASTCKVEC(6)); |
| givestress(0,ipp,sig); |
| give_full_vector(sigt,sig,ip[ipp].ssst); |
| // total mean stress |
| ret = first_invar(sigt)/3.0; |
| break; |
| } |
| default:{ |
| print_err("unknown type of quantity is required",__FILE__,__LINE__,__func__); |
| } |
| 9345,6 → 9354,7 |
| @param[in] lcid - load case id obtained from solver |
| @param[in] init - flag for initialization of actualized material model (0 = no call of initvalues, 1=initval is called) |
| @param[in] rinit - flag for initialization after restorage from hdbackup |
| @retval 0 - if no material index has been changed |
| @retval 1 - if some material indeces have been actualized |
| 9351,7 → 9361,7 |
| Created by Tomas Koudelka, 16.3.2016 |
| */ |
| long mechmat::update_actual_mat_id(long lcid, long init) |
| long mechmat::update_actual_mat_id(long lcid, long init, bool rinit) |
| { |
| long i,j,k,ipp,nip,app; |
| long oami,ami,ret = 0; |
| 9403,7 → 9413,7 |
| // if the material model has changed (ami != oami) then the |
| // initialization of the material model must be performed if required (init == 1) |
| if (init) |
| tswmat[ip[ipp].idm[0]].initvalues (lcid, ipp, 0, 0); |
| tswmat[ip[ipp].idm[0]].initvalues (lcid, ipp, 0, 0, rinit); |
| } |
| } |
| } |
| 9445,7 → 9455,7 |
| // if the material model has changed (ami != oami) then the |
| // initialization of the material model must be performed if required (init == 1) |
| if (init) |
| tswmat[aip[app].idm[0]].initvalues (lcid, app, 0, 0); |
| tswmat[aip[app].idm[0]].initvalues (lcid, app, 0, 0, rinit); |
| } |
| } |
| } |
| /trunk/SIFEL/MEFEL/SRC/mechmat.h |
|---|
| 204,13 → 204,13 |
| void allocmacrostresses (); |
| /// initailizes material models on all regular int. points before begining of main computation procedure |
| void initmaterialmodels (long lcid); |
| void initmaterialmodels (long lcid, bool rinit); |
| /// initializes material models on all auxiliary int. points before begining of main computation procedure (used in coupled problems) |
| void aip_initmaterialmodels (long lcid); |
| void aip_initmaterialmodels (long lcid, bool rinit); |
| /// initailizes material models for the given int. point |
| void initvalues (long lcid, long ipp,long im,long ido); |
| void initvalues (long lcid, long ipp,long im,long ido, bool rinit); |
| /// returns the number of the normal components in the stress/strain state of the given int. point |
| long give_num_norm_comp(long ipp); |
| 552,7 → 552,7 |
| void updateothermat (long ipp,long im,long ido); |
| /// updates actual material index and optionaly calls initialization of new activated material |
| long update_actual_mat_id(long lcid, long init); |
| long update_actual_mat_id(long lcid, long init, bool rinit); |
| /// sets all arrays on int. points to zero |
| void clean_ip (); |
| /trunk/SIFEL/MEFEL/SRC/mechtop.cpp |
|---|
| 3687,6 → 3687,10 |
| Pelq->define_meaning (i); |
| break; |
| } |
| case lineartet:{ |
| Ltet->define_meaning (i); |
| break; |
| } |
| default:{ |
| print_err("unknown element type is required",__FILE__,__LINE__,__func__); |
| } |
| /trunk/SIFEL/MEFEL/SRC/mtsolver.cpp |
|---|
| 137,7 → 137,7 |
| // inital printing of output and graphical informations |
| if (Mp->hdbcont.restore_stat()){ |
| // initiation of mechanical material models |
| Mm->initmaterialmodels(lcid); |
| Mm->initmaterialmodels(lcid, true); |
| if (Mespr==1) |
| fprintf (stdout,"\n Reading of MEFEL backup file\n"); |
| solver_restore (r,fp,i,Mp->time,dt,&Mp->timecon,n); |
| 149,7 → 149,7 |
| else |
| { |
| // initiation of mechanical material models |
| Mm->initmaterialmodels(lcid); |
| Mm->initmaterialmodels(lcid, false); |
| print_init(-1, "wt"); |
| print_step(lcid, i, Mp->time, f, fb); |
| //print_step(lcid, i, 0.0, f); |
| 640,9 → 640,9 |
| if (Mp->hdbcont.restore_stat()) |
| { |
| // initiation of mechanical material models |
| Mm->initmaterialmodels(lcid); |
| Mm->initmaterialmodels(lcid, true); |
| // initiation of material models on auxiliary integration points |
| Mm->aip_initmaterialmodels(lcid); |
| Mm->aip_initmaterialmodels(lcid, true); |
| if (Mespr==1) |
| fprintf (stdout,"\n Reading of MEFEL backup file\n"); |
| solver_restore (mt_gv.r, mt_gv.fp, mt_gv.istep, Mp->time, dt, &Mp->timecon, n); |
| 656,7 → 656,7 |
| update_elnod_stat_after_hdbrest(lcid); |
| // update actual material index for all elements |
| Mm->update_actual_mat_id(lcid, 0); |
| Mm->update_actual_mat_id(lcid, 0, true); |
| // set indicator of calculation started from restored backup data |
| rest_calc = 1; |
| // reopen output files |
| 666,15 → 666,15 |
| { |
| mt_gv.istep=0; |
| // initialize actual chain of time switch material models |
| Mm->update_actual_mat_id(lcid, 0); |
| Mm->update_actual_mat_id(lcid, 0, false); |
| // Switch on only elements active at the initial time, |
| // it is needed for the proper initialization of creep models |
| if (Mp->tprob == growing_mech_structure) |
| Gtm->update_elem(Mp->time); |
| // initiation of mechanical material models |
| Mm->initmaterialmodels(lcid); |
| Mm->initmaterialmodels(lcid, false); |
| // initiation of material models on auxiliary integration points |
| Mm->aip_initmaterialmodels(lcid); |
| Mm->aip_initmaterialmodels(lcid, false); |
| // compute stresses from eigenstrains or initial irreversible strains |
| if ((Mp->eigstrains > 0) && (Mp->eigstrains < 4)) |
| internal_forces(lcid, mt_gv.fp); |
| 766,7 → 766,7 |
| // update actual material index for all old elements and call initval for those whose indeces have changed |
| if (rest_calc == 0) |
| { |
| um = Mm->update_actual_mat_id(lcid, 1); |
| um = Mm->update_actual_mat_id(lcid, 1, false); |
| if (um) fusm = yes; |
| else fusm = no; |
| } |
| 1000,10 → 1000,10 |
| // update actual material index for all new elements |
| // but do not initialize new active material models |
| Mm->update_actual_mat_id(lcid, 0); |
| Mm->update_actual_mat_id(lcid, 0, false); |
| // initiation of mechanical material models |
| Mm->initmaterialmodels(lcid); |
| Mm->initmaterialmodels(lcid, false); |
| // |
| /trunk/SIFEL/MEFEL/SRC/nssolver.cpp |
|---|
| 86,7 → 86,7 |
| if ((Mp->nlman->hdbackupal == 2) || (Mp->hdbcont.restore_stat())) |
| { |
| // initiation of mechanical material models |
| Mm->initmaterialmodels(i); |
| Mm->initmaterialmodels(i, true); |
| if (Mespr==1) |
| fprintf (stdout,"\n Reading of backup file\n"); |
| // restorage from the backup is performed |
| 100,7 → 100,7 |
| Mp->lambda = ilambda = 0.0; |
| li=0; |
| // initiation of mechanical material models |
| Mm->initmaterialmodels(i); |
| Mm->initmaterialmodels(i, false); |
| print_init(-1, "wt"); |
| print_step(i, li, Mp->lambda, fa); |
| print_flush(); |
| /trunk/SIFEL/MEFEL/SRC/shrinkmat.cpp |
|---|
| 345,10 → 345,11 |
| @param ipp - integration point pointer |
| @param im - material index |
| @param ido - index of internal variables for given material in the ipp other array |
| @param[in] rinit - flag for initialization after restorage from hdbackup |
| Created 21. 11. 2013 by JK+TKo |
| */ |
| void shrinkmat::initvalues (long lcid, long ipp, long im, long ido) |
| void shrinkmat::initvalues (long lcid, long ipp, long im, long ido, bool rinit) |
| { |
| long ncomp = Mm->ip[ipp].ncompstr; |
| double hum; |
| 361,7 → 362,7 |
| Mm->ip[ipp].eqother[ido] = hum; |
| Mm->initvalues (lcid,ipp,im+1,ido+1+2*ncomp); |
| Mm->initvalues (lcid,ipp,im+1,ido+1+2*ncomp, rinit); |
| } |
| /trunk/SIFEL/MEFEL/SRC/shrinkmat.h |
|---|
| 28,7 → 28,7 |
| void nlstresses (long ipp, long im, long ido); |
| void nonloc_nlstresses (long ipp, long im, long ido); |
| void updateval (long ipp, long im, long ido); |
| void initvalues (long lcid, long ipp, long im, long ido); |
| void initvalues (long lcid, long ipp, long im, long ido, bool rinit); |
| void givestressincr (long lcid, long ipp, long im, long ido, long fi, vector &sig); |
| void giveirrstrains(long ipp, long im, long ido, vector &epsirr); |
| double give_actual_ft (long ipp, long im, long ido); |
| /trunk/SIFEL/MEFEL/SRC/solverm.cpp |
|---|
| 56,7 → 56,7 |
| compute_ipstrains(0); |
| //stress_initdispl(0); |
| if (Mp->tprob == linear_statics) |
| Mm->initmaterialmodels(0); |
| Mm->initmaterialmodels(0, false); |
| nullv(lhs, Ndofm); |
| nullv(rhs, Ndofm); |
| /trunk/SIFEL/MEFEL/SRC/timeswmat.cpp |
|---|
| 595,12 → 595,13 |
| @param ipp - integration point pointer |
| @param im - material index |
| @param ido - index of internal variables for given material in the ipp other array |
| @param[in] rinit - flag for initialization after restorage from hdbackup |
| @return The function deos not return anything. |
| Created by Tomas Koudelka, 21.11.2006 |
| */ |
| void timeswmat::initvalues (long lcid, long ipp, long /*im*/, long ido) |
| void timeswmat::initvalues (long lcid, long ipp, long /*im*/, long ido, bool rinit) |
| { |
| long ncompstr = 0L; |
| 608,7 → 609,7 |
| ncompstr = Mm->ip[ipp].ncompstr; |
| if (ami > 0) |
| Mm->initvalues(lcid, ipp, ami, ido+ncompstr); |
| Mm->initvalues(lcid, ipp, ami, ido+ncompstr, rinit); |
| else |
| { |
| print_err("no material is switched on (elem=%ld, ipp=%ld, time=%le)", |
| /trunk/SIFEL/MEFEL/SRC/timeswmat.h |
|---|
| 43,7 → 43,7 |
| //double give_pore_press(long ipp, long im, long ido); |
| long givencompother (long ipp, long im); |
| long givencompeqother (long ipp, long im); |
| void initvalues(long lcid, long ipp, long im, long ido); |
| void initvalues(long lcid, long ipp, long im, long ido, bool rinit); |
| void matstiff (matrix &d,long ipp, long im, long ido); |
| void givestressincr (long lcid,long ipp,long im,long ido,long fi,vector &sig); |
| void nlstressesincr (long ipp, long im, long ido); |
| /trunk/SIFEL/METR/SRC/consol_awf1c.cpp |
|---|
| 24,6 → 24,12 |
| { |
| alpha = 1.0;//Biot's constant [-] alpha = 1 - kt/ks |
| // |
| vol_strain_effectc = 0; //volumetric strain rate influence: 0=no; 1=yes |
| // |
| wrc_vol_strain_effectc = 0; //volumetric strain rate influence on water retention curve: 0=no; 1=yes |
| // |
| pore_press_effectc = 0; //pore pressure effect acording to effective stress concept: 0=no; 1=yes |
| // |
| rhos0 = 2000.0;//initial solid grain density [kg/m^3] |
| // |
| phi0 = 0.297;//initial porosity [-] |
| 59,7 → 65,7 |
| void con_awf1matc::read(XFILE *in) |
| { |
| xfscanf (in,"%k%m","waterflowmechtype",&waterflowmechtype_kwdset, &model_type); //water flow model type |
| xfscanf (in,"%lf %lf %lf %lf %d %d", &alpha, &phi0, &rhos0, &t0, &sr_type, &xi_type); |
| xfscanf (in,"%lf %d %d %d %lf %lf %lf %d %d", &alpha, &vol_strain_effectc, &pore_press_effectc, &wrc_vol_strain_effectc, &phi0, &rhos0, &t0, &sr_type, &xi_type); |
| switch (model_type){ |
| case lewis_and_schrefler_coup: |
| 1127,7 → 1133,14 |
| 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 |
| kuw = 0.0;// pore pressure effect included in MEFEL |
| if(pore_press_effectc == 1){ |
| // this below in not completed yet |
| chi = 0.0; |
| //chi = Tm->givenontransq(eff_stress_param, ipp); //actual effective stress parameter from MEFEL |
| kuw = -chi;// pore pressure effect included in METR |
| } |
| break; |
| } |
| default:{ |
| 1191,24 → 1204,31 |
| { |
| double capwu; |
| double sw,sg,rhogw; |
| //double n,dsr_depsv; |
| double n,dsr_depsv; |
| sw = get_sw(pw,ipp); |
| rhogw = get_rhogw(pw); |
| sg = 1.0 - sw; |
| capwu = (sg*rhogw + sw*rhow0)*alpha;//p. 394 or p. 95 multiplied by rhogw and plus vapour effect |
| capwu = 0.0; |
| switch (model_type){ |
| case lewis_and_schrefler_coup:{//Lewis and Schrefler's model |
| capwu = (sg*rhogw + sw*rhow0)*alpha;//p. 394 or p. 95 multiplied by rhogw and plus vapour effect |
| break; |
| } |
| case lewis_and_schrefler_mefel_coup:{//Lewis and Schrefler's model and moisture storage functions is included in MEFEL |
| //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 |
| capwu = 0.0;// volumetric strain effect included in TRFEL |
| if(vol_strain_effectc == 1){ |
| sw = get_sw(pw,ipp); |
| capwu = (sg*rhogw + sw*rhow0)*alpha;//volumetric strain effect included in METR |
| if(wrc_vol_strain_effectc == 1){ |
| // this must be rewritten for METR elements |
| n = get_phi(pw,ipp); |
| dsr_depsv = Tm->givenontransq(der_saturation_deg_depsv, ipp); //actual derivative of saturation degree with respect to volumetric strain |
| capwu = capwu + n*(rhow0 - rhogw)*dsr_depsv;//volumetric strain effect on retention curve |
| } |
| } |
| break; |
| } |
| default:{ |
| /trunk/SIFEL/METR/SRC/consol_awf1c.h |
|---|
| 95,7 → 95,7 |
| private: |
| waterflowmechtype model_type; |
| int sr_type,xi_type; |
| int vol_strain_effectc,pore_press_effectc,wrc_vol_strain_effectc,sr_type,xi_type; |
| double alpha,phi0,rhos0,rhow0,t0; |
| double gamma,lambda0,s_entry; |
| double c8,c9,c10,c11,c12,c13; |
| /trunk/SIFEL/METR/SRC/consol_hawf3c.cpp |
|---|
| 35,9 → 35,16 |
| con_hawf3matc::con_hawf3matc() |
| { |
| compress = 0; //compressible grains: 0=no; 1=yes |
| // |
| //STATE VARIABLES |
| vol_strain_effectc = 0; //volumetric strain rate influence: 0=no; 1=yes |
| wrc_vol_strain_effectc = 0; //volumetric strain rate influence on water retention curve: 0=no; 1=yes |
| pore_press_effectc = 0; //pore pressure effect acording to effective stress concept: 0=no; 1=yes |
| temper_effectc = 0; //temperature effect on displacements |
| sr_type = 1; //retention curve calculation type |
| xi_type = 1; //effective stress parameter type |
| betas_type = 0; //thermal expansion calculation type |
| //STATE VARIABLES |
| mw = 18.01528; //molar mass of water kg.mol-1 |
| ma = 28.9645; //molar mass of dry air kg.mol-1 |
| gasr = 8314.41; //universal gas constant J.mol-1.K-1 |
| 45,16 → 52,9 |
| t0 = 273.15; |
| p0 = 101325.0; |
| // PHYSICAL PROPERTIES OF DRY AIR |
| muga0 = 17.17e-6;//[Pa.s] |
| alphaa = 4.73e-8;//[Pa.s.K-1] |
| betaa = 2.222e-11;//[Pa.s.K-2] |
| // PHYSICAL PROPERTIES OF WATER VAPOUR |
| // from D.Gawin, F.Pesavento (PRVAP.f90) |
| //Hyland and Wexler equation (1983) for the saturation water vapor pressure |
| dv0 = 2.58e-5; //effective diffusion coefficient of vapour m^2/s at reference temperature 273.15 |
| bv = 1.667; |
| c8 = -5.8002206e+03; |
| c9 = 1.3914993; |
| c10 =-4.8640239e-02; |
| 61,8 → 61,6 |
| c11 = 4.1764768e-05; |
| c12 = -1.4452093e-08; |
| c13 = 6.5459673; |
| mugw0 = 8.85e-6;//[Pa.s] |
| alphaw = 3.633e-8;//[Pa.s.K-1] |
| // M. Starnoni 10-11-2010 |
| // Gaw-Maj-Sch "Numerical analysis of hygro thermal behaviour and damage of concrete at high T" |
| 73,8 → 71,6 |
| //rhow0 = 999.84;//testing??!! |
| rhow0 = 1000.0;//testing??!! |
| tcr = 647.3; |
| cwat = 0.2e9; |
| betawat = -0.414e-3; |
| hvap0 = 2.7e+5; |
| // M. Starnoni 10-11-2010 |
| 91,33 → 87,20 |
| b3 = -9.2118e-5; b4 = 3.3534e-7; b5 = -4.4034e-10; |
| pr1 = 1.0e7; prif = 2.0e7; |
| //muw0 = 0.6612;//testing??!! |
| muw0 = 1.0e-3; |
| conb = 229.0; |
| conc = -1.562; |
| cpw = 4181.0; |
| lambdaw = 0.6; |
| //kw0 = 0.43e10;//compresibility coefficient of water - this is nonsense //09/06/2017 |
| kw0 = 2.0e9;//bulk modulus of water |
| //PHYSICAL PROPERTIES OF SOIL set to zero |
| alpha0 = 0.0; |
| ks0 = 0.0; |
| kt0 = 0.0; |
| phi0 = 0.0; |
| kintr0 = 0.0; |
| betas0 = 0.0; |
| deff0 = 0.0; |
| betas_dry = 0.0; |
| betas_wet = 0.0; |
| rhos0 = 0.0; |
| cps0 = 0.0; |
| lambda0 = 0.0; |
| lambda_dry = 0.0; |
| lambda_wet = 0.0; |
| sirr0 = 0.0; |
| tau0 = 0.0; |
| emod = 1.3e6; |
| nu = 0.25; |
| gamma = 0.0; |
| lambda0 = 0.0; |
| s_entry = 0.0; |
| } |
| con_hawf3matc::~con_hawf3matc() |
| 134,28 → 117,106 |
| void con_hawf3matc::read(XFILE *in) |
| { |
| xfscanf (in,"%k%m","heatairwaterflowmechtype",&heatairwaterflowmechtype_kwdset, &model_type); |
| xfscanf (in,"%lf %d %d %d %d %lf %lf %lf %lf %d %d %d", &alpha0, &vol_strain_effectc, &pore_press_effectc, &wrc_vol_strain_effectc, &temper_effectc, &phi0, &rhos0, &emod, &nu, &sr_type, &xi_type, &betas_type); |
| switch (model_type){ |
| case lewis_and_schrefler3_coup:{//Lewis and Schrefler's book |
| xfscanf (in,"%lf %lf %lf %lf %lf %lf %lf %lf %lf %lf", &alpha0, &ks0, &phi0, &kintr0, &betas0, &rhos0, &cps0, &lambda0, &emod, &nu); |
| case lewis_and_schrefler3_coup: |
| case lewis_and_schrefler3_mefel_coup: {//Lewis and Schrefler's book |
| //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 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; |
| } |
| case masin_sr:{//extended formulation from Brooks and Correy according to Masin |
| masin_ret.read(in); |
| break; |
| } |
| case febex_granit_sr:{//FEBEX granit |
| febex_granit_ret.read(in); |
| break; |
| } |
| default:{ |
| print_err ("unknown model type is required", __FILE__, __LINE__, __func__); |
| abort(); |
| } |
| } |
| switch (xi_type){ |
| case biot_xi: |
| break; |
| case biot_reduced_xi:{//coefficient for effective stress factor reduction |
| xfscanf (in,"%le %le", &gamma,&lambda0); |
| break; |
| } |
| case biot_masin_xi:{ |
| xfscanf (in,"%le %le", &gamma,&s_entry); |
| break; |
| } |
| case masin_xi: |
| 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(); |
| } |
| } |
| break; |
| } |
| case lewis_and_schrefler3_mefel_coup:{//Lewis and Schrefler's book |
| xfscanf (in,"%lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf", &mefel_units, &alpha0, &ks0, &phi0, &kintr0, &betas0, &rhos0, &cps0, &lambda_dry, &lambda_wet, &sirr0, &tau0); //, &emod, &nu);//provisionally??!! |
| break; |
| } |
| case lewis_and_schrefler3_2_coup:{//Lewis and Schrefler's book p. 381 |
| xfscanf (in,"%lf %lf %lf %lf %lf %lf %lf %lf", &alpha0, &ks0, &phi0, &kintr0, &betas0, &rhos0, &cps0, &lambda0); |
| break; |
| } |
| default:{ |
| fprintf (stderr,"\n unknown model type is required"); |
| fprintf (stderr,"\n in function (%s, line %d).\n",__FILE__,__LINE__); |
| } |
| } |
| // incompressible grains: |
| if(compress == 0) |
| alpha0 = 1.0; |
| } |
| 168,25 → 229,26 |
| */ |
| void con_hawf3matc::print(FILE *out) |
| { |
| fprintf (out,"\n %d ", int(model_type)); |
| switch (model_type){ |
| case lewis_and_schrefler3_coup:{//Lewis and Schrefler's book |
| fprintf (out,"\n %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf\n",alpha0,ks0,phi0,kintr0,betas0,rhos0,cps0,lambda0,emod,nu); |
| break; |
| } |
| case lewis_and_schrefler3_mefel_coup:{//Lewis and Schrefler's book |
| fprintf (out,"\n %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf \n",mefel_units,alpha0,ks0,phi0,kintr0,betas0,rhos0,cps0,lambda_dry,lambda_wet,sirr0,tau0); |
| break; |
| } |
| case lewis_and_schrefler3_2_coup:{//Lewis and Schrefler's book p. 381 |
| fprintf (out,"\n %lf %lf %lf %lf %lf %lf %lf %lf \n",alpha0,ks0,phi0,kintr0,betas0,rhos0,cps0,lambda0); |
| break; |
| } |
| default:{ |
| fprintf (stderr,"\n unknown model type is required"); |
| fprintf (stderr,"\n in function (%s, line %d).\n",__FILE__,__LINE__); |
| } |
| } |
| /* fprintf (out,"\n %d ", int(model_type)); |
| switch (model_type){ |
| case lewis_and_schrefler3_coup:{//Lewis and Schrefler's book |
| fprintf (out,"\n %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf\n",alpha0,ks0,phi0,kintr0,betas0,rhos0,cps0,lambda0,emod,nu); |
| break; |
| } |
| case lewis_and_schrefler3_mefel_coup:{//Lewis and Schrefler's book |
| fprintf (out,"\n %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf \n",mefel_units,alpha0,ks0,phi0,kintr0,betas0,rhos0,cps0,lambda_dry,lambda_wet,sirr0,tau0); |
| break; |
| } |
| case lewis_and_schrefler3_2_coup:{//Lewis and Schrefler's book p. 381 |
| fprintf (out,"\n %lf %lf %lf %lf %lf %lf %lf %lf \n",alpha0,ks0,phi0,kintr0,betas0,rhos0,cps0,lambda0); |
| break; |
| } |
| default:{ |
| fprintf (stderr,"\n unknown model type is required"); |
| fprintf (stderr,"\n in function (%s, line %d).\n",__FILE__,__LINE__); |
| } |
| } |
| */ |
| } |
| /** |
| 1035,41 → 1097,161 |
| TKr 24/11/2017 |
| */ |
| double con_hawf3matc::get_sw(double pw, double pg, double /*t*/, long ipp) |
| double con_hawf3matc::get_sw(double pw, double pg, double t, long ipp) |
| { |
| double sw,pc; |
| double lam,sirr,pb,se; |
| //capillary pressure |
| pc = pg - pw; |
| sw = 0.0; |
| pc = 0.0; |
| switch (model_type){ |
| case lewis_and_schrefler3_coup:{//Lewis and Schrefler's book - Liakopoulos test |
| sw = lewis_ret.sw(pw); |
| switch (sr_type){ |
| case bazant_sr:{//Bazant |
| pc = pg - pw; |
| sw = bazant_ret.sat(pc,t); |
| break; |
| } |
| case lewis_and_schrefler3_2_coup:{//Lewis and Schrefler's book p. 381 |
| lam = 2.308; |
| sirr = 0.3216; |
| pb = 133813.0; |
| se = pow(pb/pc,lam); |
| sw = se*(1 - sirr) + sirr;//p. 381 |
| case lewis_and_schrefler_sr:{//Lewis and Schrefler's book |
| pg = pg - p0;//pore gas pressure without atmospheric pressure |
| pc = pg - pw; |
| sw = lewis_ret.sw(-pc); |
| break; |
| } |
| case lewis_and_schrefler3_mefel_coup:{//saturation degree ind its derivative are obtained from mefel; |
| sw = Tm->givenontransq(saturation_deg, ipp); //actual saturation degree from models included in TRFEL |
| 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 van_genuchten_sr:{//partially saturated medium = Van Genuchten model |
| //pg = pg - p0;//pore gas pressure without atmospheric pressure//debug?? |
| pc = pg - pw; |
| sw = van_genuchten_ret.sw(-pc,t);//suction = capillary pressure |
| break; |
| } |
| case van_genuchten2_sr:{//partially saturated medium = Van Genuchten model type 2 |
| //sw = van_genuchten_ret.sw2(pw); |
| break; |
| } |
| case masin_sr:{//extended formulation from Brooks and Correy according to Masin |
| double e=0.0,dpw=0.0,dpg=0.0,dpc=0.0,pc=0.0,por=0.0; |
| if(Tm->nontransq != NULL){ |
| por = Tm->givenontransq(porosity, ipp);// from mefel |
| e = por/(1-por); |
| } |
| else{ |
| e = phi0/(1-phi0); |
| } |
| dpw = Tm->ip[ipp].av[0]-Tm->ip[ipp].pv[0]; |
| dpg = Tm->ip[ipp].av[1]-Tm->ip[ipp].pv[1]; |
| dpc = dpg - dpw; |
| pg = pg - p0;//pore gas pressure without atmospheric pressure |
| pc = pg - pw; |
| sw = masin_ret.sw(pc,dpc,e,t);//positive value of suction |
| break; |
| } |
| case febex_granit_sr:{//FEBEX granit |
| //pg = pg - p0;//pore gas pressure without atmospheric pressure//debug?? |
| pc = pg - pw; |
| sw = febex_granit_ret.sw(pc); |
| break; |
| } |
| default:{ |
| fprintf (stderr,"\n unknown model type is required"); |
| fprintf (stderr,"\n in function (%s, line %d).\n",__FILE__,__LINE__); |
| print_err ("unknown model type is required", __FILE__, __LINE__, __func__); |
| abort(); |
| } |
| } |
| } |
| return(sw); |
| } |
| /** |
| function computes effective stress factor xi |
| @param pw - water pressure |
| @param ipp - number of integration point |
| @retval xi - factor xi |
| 13/11/2023, TKr |
| */ |
| double con_hawf3matc::get_xi(double pw, double pg, double t, long ipp) |
| { |
| double suc=0.0,sr=0.0,xi=0.0; |
| switch (xi_type){ |
| case biot_xi: |
| xi = get_sw(pw,pg,t,ipp); |
| break; |
| case biot_reduced_xi: |
| //xi = gamma*get_sw(pw,ipp); |
| sr = get_sw(pw,pg,t,ipp); |
| //xi = pow(sr,(gamma/lambda0)); |
| xi = pow(sr,gamma); |
| xi = (1-gamma)*xi; |
| break; |
| case biot_masin_xi:{//according to masin for testing |
| sr = get_sw(pw,pg,t,ipp); |
| suc = pg - p0 - pw;//pore gas pressure without atmospheric pressure |
| if (suc>=s_entry) |
| xi = pow((s_entry/suc),gamma); |
| else |
| xi = 1.0; |
| if (suc>=s_entry) |
| xi = (1-gamma)*xi; |
| else |
| xi = 1.0; |
| break; |
| } |
| case masin_xi:{ |
| double e=0.0,dpw=0.0,dpg=0.0,dpc=0.0,pc=0.0,por=0.0; |
| if(Tm->nontransq != NULL){ |
| por = Tm->givenontransq(porosity, ipp);// from mefel |
| e = por/(1-por); |
| } |
| else{ |
| e = phi0/(1-phi0); |
| } |
| dpw = Tm->ip[ipp].av[0]-Tm->ip[ipp].pv[0]; |
| dpg = Tm->ip[ipp].av[1]-Tm->ip[ipp].pv[1]; |
| dpc = dpg - dpw; |
| pg = pg - p0;//pore gas pressure without atmospheric pressure |
| pc = pg - pw; |
| xi = masin_ret.psi(pc,dpc,e,t);//positive value of suction |
| break; |
| } |
| default:{ |
| print_err ("unknown model type is required", __FILE__, __LINE__, __func__); |
| abort(); |
| } |
| } |
| return(xi); |
| } |
| /** |
| function computes Biot's constant |
| @param pw - pore water pressure |
| @param pg - pore gas pressure |
| 1089,7 → 1271,7 |
| //alpha = 1.0 - kt/ks; |
| alpha = alpha0;//provisionally??!! |
| alpha = alpha0; |
| return(alpha); |
| } |
| 1105,21 → 1287,24 |
| */ |
| double con_hawf3matc::get_kuw(double pw, double pg, double t,long ipp) |
| { |
| double kuw; |
| double sw,alpha; |
| double chi,kuw; |
| kuw = 0.0; |
| switch (model_type){ |
| case lewis_and_schrefler3_coup:{//Lewis and Schrefler's book |
| sw = get_sw(pw,pg,t,ipp); |
| alpha = get_alpha(pw,pg,t,ipp); |
| kuw = -sw*alpha;//p. 362, but non-incremental eq. for displ. |
| chi = get_xi(pw,pg,t,ipp); |
| kuw = -chi;//p. 362, but non-incremental eq. for displ.; effective stress param. |
| break; |
| } |
| case lewis_and_schrefler3_mefel_coup:{//Lewis and Schrefler's book |
| sw = get_sw(pw,pg,t,ipp); |
| alpha = get_alpha(pw,pg,t,ipp); |
| kuw = -sw*alpha;//p. 362, but non-incremental eq. for displ. |
| kuw = 0.0;// this part is in MEFEL already included |
| case lewis_and_schrefler3_mefel_coup:{//Lewis and Schrefler's model and moisture storage functions is included in MEFEL |
| kuw = 0.0;// pore pressure effect included in MEFEL |
| if(pore_press_effectc == 1){ |
| // this below in not completed yet |
| chi = 0.0; |
| //chi = Tm->givenontransq(eff_stress_param, ipp); //actual effective stress parameter from MEFEL |
| kuw = -chi;// pore pressure effect included in METR |
| } |
| break; |
| } |
| default:{ |
| 1144,23 → 1329,26 |
| */ |
| double con_hawf3matc::get_kug(double pw, double pg, double t,long ipp) |
| { |
| double kug; |
| double sw,alpha,sg; |
| double chi,kug; |
| double sg; |
| kug = 0.0; |
| switch (model_type){ |
| case lewis_and_schrefler3_coup:{//Lewis and Schrefler's book |
| sw = get_sw(pw,pg,t,ipp); |
| alpha = get_alpha(pw,pg,t,ipp); |
| sg = 1.0 - sw; |
| kug = -sg*alpha;//p. 362, but non-incremental eq. for displ. |
| case lewis_and_schrefler3_coup:{//Lewis and Schrefler's book; effective stress param. |
| chi = get_xi(pw,pg,t,ipp); |
| sg = 1.0 - chi; |
| kug = -sg;//p. 362, but non-incremental eq. for displ. |
| break; |
| } |
| case lewis_and_schrefler3_mefel_coup:{//Lewis and Schrefler's book |
| sw = get_sw(pw,pg,t,ipp); |
| alpha = get_alpha(pw,pg,t,ipp); |
| sg = 1.0 - sw; |
| kug = -sg*alpha;//p. 362, but non-incremental eq. for displ. |
| kug = 0.0;// this part is in MEFEL already included |
| case lewis_and_schrefler3_mefel_coup:{//Lewis and Schrefler's model and moisture storage functions is included in MEFEL |
| kug = 0.0;// pore pressure effect included in MEFEL |
| if(pore_press_effectc == 1){ |
| // this below in not completed yet |
| chi = 0.0; |
| //chi = Tm->givenontransq(eff_stress_param, ipp); //actual effective stress parameter from MEFEL |
| kug = -chi;// pore pressure effect included in METR |
| } |
| break; |
| } |
| default:{ |
| 1168,7 → 1356,7 |
| fprintf (stderr,"\n in function (%s, line %d).\n",__FILE__,__LINE__); |
| } |
| } |
| return(kug); |
| } |
| 1246,26 → 1434,40 |
| { |
| double capwu; |
| double sw,sg,alpha,rhow,rhogw; |
| //double n,dsr_depsv; |
| double n,dsr_depsv; |
| sw = get_sw(pw,pg,t,ipp); |
| alpha = get_alpha(pw,pg,t,ipp); |
| rhogw = get_rhogw(pw,pg,t); |
| rhow = get_rhow(t); |
| sg = 1.0 - sw; |
| capwu = 0.0; |
| capwu = (sg*rhogw + sw*rhow)*alpha;//p. 394 or p. 95 multiplied by rhogw and plus vapour effect |
| switch (model_type){ |
| case lewis_and_schrefler3_coup:{//Lewis and Schrefler's book |
| sw = get_sw(pw,pg,t,ipp); |
| alpha = get_alpha(pw,pg,t,ipp); |
| rhogw = get_rhogw(pw,pg,t); |
| rhow = get_rhow(t); |
| sg = 1.0 - sw; |
| capwu = (sg*rhogw + sw*rhow)*alpha;//p. 394 or p. 95 multiplied by rhogw and plus vapour effect |
| 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 |
| capwu = 0.0;// this part is in TRFEL already included |
| case lewis_and_schrefler3_mefel_coup:{//Lewis and Schrefler's model and moisture storage functions is included in MEFEL |
| capwu = 0.0;// volumetric strain effect included in TRFEL |
| if(vol_strain_effectc == 1){ |
| sw = get_sw(pw,pg,t,ipp); |
| alpha = get_alpha(pw,pg,t,ipp); |
| rhogw = get_rhogw(pw,pg,t); |
| rhow = get_rhow(t); |
| sg = 1.0 - sw; |
| capwu = (sg*rhogw + sw*rhow)*alpha;//p. 394 or p. 95 multiplied by rhogw and plus vapour effect |
| if(wrc_vol_strain_effectc == 1){ |
| // this must be rewritten for METR elements |
| 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 |
| } |
| } |
| break; |
| } |
| default:{ |
| 1312,27 → 1514,42 |
| { |
| double capgu; |
| double sg,sw,alpha,rhoga; |
| //double n,dsr_depsv; |
| sw = get_sw(pw,pg,t,ipp); |
| alpha = get_alpha(pw,pg,t,ipp); |
| rhoga = get_rhoga(pw,pg,t); |
| double n,dsr_depsv; |
| sg = 1.0-sw; |
| capgu = sg*rhoga*alpha;//p. 95 multiplied by rhoga or p. 395 |
| capgu = 0.0; |
| switch (model_type){ |
| case lewis_and_schrefler3_coup:{//Lewis and Schrefler's book |
| sw = get_sw(pw,pg,t,ipp); |
| alpha = get_alpha(pw,pg,t,ipp); |
| rhoga = get_rhoga(pw,pg,t); |
| sg = 1.0-sw; |
| capgu = sg*rhoga*alpha;//p. 95 multiplied by rhoga or p. 395 |
| 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 |
| case lewis_and_schrefler3_mefel_coup:{//Lewis and Schrefler's model and moisture storage functions is included in MEFEL |
| capgu = 0.0;// this part is in TRFEL already included |
| if(vol_strain_effectc == 1){ |
| sw = get_sw(pw,pg,t,ipp); |
| alpha = get_alpha(pw,pg,t,ipp); |
| rhoga = get_rhoga(pw,pg,t); |
| sg = 1.0-sw; |
| capgu = sg*rhoga*alpha;//p. 95 multiplied by rhoga or p. 395 |
| if(wrc_vol_strain_effectc == 1){ |
| // this must be rewritten for METR elements |
| 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 |
| } |
| } |
| break; |
| } |
| default:{ |
| 1358,6 → 1575,8 |
| double kut; |
| double betas; |
| kut = 0.0; |
| switch (model_type){ |
| case lewis_and_schrefler3_coup:{//Lewis and Schrefler's book |
| betas = get_betas(ipp); |
| 1364,11 → 1583,13 |
| kut = -1.0*betas/3.0;//p. 348 |
| break; |
| } |
| case lewis_and_schrefler3_mefel_coup:{//Lewis and Schrefler's book |
| betas = get_betas(ipp); |
| kut = -1.0*betas/3.0;//p. 348 |
| case lewis_and_schrefler3_mefel_coup:{//Lewis and Schrefler's model and moisture storage functions is included in MEFEL |
| kut = 0.0;// this part is in MEFEL already included |
| if(temper_effectc == 1){ |
| betas = get_betas(ipp); |
| kut = -1.0*betas/3.0;//p. 348//temperature effect included in METR |
| } |
| break; |
| } |
| default:{ |
| 1430,25 → 1651,37 |
| { |
| double captu; |
| double dhvap,sw,rhow,alpha; |
| //double n,dsr_depsv; |
| double n,dsr_depsv; |
| captu = 0.0; |
| dhvap = get_dhvap(t); |
| sw = get_sw(pw,pg,t,ipp); |
| rhow = get_rhow(t); |
| alpha = get_alpha(pw,pg,t,ipp); |
| captu = -1.0*dhvap*sw*rhow*alpha;//p. 395 but non-incremental eq. for displ. |
| switch (model_type){ |
| case lewis_and_schrefler3_coup:{//Lewis and Schrefler's book |
| dhvap = get_dhvap(t); |
| sw = get_sw(pw,pg,t,ipp); |
| rhow = get_rhow(t); |
| alpha = get_alpha(pw,pg,t,ipp); |
| captu = -1.0*dhvap*sw*rhow*alpha;//p. 395 but non-incremental eq. for displ. |
| 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 |
| captu = 0.0;// this part is in TRFEL already included |
| case lewis_and_schrefler3_mefel_coup:{//Lewis and Schrefler's model and moisture storage functions is included in MEFEL |
| captu = 0.0;// volumetric strain effect included in TRFEL |
| if(vol_strain_effectc == 1){ |
| dhvap = get_dhvap(t); |
| sw = get_sw(pw,pg,t,ipp); |
| rhow = get_rhow(t); |
| alpha = get_alpha(pw,pg,t,ipp); |
| captu = -1.0*dhvap*sw*rhow*alpha;//p. 395 but non-incremental eq. for displ. |
| if(wrc_vol_strain_effectc == 1){ |
| // this must be rewritten for METR elements |
| 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 |
| } |
| } |
| break; |
| } |
| default:{ |
| 1457,7 → 1690,6 |
| } |
| } |
| captu = 0.0;//debug??!! |
| return(captu); |
| } |
| 1822,6 → 2054,8 |
| double fuw1; |
| double phi,rhos,sw,rhow,rhog; |
| fuw1 = 0.0; |
| phi = get_porosity(pw,pg,t,ipp); |
| rhos = get_rhos(t); |
| sw = get_sw(pw,pg,t,ipp); |
| 1839,6 → 2073,7 |
| } |
| case lewis_and_schrefler3_mefel_coup:{//saturation degree ind its derivative are obtained from mefel; |
| fuw1 = -1.0*((1.0 - phi)*rhos + phi*sw*rhow + phi*(1.0 - sw)*rhog); |
| fuw1 = 0.0;//temporarily |
| break; |
| } |
| default:{ |
| 1901,18 → 2136,21 |
| { |
| double fut2; |
| double betas; |
| fut2 = 0.0; |
| betas = get_betas(ipp); |
| switch (model_type){ |
| case lewis_and_schrefler3_coup:{//Lewis and Schrefler's book |
| betas = get_betas(ipp); |
| fut2 = -1.0*betas/3.0; |
| break; |
| } |
| case lewis_and_schrefler3_mefel_coup:{//Lewis and Schrefler's model |
| fut2 = -1.0*betas/3.0; |
| case lewis_and_schrefler3_mefel_coup:{//Lewis and Schrefler's model and moisture storage functions is included in MEFEL |
| fut2 = 0.0;// this part is in MEFEL already included |
| if(temper_effectc == 1){ |
| betas = get_betas(ipp); |
| fut2 = -1.0*betas/3.0;//p. 348//temperature effect included in METR |
| } |
| break; |
| } |
| default:{ |
| 1984,28 → 2222,6 |
| */ |
| double con_hawf3matc::get_pgws(double t) |
| { |
| /* //Clausius-Clapeyron equation |
| double pgws,pgws0,t0,dhvap; |
| switch (Tm->ip[ipp].tm){//type of material |
| case concreteB:{ |
| Tm->concrete[i].concreteB_pgws0(pgws0,t0);//funkce z concreteB.cpp |
| break; |
| } |
| default:{ |
| print_err ("unknown material type is required", __FILE__, __LINE__, __func__); |
| abort(); |
| } |
| } |
| dhvap = get_dhvap(pc,pg,t); |
| pgws = pgws0*exp(-mw*dhvap/gasr*(1.0/t - 1.0/t0)); |
| */ |
| double t1,t2,t3,pgws,psl; |
| t1 = 1.0/t; |
| 2034,55 → 2250,15 |
| double con_hawf3matc::get_rhow(double t) |
| { |
| double rhow,tem; |
| //double pw; |
| /* |
| //nuw = 1.0/rhow; //water specific volume |
| //betaw = 1.0/nuw*dnuw_dt;//volume thermal expansion coefficient |
| //alphaw = 1.0/nuw*dnuw_dp;//isothermal compression modulus of water |
| betaw = 0.68e-4;//[K-1] at t=273.15 |
| //betaw = 10.1e-4;//[K-1] at t=420.0 |
| alphaw = 4.3e-9;//[Pa-1] |
| pw = pc - pg; |
| rhow = rhow0*(1.0 - betaw*(t-t0) + alphaw*(pw - p0)); |
| */ |
| if (t < tcr){ |
| //rhow = rhow0 *(1.+(pg - p0)/cwat + betawat*(t-t0)); |
| //drhow_dpg = rhow0/cwat; |
| //drhow_dpc = 0.0; |
| //drhow_dt = rhow0*betawat; |
| //am_rhow = rhow0 *(1.+ betawat*(t-t0)); |
| tem = t - t0; |
| rhow = (b0+(b1+(b2+(b3+(b4+b5*tem)*tem)*tem)*tem)*tem) + (pr1-prif)* |
| (a0+(a1+(a2+(a3+(a4+a5*tem)*tem)*tem)*tem)*tem); |
| //drhow_dpg = 0.0; |
| //drhow_dpc = 0.0; |
| //am_drhow_dt = rhow0*betawat; |
| //drhow_dt = (b1+3(2*b2+(3*b3+(4*b4+5*b5*tem)*tem)*tem)*tem) + (pr1-prif)* |
| //(a1+(2*a2+(3*a3+(4*a4+5*a5*tem)*tem)*tem)*tem); |
| //am_d2rhow_d2t = 0.0; |
| //d2rhow_d2t = (2*b2+(6*b3+(12*b4+20*b5*tem)*tem)*tem) + (pr1-prif)* |
| //(2*a2+(6*a3+(12*a4+20*a5*tem)*tem)*tem); |
| //d2rhow_d2pc = 0.0; |
| //d2rhow_dtdpc = 0.0; |
| } |
| else{ |
| //rhow = rhow0 *(1.+(pg-p0)/cwat+betawat*(tcr-t0)); |
| //drhow_dpg = rhow0/cwat; |
| //am_rhow = rhow0 *(1.+ betawat*(t-t0)); |
| tem = tcr - t0; |
| rhow = (b0+(b1+(b2+(b3+(b4+b5*tem)*tem)*tem)*tem)*tem) + (pr1-prif)* |
| (a0+(a1+(a2+(a3+(a4+a5*tem)*tem)*tem)*tem)*tem); |
| //drhow_dpg = 0.0; |
| //drhow_dpc = 0.0; |
| //drhow_dt = 0.0; |
| //d2rhow_d2t = 0.0; |
| //d2rhow_d2pc = 0.0; |
| //d2rhow_dtdpc = 0.0; |
| } |
| rhow = rhow0;//testing??!! |
| 2207,16 → 2383,14 |
| betas = betas0; |
| switch (model_type){ |
| case lewis_and_schrefler3:{//Lewis and Schrefler's book |
| switch (betas_type){ |
| case 0:{ //constant |
| betas = betas0; |
| break; |
| } |
| case lewis_and_schrefler3_mefel:{//from mefel; |
| case 1:{ //moisture dependent - not finished |
| break; |
| } |
| case lewis_and_schrefler3_2:{//Lewis and Schrefler's book |
| break; |
| } |
| default:{ |
| print_err ("unknown model type is required", __FILE__, __LINE__, __func__); |
| abort(); |
| 2228,6 → 2402,30 |
| /** |
| function returns Young's modulus of solid [Pa] |
| @retval e - Young's modulus of solid [Pa] |
| */ |
| double con_hawf3matc::get_e(long /*ipp*/) |
| { |
| return(emod); //provisionally |
| } |
| /** |
| function returns Poisson's ratio [-] |
| @retval nu - Poisson's ratio [-] |
| */ |
| double con_hawf3matc::get_nu(long /*ipp*/) |
| { |
| return(nu); //provisionally |
| } |
| /** |
| function computes enthalpy of evaporation (latent heat of vaporization) |
| @param t - temperature |
| 2251,28 → 2449,3 |
| return(dhvap); |
| } |
| /** |
| function returns Young's modulus of solid [Pa] |
| @retval e - Young's modulus of solid [Pa] |
| */ |
| double con_hawf3matc::get_e(long /*ipp*/) |
| { |
| return(emod); //provisionally |
| } |
| /** |
| function returns Poisson's ratio [-] |
| @retval nu - Poisson's ratio [-] |
| */ |
| double con_hawf3matc::get_nu(long /*ipp*/) |
| { |
| return(nu); //provisionally |
| } |
| /trunk/SIFEL/METR/SRC/consol_hawf3c.h |
|---|
| 2,7 → 2,14 |
| #define CONSOLHAWF3MATC_H |
| #include "genfile.h" |
| #include "baroghel_reten.h" |
| #include "bazant_reten.h" |
| #include "lewis_schrefler.h" |
| #include "gardner.h" |
| #include "potts.h" |
| #include "van_genuchten.h" |
| #include "masin_reten.h" |
| #include "febex_granit_reten.h" |
| class con_hawf3matc |
| { |
| 78,6 → 85,7 |
| double get_fut2(double pw,double pg,double t,long ipp); |
| double get_sw(double pw, double pg, double t,long ipp); |
| double get_xi(double pw, double pg, double t, long ipp); |
| double get_alpha(double pw, double pg, double t,long ipp); |
| double get_rhogw(double pw,double pg,double t); |
| 89,18 → 97,35 |
| double get_rhoga(double pw,double pg,double t); |
| double get_rhog(double pw,double pg,double t); |
| double get_betas(long ipp); |
| double get_dhvap(double t); |
| double get_e(long ipp); |
| double get_nu(long ipp); |
| double get_dhvap(double t); |
| /// Baroghel retention curve: |
| baroghel_reten baroghel_ret; |
| /// Bazant retention curve: |
| bazant_reten bazant_ret; |
| /// 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; |
| /// Masin's retention curve for bentonite |
| masin_reten masin_ret; |
| /// general function for retention curve given by set of data |
| gfunct data; |
| /// FEBEX retention curve |
| febex_granit_reten febex_granit_ret; |
| private: |
| heatairwaterflowmechtype model_type; |
| int compress; |
| int compress, vol_strain_effectc, wrc_vol_strain_effectc, pore_press_effectc, temper_effectc, sr_type, xi_type, betas_type; |
| double ma;//molar mass of dry air |
| double mw;//molar mass of water |
| 109,21 → 134,14 |
| double t0;//[K] reference temperature |
| double p0;//Pa reference athmospheric pressure |
| // PHYSICAL PROPERTIES OF DRY AIR |
| double muga0,alphaa,betaa; |
| // PHYSICAL PROPERTIES OF WATER VAPOUR |
| // from D.Gawin, F.Pesavento (PRVAP.f90) |
| double dv0,bv; |
| double c8,c9,c10,c11,c12,c13; |
| double mugw0,alphaw; |
| // PHYSICAL PROPERTIES OF WATER |
| // from Dariusz Gawin (WATPROP.f90) |
| double rhow0;//density at refernce temperature and pressure |
| double tcr;//[K] critical temperature of water |
| double cwat; |
| double betawat; |
| double hvap0;//latent heat of vaporization at reference temperature |
| double a0,a1,a2; |
| double a3,a4,a5; |
| 130,33 → 148,16 |
| double b0,b1,b2; |
| double b3,b4,b5; |
| double pr1,prif; |
| double muw0;//reference dynamic viscosity of water |
| double conb; |
| double conc; |
| double cpw;//water specific heat at const. pressure = 4181 J/kg.K |
| double lambdaw;//water heat conductivity = 0.6 W/(m.K) |
| double kw0;//compresibility coefficient of water |
| //PHYSICAL PROPERTIES OF SOIL |
| double alpha0; //initial Boit's constant |
| double ks0; //inital bulk modulus of solid phase |
| double kt0; //inital bulk modulus of porous medium |
| double phi0; //inital porosity |
| double kintr0; //intial intrinsic permeability |
| double betas0; //inital cubic thermal dilatation coefficient |
| double deff0; //initial effective diffusivity |
| double betas_dry;//cubic thermal dilatation coefficient of dry soil |
| double betas_wet;//cubic thermal dilatation coefficient of saturated soil |
| double rhos0; //inital volume density of soil skeleton |
| double cps0; //inital specific heat of soil skeleton |
| double lambda0;//effective thermal conductivity of soil skeleton |
| double lambda_dry;//effective thermal conductivity of dry soil |
| double lambda_wet;//effective thermal conductivity of saturated soil |
| double sirr0; //irreversible (residual) saturation degree |
| double tau0; //tortuosity factor 0.4 to 0.6 |
| double gamma,lambda0,s_entry; |
| double emod,nu; |
| double mefel_units; //basic units for pressures = Pa (Pascals) |
| }; |
| #endif |
| /trunk/SIFEL/METR/SRC/consol_wf1c.cpp |
|---|
| 25,6 → 25,12 |
| { |
| alpha = 1.0;//Biot's constant [-] alpha = 1 - kt/ks |
| // |
| vol_strain_effectc = 0; //volumetric strain rate influence: 0=no; 1=yes |
| // |
| wrc_vol_strain_effectc = 0; //volumetric strain rate influence on water retention curve: 0=no; 1=yes |
| // |
| pore_press_effectc = 0; //pore pressure effect acording to effective stress concept: 0=no; 1=yes |
| // |
| rhos0 = 2000.0;//initial solid grain density [kg/m^3] |
| // |
| phi0 = 0.297;//initial porosity [-] |
| 49,7 → 55,7 |
| { |
| xfscanf (in,"%k%m","waterflowmechtype",&waterflowmechtype_kwdset, &model_type); //water flow model type |
| xfscanf (in,"%lf %lf %lf %lf %d %d", &alpha, &phi0, &rhos0, &t0, &sr_type, &xi_type); |
| xfscanf (in,"%lf %d %d %d %lf %lf %lf %d %d", &alpha, &vol_strain_effectc, &pore_press_effectc, &wrc_vol_strain_effectc, &phi0, &rhos0, &t0, &sr_type, &xi_type); |
| switch (model_type){ |
| case lewis_and_schrefler_coup: |
| 1063,7 → 1069,7 |
| 12/9/2008, TKr |
| */ |
| double con_wf1matc::get_phi(double /*pw*/,long ipp) |
| double con_wf1matc::get_phi(double pw,long ipp) |
| { |
| double phi; |
| 1118,7 → 1124,14 |
| 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 |
| kuw = 0.0;// pore pressure effect included in MEFEL |
| if(pore_press_effectc == 1){ |
| // this below in not completed yet |
| chi = 0.0; |
| //chi = Tm->givenontransq(eff_stress_param, ipp); //actual effective stress parameter from MEFEL |
| kuw = -chi;// pore pressure effect included in METR |
| } |
| break; |
| } |
| default:{ |
| 1180,7 → 1193,7 |
| */ |
| double con_wf1matc::get_capwu(double pw,long ipp) |
| { |
| double sw,capwu; |
| double sw,capwu,n,dsr_depsv; |
| capwu = 0.0; |
| 1191,7 → 1204,17 |
| 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 |
| capwu = 0.0;// volumetric strain effect included in TRFEL |
| if(vol_strain_effectc == 1){ |
| sw = get_sw(pw,ipp); |
| capwu = alpha*sw;//volumetric strain effect included in METR |
| if(wrc_vol_strain_effectc == 1){ |
| // this must be rewritten for METR elements |
| dsr_depsv = Tm->givenontransq(der_saturation_deg_depsv, ipp); //actual derivative of saturation degree with respect to volumetric strain |
| n = get_phi(pw,ipp); |
| capwu = capwu + n*dsr_depsv;//volumetric strain effect on retention curve |
| } |
| } |
| break; |
| } |
| default:{ |
| /trunk/SIFEL/METR/SRC/consol_wf1c.h |
|---|
| 91,7 → 91,7 |
| private: |
| waterflowmechtype model_type; |
| int sr_type,xi_type; |
| int vol_strain_effectc,pore_press_effectc,wrc_vol_strain_effectc,sr_type,xi_type; |
| double alpha,phi0,rhos0,rhow0,t0; |
| double gamma,lambda0,s_entry; |
| }; |
| /trunk/SIFEL/METR/SRC/cpcsolver.cpp |
|---|
| 203,7 → 203,7 |
| solver_restore (r,fp,ii,Mp->time,dt,&Mp->timecon,Ndofm); |
| // Mt->rhs_save (fp); |
| // initiation of mechanical material models |
| Mm->initmaterialmodels(lcid); |
| Mm->initmaterialmodels(lcid, true); |
| print_init(-1, "at"); |
| // print_step(lcid,i,Mp->time,f); |
| // print_flush(); |
| 211,7 → 211,7 |
| else |
| { |
| // initiation of mechanical material models |
| Mm->initmaterialmodels(lcid); |
| Mm->initmaterialmodels(lcid, false); |
| // compute stresses from eigenstrains or initial irreversible strains |
| internal_forces(lcid, fp); |
| print_init(-1, "wt"); |
| 876,7 → 876,7 |
| // update status indicators of prescribed displacements on elements |
| Mt->elemprescdisp(); |
| // initiation of mechanical material models |
| Mm->initmaterialmodels(lcid); |
| Mm->initmaterialmodels(lcid, true); |
| print_init(-1, "at"); |
| // print_step(lcid,i,Mp->time,f); |
| 885,7 → 885,7 |
| else |
| { |
| // initiation of mechanical material models |
| Mm->initmaterialmodels(lcid); |
| Mm->initmaterialmodels(lcid, false); |
| // compute stresses from eigenstrains or initial irreversible strains |
| internal_forces(lcid, fp); |
| print_init(-1, "wt"); |
| 1231,7 → 1231,7 |
| Mt->restore_nodval(r, Ndofm); |
| // initiation of mechanical material models |
| Mm->initmaterialmodels(lcid); |
| Mm->initmaterialmodels(lcid, false); |
| // |
| /trunk/SIFEL/METR/SRC/fcsolver.cpp |
|---|
| 230,8 → 230,8 |
| // inital printing of output and graphical informations |
| if (Mp->hdbcont.restore_stat()){ |
| init_trfel_mefel(); |
| Mm->initmaterialmodels(lcid); |
| Mm->aip_initmaterialmodels(lcid); |
| Mm->initmaterialmodels(lcid, true); |
| Mm->aip_initmaterialmodels(lcid, true); |
| if (Mespr==1) |
| fprintf (stdout,"\n Reading of backup file for MEFEL\n"); |
| solver_restore (lhs,rhs,i,Cp->time,dt,&Cp->timecon,n); |
| 246,8 → 246,8 |
| // approximation of temperatures and humidities from TRFEL nodes to MEFEL integration points |
| init_trfel_mefel (); |
| //init_trfel_mefel_by_nodes_comp (); |
| Mm->initmaterialmodels(lcid); |
| Mm->aip_initmaterialmodels(lcid); |
| Mm->initmaterialmodels(lcid, false); |
| Mm->aip_initmaterialmodels(lcid, false); |
| // compute stresses from eigenstrains or initial irreversible strains |
| if (Mp->eigstrains > 0) |
| internal_forces(lcid, rhs); |
| 564,8 → 564,8 |
| // passes initial values from TRFEL AUX. points to nonmechq array in MEFEL (associated with the regular MEFEL INT. points) |
| // initial values of TRFEL AUX. points are calculated in the body of init_trfel_mefel() by aip_initapproximation |
| init_trfel_mefel(); |
| Mm->initmaterialmodels(lcid); |
| Mm->aip_initmaterialmodels(lcid); |
| Mm->initmaterialmodels(lcid, true); |
| Mm->aip_initmaterialmodels(lcid, true); |
| if (Mespr==1) |
| fprintf (stdout,"\n Reading of backup file for MEFEL\n"); |
| solver_restore (lhs,rhs,i,Cp->time,dt,&Cp->timecon,n); |
| 582,8 → 582,8 |
| // passes initial values from TRFEL AUX. points to nonmechq array in MEFEL (associated with the regular MEFEL INT. points) |
| // initial values of TRFEL AUX. points are calculated in the body of init_trfel_mefel() by aip_initapproximation |
| init_trfel_mefel (); |
| Mm->initmaterialmodels(lcid); |
| Mm->aip_initmaterialmodels(lcid); |
| Mm->initmaterialmodels(lcid, false); |
| Mm->aip_initmaterialmodels(lcid, false); |
| // compute stresses from eigenstrains or initial irreversible strains |
| if ((Mp->eigstrains > 0) && (Mp->eigstrains < 4)) |
| internal_forces(lcid, rhs); |
| 1157,8 → 1157,8 |
| // inital printing of output and graphical informations |
| if (Mp->hdbcont.restore_stat()){ |
| init_trfel_mefel(); |
| Mm->initmaterialmodels(lcid); |
| Mm->aip_initmaterialmodels(lcid); |
| Mm->initmaterialmodels(lcid, true); |
| Mm->aip_initmaterialmodels(lcid, true); |
| if (Mespr==1) |
| fprintf (stdout,"\n Reading of backup file for MEFEL\n"); |
| solver_restore (lhs,rhs,i,Cp->time,dt,&Cp->timecon,n); |
| 1173,8 → 1173,8 |
| // approximation of temperatures and humidities from TRFEL nodes to MEFEL integration points |
| //init_trfel_mefel (); |
| init_trfel_mefel_by_nodes_comp (); |
| Mm->initmaterialmodels(lcid); |
| Mm->aip_initmaterialmodels(lcid); |
| Mm->initmaterialmodels(lcid, false); |
| Mm->aip_initmaterialmodels(lcid, false); |
| // compute stresses from eigenstrains or initial irreversible strains |
| if (Mp->eigstrains > 0) |
| internal_forces(lcid, rhs); |
| /trunk/SIFEL/METR/SRC/pcsolver.cpp |
|---|
| 246,9 → 246,10 |
| { |
| if (tret >= 0){ |
| // equilibrium has been attained |
| // approximation of nodal values into integration points |
| approximation (); |
| if (Cp->dpt == pass_by_aux_ip) |
| aip_approximation(Mm->tnip, MTipmap); |
| Tm->updateipval(); |
| 1146,7 → 1147,7 |
| // initiation of mechanical material models |
| // approximation of temperatures and humidities from TRFEL nodes to MEFEL integration points |
| init_trfel_mefel (); |
| Mm->initmaterialmodels(lcid); |
| Mm->initmaterialmodels(lcid, true); |
| if (Mespr==1) |
| fprintf (stdout,"\n Reading of backup file for MEFEL\n"); |
| solver_restore (r,fp,ii,Mp->time,dt,&Mp->timecon,Ndofm); |
| 1159,7 → 1160,7 |
| // initiation of mechanical material models |
| // approximation of temperatures and humidities from TRFEL nodes to MEFEL integration points |
| init_trfel_mefel (); |
| Mm->initmaterialmodels(lcid); |
| Mm->initmaterialmodels(lcid, false); |
| // compute stresses from eigenstrains or initial irreversible strains |
| if (Mp->eigstrains > 0) |
| internal_forces(lcid, fp); |
| 1683,7 → 1684,7 |
| fprintf (stdout,"\n Reading of backup file for MEFEL\n"); |
| solver_restore (r,fp,ii,Mp->time,dt,&Mp->timecon,Ndofm); |
| // initiation of mechanical material models |
| Mm->initmaterialmodels(lcid); |
| Mm->initmaterialmodels(lcid, true); |
| print_init(-1, "at"); |
| // print_step(lcid,i,Mp->time,f); |
| // print_flush(); |
| 1691,7 → 1692,7 |
| else |
| { |
| // initiation of mechanical material models |
| Mm->initmaterialmodels(lcid); |
| Mm->initmaterialmodels(lcid, false); |
| // compute stresses from eigenstrains or initial irreversible strains |
| internal_forces(lcid, fp); |
| print_init(-1, "wt"); |
| 2290,7 → 2291,7 |
| approximation_humid (); |
| // initiation of mechanical material models |
| Mm->initmaterialmodels(lcid); |
| Mm->initmaterialmodels(lcid, false); |
| // ********************************************************* |
| // transport part |
| 2601,7 → 2602,7 |
| approximation_humid (); |
| // initiation of mechanical material models |
| Mm->initmaterialmodels(lcid); |
| Mm->initmaterialmodels(lcid, false); |
| // ********************************************************* |
| // transport part |
| /trunk/SIFEL/PARMEF/SRC/pcpsolver.cpp |
|---|
| 100,7 → 100,7 |
| solver_restore (r,fp,i,Mp->time,dt,&Mp->timecon,n); |
| Mt->rhs_save (fp); |
| // initiation of mechanical material models |
| Mm->initmaterialmodels(lcid); |
| Mm->initmaterialmodels(lcid, true); |
| print_init(-1, "at",Pmp->fni,Pmp->fei); |
| //print_step(lcid,i,Mp->time,f); |
| //print_flush(); |
| 108,7 → 108,7 |
| else |
| { |
| // initiation of mechanical material models |
| Mm->initmaterialmodels(lcid); |
| Mm->initmaterialmodels(lcid, false); |
| print_init(-1, "wt",Pmp->fni,Pmp->fei); |
| print_step(lcid, i, Mp->time, f); |
| print_flush(); |
| /trunk/SIFEL/PARMEF/SRC/pmtsolver.cpp |
|---|
| 217,9 → 217,9 |
| if (Mp->hdbcont.restore_stat()) |
| { |
| // initiation of mechanical material models |
| Mm->initmaterialmodels(lcid); |
| Mm->initmaterialmodels(lcid, true); |
| // initiation of material models on auxiliary integration points |
| Mm->aip_initmaterialmodels(lcid); |
| Mm->aip_initmaterialmodels(lcid, true); |
| if (Mespr==1) |
| fprintf (stdout,"\n Reading of backup file\n"); |
| solver_restore (mt_gv.r, mt_gv.fp, mt_gv.istep, Mp->time, dt, &Mp->timecon, n); |
| 233,7 → 233,7 |
| update_elnod_stat_after_hdbrest(lcid); |
| // update actual material index for all elements |
| Mm->update_actual_mat_id(lcid, 0); |
| Mm->update_actual_mat_id(lcid, 0, true); |
| // set indicator of calculation started from restored backup data |
| rest_calc = 1; |
| print_init(-1, "at",Pmp->fni,Pmp->fei); |
| 242,9 → 242,9 |
| { |
| mt_gv.istep=0; |
| // initiation of mechanical material models |
| Mm->initmaterialmodels(lcid); |
| Mm->initmaterialmodels(lcid, false); |
| // initiation of material models on auxiliary integration points |
| Mm->aip_initmaterialmodels(lcid); |
| Mm->aip_initmaterialmodels(lcid, false); |
| // compute stresses from eigenstrains or initial irreversible strains |
| // compute stresses from eigenstrains or initial irreversible strains |
| if ((Mp->eigstrains > 0) && (Mp->eigstrains < 4)) |
| 330,7 → 330,7 |
| // update actual material index for all old elements and call initval for those whose indeces have changed |
| if (rest_calc == 0) |
| { |
| um = Mm->update_actual_mat_id(lcid, 1); |
| um = Mm->update_actual_mat_id(lcid, 1, false); |
| if (um) fusm = yes; |
| else fusm = no; |
| } |
| 499,7 → 499,7 |
| fprintf (stdout,"\n Reading of backup file\n"); |
| solver_restore (r,fp,i,Mp->time,dt,&Mp->timecon,n); |
| // initiation of mechanical material models |
| Mm->initmaterialmodels(lcid); |
| Mm->initmaterialmodels(lcid, true); |
| print_init(-1, "at",Pmp->fni,Pmp->fei); |
| print_step(lcid, i, Mp->time, f); |
| print_flush(); |
| 507,7 → 507,7 |
| else |
| { |
| // initiation of mechanical material models |
| Mm->initmaterialmodels(lcid); |
| Mm->initmaterialmodels(lcid, false); |
| print_init(-1, "wt",Pmp->fni,Pmp->fei); |
| print_step(lcid, i, Mp->time, f); |
| print_flush(); |
| /trunk/SIFEL/PARMETR/SRC/pcpcsolver.cpp |
|---|
| 217,7 → 217,7 |
| // approximation of humidities from TRFEL nodes to MEFEL integration points |
| //approximation_inithumid (); |
| // initiation of mechanical material models |
| Mm->initmaterialmodels(lcid); |
| Mm->initmaterialmodels(lcid, true); |
| print_init(-1, "at",Pcp->fnim,Pcp->feim); |
| // print_step(lcid,i,Mp->time,f); |
| // print_flush(); |
| 229,7 → 229,7 |
| // approximation of humidities from TRFEL nodes to MEFEL integration points |
| //approximation_inithumid (); |
| // initiation of mechanical material models |
| Mm->initmaterialmodels(lcid); |
| Mm->initmaterialmodels(lcid, false); |
| print_init(-1, "wt",Pcp->fnim,Pcp->feim); |
| print_step(lcid, i, Mp->time, f); |
| print_flush(); |
| 984,7 → 984,7 |
| solver_restore (r,fp,ii,Mp->time,dt,&Mp->timecon,Ndofm); |
| Mt->rhs_save (fp); |
| // initiation of mechanical material models |
| Mm->initmaterialmodels(lcid); |
| Mm->initmaterialmodels(lcid, true); |
| print_init(-1, "at",Pcp->fnim,Pcp->feim); |
| // print_step(lcid,i,Mp->time,f); |
| // print_flush(); |
| 992,7 → 992,7 |
| else |
| { |
| // initiation of mechanical material models |
| Mm->initmaterialmodels(lcid); |
| Mm->initmaterialmodels(lcid, false); |
| print_init(-1, "wt",Pcp->fnim,Pcp->feim); |
| print_step(lcid, i, Mp->time, f); |
| print_flush(); |
| /trunk/SIFEL/PARMETR/SRC/ppcsolver.cpp |
|---|
| 453,7 → 453,7 |
| fprintf (stdout,"\n Reading of backup file\n"); |
| solver_restore (mt_gv.r, mt_gv.fp, mt_gv.istep, Mp->time, dt, &Mp->timecon, n); |
| // initiation of mechanical material models |
| Mm->initmaterialmodels(lcid); |
| Mm->initmaterialmodels(lcid, true); |
| print_init(-1, "at",Pcp->fnim,Pcp->feim); |
| //print_step(lcid,i,Mp->time,f); |
| print_flush(); |
| 462,7 → 462,7 |
| { |
| mt_gv.istep=0; |
| // initiation of mechanical material models |
| Mm->initmaterialmodels(lcid); |
| Mm->initmaterialmodels(lcid, false); |
| // compute stresses from eigenstrains or initial irreversible strains |
| if (Mp->eigstrains > 0) |
| internal_forces(lcid, mt_gv.fp); |
| 1224,7 → 1224,7 |
| fprintf (stdout,"\n Reading of backup file for MEFEL\n"); |
| solver_restore (r,fp,ii,Mp->time,dt,&Mp->timecon,Ndofm); |
| // initiation of mechanical material models |
| Mm->initmaterialmodels(lcid); |
| Mm->initmaterialmodels(lcid, true); |
| print_init(-1, "at",Pcp->fnim,Pcp->feim); |
| // print_step(lcid,i,Mp->time,f); |
| // print_flush(); |
| 1232,7 → 1232,7 |
| else |
| { |
| // initiation of mechanical material models |
| Mm->initmaterialmodels(lcid); |
| Mm->initmaterialmodels(lcid, false); |
| print_init(-1, "wt",Pcp->fnim,Pcp->feim); |
| print_step(lcid, i, Mp->time, f); |
| print_flush(); |
| 1785,7 → 1785,7 |
| fprintf (stdout,"\n Reading of backup file for MEFEL\n"); |
| solver_restore (r,fp,ii,Mp->time,dt,&Mp->timecon,Ndofm); |
| // initiation of mechanical material models |
| Mm->initmaterialmodels(lcid); |
| Mm->initmaterialmodels(lcid, true); |
| print_init(-1, "at",Pcp->fnim,Pcp->feim); |
| // print_step(lcid,i,Mp->time,f); |
| // print_flush(); |
| 1793,7 → 1793,7 |
| else |
| { |
| // initiation of mechanical material models |
| Mm->initmaterialmodels(lcid); |
| Mm->initmaterialmodels(lcid, false); |
| print_init(-1, "wt",Pcp->fnim,Pcp->feim); |
| print_step(lcid, i, Mp->time, f); |
| print_flush(); |
| /trunk/SIFEL/TRFEL/SRC/aliast.h |
|---|
| 200,8 → 200,8 |
| const kwdset sr_type_kwdset(sizeof(sr_typestr)/sizeof(*sr_typestr),sr_typestr); |
| // aliases for models of effective stress factor (parameter) type |
| enum xi_type {biot_xi=1,masin_xi=2,biot_reduced_xi=3,biot_masin_xi=4}; |
| const enumstr xi_typestr[] = {{"biot_xi",1}, {"masin_xi",2},{"biot_reduced_xi",3},{"biot_masin_xi",4}}; |
| enum xi_type {biot_xi=1,masin_xi=2,biot_reduced_xi=3,biot_masin_xi=4,masin_hypopl_xi=5}; |
| const enumstr xi_typestr[] = {{"biot_xi",1}, {"masin_xi",2},{"biot_reduced_xi",3},{"biot_masin_xi",4},{"masin_hypopl_xi",5}}; |
| const kwdset xi_type_kwdset(sizeof(xi_typestr)/sizeof(*xi_typestr),xi_typestr); |
| /trunk/SIFEL/TRFEL/SRC/consol_hawf3.cpp |
|---|
| 30,6 → 30,8 |
| #include "bazant_reten.h" |
| #include "globmatt.h" |
| #include <errno.h> |
| con_hawf3mat::con_hawf3mat() |
| { |
| compress = 0; //compressible grains: 0=no; 1=yes |
| 42,10 → 44,12 |
| xi_type = 1; //effective stress parameter type |
| lambda_type =0; //heat conduction calculation type |
| cps_type = 0; //specific heat calculation type |
| thermal_capacity_type = 0; //thermal capacity type |
| betas_type = 0; //thermal expansion calculation type |
| vol_strain_effect = 0; //volumetric strain rate influence: 0=no; 1=yes |
| mefel_units = 1.0; //basic units for pressures = Pa (Pascals) |
| wrc_vol_strain_effect = 0; //volumetric strain rate influence on water retention curve: 0=no; 1=yes |
| //STATE VARIABLES |
| 102,8 → 106,8 |
| b3 = -9.2118e-5; b4 = 3.3534e-7; b5 = -4.4034e-10; |
| pr1 = 1.0e7; prif = 2.0e7; |
| //muw0 = 0.6612;//testing??!! |
| muw0 = 1.0e-3; |
| muw0 = 0.6612; |
| muw0const = 1.0e-3; |
| conb = 229.0; |
| conc = -1.562; |
| cpw0 = 4181.0; |
| 117,7 → 121,11 |
| kt0 = 0.0; |
| phi0 = 0.0; |
| kintr0 = 0.0; |
| kintrw0 = 0.0; |
| kintrg0 = 0.0; |
| betas0 = 0.0; |
| betas_dry = 0.0; |
| betas_wet = 0.0; |
| deff0 = 0.0; |
| cdiff0 = 25.4e-6;//initial water vapour diffusivity in air |
| rhos0 = 0.0; |
| 127,6 → 135,8 |
| lambda_wet = 0.0; |
| sr_dry = 0.0; |
| sr_wet = 0.0; |
| rhocp_dry = 0.0; |
| rhocp_wet = 0.0; |
| tau0 = 0.0; |
| gamma = 0.0; |
| 142,6 → 152,7 |
| bb1 = 0.0; |
| phi01 = 0.0; |
| ng = 0.0; |
| //gas relative permeability parameters: |
| krg0 = 1.0; |
| 181,7 → 192,7 |
| xfscanf (in,"%d", &compress); |
| // common material parameters |
| xfscanf (in,"%le %le %le %le %le %le %d %d %d %d %d %d %d %d %d %d", &alpha0, &ks0, &rhos0, &pw_bc, &tau0, &kintr0, &por_type, &kintr_type, &krw_type, &krg_type, &deff_type, &sr_type, &xi_type, &lambda_type, &cps_type, &betas_type); |
| xfscanf (in,"%le %le %le %le %le %le %d %d %d %d %d %d %d %d %d %d %d", &alpha0, &ks0, &rhos0, &pw_bc, &tau0, &kintr0, &por_type, &kintr_type, &krw_type, &krg_type, &deff_type, &sr_type, &xi_type, &lambda_type, &cps_type, &thermal_capacity_type, &betas_type); |
| switch (model_type){ |
| case artificial3:{//artificial isotropic material for non-isotherma air-water flow |
| 192,7 → 203,7 |
| break; |
| } |
| case lewis_and_schrefler3_mefel:{//Lewis and Schrefler's model approach coupled with mechanics |
| xfscanf (in,"%le %d", &mefel_units, &vol_strain_effect); |
| xfscanf (in,"%le %d %d", &mefel_units, &vol_strain_effect, &wrc_vol_strain_effect); |
| break; |
| } |
| default:{ |
| 219,17 → 230,27 |
| //intrinsic permebility calculation: |
| switch (kintr_type){ |
| case 0:{//constant |
| case 0:{//constant - for permeability of both phases |
| break; |
| } |
| case 1:{//dependent on porosity |
| case 1:{//dependent on porosity - for permeability of both phases |
| xfscanf (in,"%le %le", &bb1, &phi01); |
| break; |
| } |
| case 2:{//dependent on porosity - cubic and quadratic |
| case 2:{//dependent on porosity - cubic and quadratic - for permeability of both phases |
| xfscanf (in,"%le", &phi01); |
| break; |
| } |
| case 3:{//constant separate permeabilities of liquid water and gas |
| xfscanf (in,"%le", &kintrw0); |
| xfscanf (in,"%le", &kintrg0); |
| break; |
| } |
| case 4:{//dependent on porosity - separate permeabilities of liquid water and gas |
| xfscanf (in,"%le %le %le", &kintrw0, &bb1, &phi01); |
| xfscanf (in,"%le %le", &kintrg0, &ng); |
| break; |
| } |
| default:{ |
| print_err ("unknown model type is required", __FILE__, __LINE__, __func__); |
| abort(); |
| 375,6 → 396,10 |
| xfscanf (in,"%le %le %le %le", &lambda_dry, &lambda_wet, &sr_dry, &sr_wet); |
| break; |
| } |
| case 3:{//dependent on moisture |
| xfscanf (in,"%le %le", &lambda_dry, &lambda_wet); |
| break; |
| } |
| default:{ |
| print_err ("unknown model type is required", __FILE__, __LINE__, __func__); |
| abort(); |
| 388,10 → 413,14 |
| break; |
| } |
| case 1:{//dependent on moisture |
| xfscanf (in,"%le %le %le %le", &cps_dry, &cps_wet, &sr_dry, &sr_wet); |
| break; |
| } |
| case 2:{//dependent on moisture |
| xfscanf (in,"%le %le", &cps_dry, &cps_wet); |
| break; |
| } |
| case 2:{//dependent on temperature |
| case 3:{//dependent on temperature |
| xfscanf (in,"%le %le", &cps0, &cps_lin); |
| break; |
| } |
| 401,6 → 430,25 |
| } |
| } |
| //thermal capacity calculation: |
| switch (thermal_capacity_type){ |
| case 0: |
| case 1: |
| case 2:{ |
| //nothing to read |
| break; |
| } |
| case 3:{//effective heat capacity is directly measured with respec to saturation degree |
| //linear approximation |
| xfscanf (in,"%le %le %le %le", &rhocp_dry, &rhocp_wet, &sr_dry, &sr_wet); |
| break; |
| } |
| default:{ |
| print_err ("unknown model type is required", __FILE__, __LINE__, __func__); |
| abort(); |
| } |
| } |
| //thermal expansion calculation: |
| switch (betas_type){ |
| case 0:{//constant |
| 521,14 → 569,17 |
| cps = cps0; |
| break; |
| } |
| case 1: |
| case 3:{ //moisture dependent |
| case 1:{ //moisture dependent |
| sw = get_sw(pw,pg,t,ipp); |
| cps = cps_dry + (cps_wet - cps_dry)*(sw - sr_dry)/(sr_wet - sr_dry); |
| break; |
| } |
| case 2:{ //moisture dependent |
| sw = get_sw(pw,pg,t,ipp); |
| cps = cps_dry + (cps_wet - cps_dry)*sw; //this is maybe linear approximation for bentonite: |
| break; |
| } |
| case 2: |
| case 4:{ //temperature dependent |
| case 3:{ //temperature dependent |
| cps = cps_lin*t + cps_0; //this is maybe linear approximation for bentonite: |
| break; |
| } |
| 760,7 → 811,7 |
| } |
| case van_genuchten_sr:{//partially saturated medium = Van Genuchten model |
| //pg = pg - p0;//pore gas pressure without atmospheric pressure//debug?? |
| pg = pg - p0;//pore gas pressure without atmospheric pressure//debug?? |
| pc = pg - pw; |
| sw = van_genuchten_ret.sw(-pc,t);//suction = capillary pressure |
| 849,7 → 900,7 |
| } |
| case van_genuchten_sr:{//partially saturated medium = Van Genuchten model |
| //pg = pg - p0;//pore gas pressure without atmospheric pressure//debug?? |
| pg = pg - p0;//pore gas pressure without atmospheric pressure//debug?? |
| pc = pg - pw; |
| dsw_dpc = -van_genuchten_ret.dsw_dpw(-pc,t);//suction = capillary pressure |
| 938,7 → 989,7 |
| } |
| case van_genuchten_sr:{//partially saturated medium = Van Genuchten model |
| //pg = pg - p0;//pore gas pressure without atmospheric pressure//debug?? |
| pg = pg - p0;//pore gas pressure without atmospheric pressure//debug?? |
| pc = pg - pw; |
| dsw_dt = van_genuchten_ret.dsw_dt(-pc,t); |
| 1062,7 → 1113,103 |
| return(kintr); |
| } |
| /** |
| function computes intrinsic permeability assumed only for liquid (some models have two separate values) |
| @param pw - water pressure |
| @param pg - air pressure |
| @param t - temperature |
| @param ipp - number of integration point |
| @retval kintr - intrinsic permeability of water |
| 02/02/2024, TKr |
| */ |
| double con_hawf3mat::get_kintrw(double pw, double pg, double t, long ipp) |
| { |
| double kintr; |
| double phi; |
| kintr = kintrw0; |
| switch (kintr_type){ |
| case 0: |
| case 3:{//constant |
| kintr = kintrw0; |
| break; |
| } |
| case 1: |
| case 4:{//dependent on porosity |
| phi = get_porosity(ipp); |
| kintr = kintrw0*exp(bb1*(phi - phi01)); |
| break; |
| } |
| case 2:{//dependent on porosity - cubic and quadratic |
| //Kozeny's approach for bentonite: |
| phi = get_porosity(ipp); |
| kintr = kintrw0*phi*phi*phi/(1 - phi01)*(1 - phi01)*(1 - phi)*(1 - phi)/(phi01*phi01*phi01); |
| break; |
| } |
| default:{ |
| print_err ("unknown model type is required", __FILE__, __LINE__, __func__); |
| abort(); |
| } |
| } |
| return(kintr); |
| } |
| /** |
| function computes intrinsic permeability assumed only for gas (some models have two separate values) |
| @param pw - water pressure |
| @param pg - air pressure |
| @param t - temperature |
| @param ipp - number of integration point |
| @retval kintr - intrinsic permeability of gas |
| 02/02/2024, TKr |
| */ |
| double con_hawf3mat::get_kintrg(double pw, double pg, double t, long ipp) |
| { |
| double kintr; |
| double phi,e; |
| kintr = kintrg0; |
| switch (kintr_type){ |
| case 0: |
| case 3:{//constant |
| kintr = kintrg0; |
| break; |
| } |
| case 1: |
| case 4:{//dependent on porosity |
| phi = get_porosity(ipp); |
| e = phi/(1 - phi); |
| kintr = kintrg0*pow(e,ng); |
| break; |
| } |
| case 2:{//dependent on porosity - cubic and quadratic |
| //Kozeny's approach for bentonite: |
| //phi = get_porosity(ipp); |
| //kintr = kintrw0*phi*phi*phi/(1 - phi01)*(1 - phi01)*(1 - phi)*(1 - phi)/(phi01*phi01*phi01); |
| kintr = kintrg0;//temp. |
| break; |
| } |
| default:{ |
| print_err ("unknown model type is required", __FILE__, __LINE__, __func__); |
| abort(); |
| } |
| } |
| return(kintr); |
| } |
| /** |
| function computes water relative permeability |
| @param pw - water pressure |
| @param pg - air pressure |
| 1311,7 → 1458,7 |
| */ |
| double con_hawf3mat::get_pgw(double pw,double pg,double t) |
| { |
| double pgw,rhow,pgws,tt,pc; |
| double pgw,rhow,pgws,tt,arg,pc; |
| //critical point of water check |
| if(t >= tcr) |
| 1324,13 → 1471,24 |
| if(rel_gas_press == 1) |
| pg = pg + p0; |
| pc = pg - pw; |
| pgw = pgws*exp(-1.0*pc*mw/rhow/gasr/tt); |
| arg = -1.0*pc*mw/rhow/gasr/tt; |
| if(arg < -7){ |
| pgw = pgws*0.001; |
| //printf("relative humidity correction to limit value 0.001",__FILE__,__LINE__,__func__); |
| //fprintf(Outt,"relative humidity correction to limit value 0.001",__FILE__,__LINE__,__func__); |
| //print_err("relative humidity correction to limit value 0.001",__FILE__,__LINE__,__func__); |
| } |
| else{ |
| //pgw = pgws*exp(-1.0*pc*mw/rhow/gasr/tt); |
| pgw = pgws*exp(arg); |
| } |
| check_math_err(); |
| return(pgw); |
| } |
| 1341,29 → 1499,40 |
| @param pw - pore water pressure |
| @param pg - pore gas pressure |
| @param t - temperature |
| @retval dpgw_dpc - partial derivative of pgw with respect to pc (Kelvin equation) |
| */ |
| double con_hawf3mat::get_dpgw_dpc(double pw,double pg,double t) |
| { |
| double dpgw_dpc,pgws,rhow,tt,pc; |
| double dpgw_dpc,pgws,rhow,tt,pc,arg; |
| //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; |
| arg = -1.0*pc*mw/rhow/gasr/tt; |
| if(arg < -7){ |
| dpgw_dpc = 0.0; |
| //printf("relative humidity correction to limit value 0.001",__FILE__,__LINE__,__func__); |
| //fprintf(Outt,"relative humidity correction to limit value 0.001",__FILE__,__LINE__,__func__); |
| //print_err("relative humidity correction to limit value 0.001",__FILE__,__LINE__,__func__); |
| } |
| else{ |
| //dpgw_dpc = -pgws*exp(-pc*mw/rhow/gasr/tt)*mw/rhow/gasr/tt; |
| dpgw_dpc = -pgws*exp(arg)*mw/rhow/gasr/tt; |
| } |
| check_math_err(); |
| return(dpgw_dpc); |
| 1380,7 → 1549,7 |
| */ |
| double con_hawf3mat::get_dpgw_dt(double pw,double pg,double t) |
| { |
| double dpgw_dt,dpgws_dt,tt,pgws,rhow,pc,drhow_dt; |
| double dpgw_dt,dpgws_dt,tt,pgws,rhow,pc,drhow_dt,arg; |
| //critical point of water check |
| if(t >= tcr) |
| 1394,14 → 1563,31 |
| dpgws_dt = get_dpgws_dt(tt); |
| drhow_dt = get_drhow_dt(tt); |
| arg = -1.0*pc*mw/rhow/gasr/tt; |
| if(arg < -7){ |
| dpgw_dt = 0.0; |
| if(t < tcr){ |
| dpgw_dt = dpgws_dt*0.0001 + pgws*0.0001*(7/tt); |
| } |
| else{ |
| dpgw_dt = dpgws_dt*0.0001; |
| } |
| //printf("relative humidity correction to limit value 0.001",__FILE__,__LINE__,__func__); |
| //fprintf(Outt,"relative humidity correction to limit value 0.001",__FILE__,__LINE__,__func__); |
| //print_err("relative humidity correction to limit value 0.001",__FILE__,__LINE__,__func__); |
| } |
| else{ |
| 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); |
| } |
| 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); |
| return(dpgw_dt); |
| } |
| 1698,15 → 1884,15 |
| 03/10/2023 TKr |
| */ |
| double con_hawf3mat::get_muw(double /*t*/) |
| double con_hawf3mat::get_muw(double t) |
| { |
| double muw; |
| //muw = muw0*pow((t - conb),conc); |
| muw = muw0*pow((t - conb),conc); |
| //check_math_errel(0); |
| muw = muw0;//temp. |
| //muw = muw0const;//constant value |
| return(muw); |
| } |
| 1753,16 → 1939,17 |
| 03/10/2023 TKr |
| */ |
| double con_hawf3mat::get_betaw(double /*t*/) |
| double con_hawf3mat::get_betaw(double t) |
| { |
| double betaw; |
| betaw = 0.68e-4;//[K-1] at t=273.15 |
| //betaw = 0.68e-4;//[K-1] at t=273.15 |
| //betaw = 10.1e-4;//[K-1] at t=420.0 |
| betaw = 0.63e-5;//temp. |
| //betaw = 0.68 + (10.1 - 0.68)/(420.0 - 273.15)*(t - 273.15);//linear expression (rough) |
| //betaw = 0.63e-5;//temp. |
| betaw = 0.68 + (10.1 - 0.68)/(420.0 - 273.15)*(t - 273.15);//linear expression (rough) |
| betaw = betaw*1.0e-4; |
| return(betaw); |
| } |
| 1933,9 → 2120,8 |
| */ |
| double con_hawf3mat::get_rhocp(double pw,double pg,double t,long ipp) |
| { |
| double sw,n,rhocp,rhow,rhog,rhogw,cps,cpw,cpga,cpgw,rhos; |
| double sw,n,rhocp,rhow,rhog,rhogw,cps,cpw,cpga,cpgw,rhos,cp; |
| sw = get_sw(pw,pg,t,ipp); |
| n = get_porosity(ipp); |
| 1950,14 → 2136,36 |
| cpgw = get_cpgw(); |
| rhos = get_rhos(t); |
| if((cps_type == 3) || (cps_type == 4)){//effective heat capacity without water and gas components |
| rhocp = rhos*cps; |
| switch (thermal_capacity_type){ |
| case 0:{//constant value without water and gas phases |
| rhocp = (1.0-n)*rhos*cps;// constant effective heat capacity only for solid phase |
| break; |
| } |
| else{ //effective heat capacity from all components |
| rhocp = (1.0-n)*rhos*cps;// solid phase |
| case 1:{//effective heat capacity from all components |
| rhocp = (1.0-n)*rhos*cps;// solid phase = rhod*cps |
| rhocp = rhocp + n*(sw*rhow*cpw + (1.0-sw)*(rhog*cpga + rhogw*(cpgw-cpga)));//liquid and gas phases |
| break; |
| } |
| case 2:{ |
| rhocp = rhos*cps;// special case - constant heat capacity only for solid phase |
| break; |
| } |
| case 3:{//effective heat capacity is directly measured with respec to saturation degree |
| sw = get_sw(pw,pg,t,ipp); |
| //linear approximation |
| rhocp = rhocp_dry + (rhocp_wet - rhocp_dry)*(sw - sr_dry)/(sr_wet - sr_dry); |
| if(rhocp < rhocp_dry) |
| rhocp = rhocp_dry; |
| if(rhocp > rhocp_wet) |
| rhocp = rhocp_wet; |
| break; |
| } |
| default:{ |
| print_err ("unknown model type is required", __FILE__, __LINE__, __func__); |
| abort(); |
| } |
| } |
| //this below is temporarilly for mock_up experiment testing |
| switch (model_type){ |
| 2034,6 → 2242,12 |
| case 2:{// linear approximation |
| sw = get_sw(pw,pg,t,ipp); |
| lambdaeff = lambda_dry + (lambda_wet - lambda_dry)*(sw - sr_dry)/(sr_wet - sr_dry); |
| if(lambdaeff < lambda_dry) |
| lambdaeff = lambda_dry; |
| if(lambdaeff > lambda_wet) |
| lambdaeff = lambda_wet; |
| //fprintf(Outt,"lambda = %e, sr = %e\n",lambdaeff,sw); |
| //fflush(Outt); |
| break; |
| } |
| case 3:{// |
| 2712,8 → 2926,9 |
| alpha = get_alpha(pw,pg,t,ipp); |
| fwu = -1.0*(sg*rhogw + sw*rhow)*alpha*depsv_r;//volumetric strain effect |
| fwu = fwu - n*(rhow - rhogw)*dsr_depsv*depsv_r;//volumetric strain effect on retention curve |
| if(wrc_vol_strain_effect == 1) |
| fwu = fwu - n*(rhow - rhogw)*dsr_depsv*depsv_r;//volumetric strain effect on retention curve |
| } |
| break; |
| } |
| 2817,8 → 3032,9 |
| rhoga = get_rhoga(pw,pg,t); |
| alpha = get_alpha(pw,pg,t,ipp); |
| fgu = -alpha*sg*rhoga*depsv_r; |
| fgu = fgu + n*rhoga*dsr_depsv*depsv_r;//volumetric strain rate effect |
| if(wrc_vol_strain_effect == 1) |
| fgu = fgu + n*rhoga*dsr_depsv*depsv_r;//volumetric strain rate effect |
| } |
| break; |
| } |
| 2924,8 → 3140,8 |
| alpha = get_alpha(pw,pg,t,ipp); |
| ftu = dhvap*sw*rhow*alpha*depsv_r; |
| ftu = ftu + dhvap*rhow*n*dsr_depsv*depsv_r;//volumetric strain rate effect |
| if(wrc_vol_strain_effect == 1) |
| ftu = ftu + dhvap*rhow*n*dsr_depsv*depsv_r;//volumetric strain rate effect |
| } |
| break; |
| } |
| 3162,6 → 3378,8 |
| print_err("gas pressure was modified in integration point No. %ld",__FILE__,__LINE__,__func__,ipp); |
| pg = 0.0; |
| } |
| check_math_err(); |
| } |
| 3177,8 → 3395,16 |
| TKr, 18.5.2005 |
| TKr 06/12/2022 - actualized |
| */ |
| void con_hawf3mat::waterpress_check(double &/*pw*/,double /*pg*/,double /*t*/,long /*ipp*/) |
| void con_hawf3mat::waterpress_check(double &pw,double pg,double t,long ipp) |
| { |
| if(pw < -5.0e8){ |
| //printf("water pressure was modified from value pw = %e to -5.0e8 in integration point No. %ld\n",pw,ipp); |
| //fprintf(Outt,"water pressure was modified from value pw = %e to -5.0e8 in integration point No. %ld\n",pw,ipp); |
| //print_err("water pressure was modified from value pw = %e to -5.0e8 in integration point No. %ld",__FILE__,__LINE__,__func__,pw,ipp); |
| pw = -5.0e8; |
| } |
| check_math_err(); |
| } |
| /** |
| 3394,7 → 3620,10 |
| 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(kintr_type < 3) |
| kintr = get_kintr(pw,pg,t,ipp); |
| else |
| kintr = get_kintrw(pw,pg,t,ipp); |
| if(rel_gas_press == 1) |
| pg = pg + p0; |
| 3441,8 → 3670,12 |
| 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); |
| //kintr = get_kintr(pw,pg,t,ipp); |
| if(kintr_type < 3) |
| kintr = get_kintr(pw,pg,t,ipp); |
| else |
| kintr = get_kintrg(pw,pg,t,ipp); |
| if(rel_gas_press == 1) |
| pg = pg + p0; |
| 3689,8 → 3922,12 |
| 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); |
| //kintr = get_kintr(pw,pg,t,ipp); |
| if(kintr_type < 3) |
| kintr = get_kintr(pw,pg,t,ipp); |
| else |
| kintr = get_kintrg(pw,pg,t,ipp); |
| if(rel_gas_press == 1) |
| pg = pg + p0; |
| 3979,8 → 4216,13 |
| rhow = get_rhow(t); |
| krw = get_krw(pw,pg,t,ipp); |
| muw = get_muw(t); |
| kintr = get_kintr(pw,pg,t,ipp); |
| //kintr = get_kintr(pw,pg,t,ipp); |
| if(kintr_type < 3) |
| kintr = get_kintr(pw,pg,t,ipp); |
| else |
| kintr = get_kintrw(pw,pg,t,ipp); |
| //correction 02/12/2022 |
| ktw = -dhvap*(rhow*kintr*krw/muw);//water |
| } |
| 4052,8 → 4294,12 |
| cpw = get_cpw(); |
| krw = get_krw(pw,pg,t,ipp); |
| muw = get_muw(t); |
| kintr = get_kintr(pw,pg,t,ipp); |
| //kintr = get_kintr(pw,pg,t,ipp); |
| if(kintr_type < 3) |
| kintr = get_kintr(pw,pg,t,ipp); |
| else |
| kintr = get_kintrw(pw,pg,t,ipp); |
| ktt2a = n*sw*rhow*cpw*kintr*krw/muw; |
| } |
| return(ktt2a); |
| 4090,8 → 4336,12 |
| 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); |
| //kintr = get_kintr(pw,pg,t,ipp); |
| if(kintr_type < 3) |
| kintr = get_kintr(pw,pg,t,ipp); |
| else |
| kintr = get_kintrg(pw,pg,t,ipp); |
| ktt2b = n*sg*rhog*cpg*kintr*krg/mug; |
| } |
| 4176,7 → 4426,7 |
| double con_hawf3mat::get_fw1(double pw,double pg,double t,long ipp) |
| { |
| double fc1; |
| double rhow,rhogw,rhog,krw,muw,krg,mug,kintr; |
| double rhow,rhogw,rhog,krw,muw,krg,mug,kintr,kintrg,kintrw; |
| if(model_type == artificial3) |
| fc1 = 0.0; |
| 4193,9 → 4443,16 |
| 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; |
| //kintr = get_kintr(pw,pg,t,ipp); |
| if(kintr_type < 3){ |
| kintr = get_kintr(pw,pg,t,ipp); |
| fc1 = rhogw*kintr*krg*rhog/mug + rhow*kintr*krw*rhow/muw; |
| } |
| else{ |
| kintrw = get_kintrw(pw,pg,t,ipp); |
| kintrg = get_kintrg(pw,pg,t,ipp); |
| fc1 = rhogw*kintrg*krg*rhog/mug + rhow*kintrw*krw*rhow/muw; |
| } |
| } |
| return(fc1); |
| } |
| 4227,8 → 4484,12 |
| 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); |
| //kintr = get_kintr(pw,pg,t,ipp); |
| if(kintr_type < 3) |
| kintr = get_kintr(pw,pg,t,ipp); |
| else |
| kintr = get_kintrg(pw,pg,t,ipp); |
| fg = rhoga*kintr*krg*rhog/mug; |
| fg = 0.0;//no gravity force is included for the moist air |
| } |
| 4263,7 → 4524,11 |
| rhow = get_rhow(t); |
| krw = get_krw(pw,pg,t,ipp); |
| muw = get_muw(t); |
| kintr = get_kintr(pw,pg,t,ipp); |
| //kintr = get_kintr(pw,pg,t,ipp); |
| if(kintr_type < 3) |
| kintr = get_kintr(pw,pg,t,ipp); |
| else |
| kintr = get_kintrw(pw,pg,t,ipp); |
| ft1 = -dhvap*rhow*kintr*krw*rhow/muw; |
| } |
| /trunk/SIFEL/TRFEL/SRC/consol_hawf3.h |
|---|
| 30,6 → 30,8 |
| double get_ds_dt(double pw, double pg, double t, long ipp); |
| double get_porosity(long ipp); |
| double get_kintr(double pw, double pg, double t, long ipp); |
| double get_kintrg(double pw, double pg, double t, long ipp); |
| double get_kintrw(double pw, double pg, double t, long ipp); |
| double get_krw(double pw, double pg, double t, long ipp); |
| double get_krg(double pw, double pg, double t, long ipp); |
| 179,7 → 181,8 |
| heatairwaterflowtype model_type; |
| int compress; |
| int vol_strain_effect,por_type,kintr_type,krw_type,krg_type,sr_type,xi_type,lambda_type,cps_type,betas_type,deff_type; |
| int vol_strain_effect,wrc_vol_strain_effect,por_type,kintr_type,krw_type,krg_type,sr_type,xi_type,lambda_type,cps_type,betas_type,deff_type; |
| int thermal_capacity_type; |
| double ma;//molar mass of dry air |
| double mw;//molar mass of water |
| double gasr;//universal gas constant |
| 209,6 → 212,7 |
| double b3,b4,b5; |
| double pr1,prif; |
| double muw0;//reference dynamic viscosity of water |
| double muw0const;//constant dynamic viscosity of water in 20 deg.C |
| double conb; |
| double conc; |
| double cpw0;//water specific heat at const. pressure = 4181 J/kg.K |
| 220,7 → 224,9 |
| double ks0; //inital bulk modulus of solid phase |
| double kt0; //inital bulk modulus of porous medium |
| double phi0; //inital porosity |
| double kintr0; //intial intrinsic permeability |
| double kintr0; //general intial intrinsic permeability (for both phases) |
| double kintrw0; //intial intrinsic permeability for liquid water, if it is separated from gas |
| double kintrg0; //intial intrinsic permeability for gas, if it is separated from water |
| double betas0; //inital cubic thermal dilatation coefficient |
| double betas_dry;//cubic thermal dilatation coefficient of dry soil |
| double betas_wet;//cubic thermal dilatation coefficient of saturated soil |
| 237,6 → 243,8 |
| double lambda_wet;//effective thermal conductivity of saturated soil |
| double sr_dry; //saturation degree according to lambda_dry |
| double sr_wet; //saturation degree according to lambda_wet |
| double rhocp_dry; //thermal capacity according to lambda_dry |
| double rhocp_wet; //thermal capacity according to lambda_wet |
| double tau0; //tortuosity factor 0.4 to 0.6 |
| double mefel_units; //basic units for pressures = Pa (Pascals) |
| 244,7 → 252,7 |
| double scale_pw,scale_pg,scale_t; |
| int rel_gas_press; //relative gas pressure according to ambient air; yes=1=relative pg, no=0=absolute pg |
| double krw0,sirr,ssat,lambda_krw,beta_krw; |
| double bb1,phi01; //intrinsic permeability parameters |
| double bb1,phi01,ng; //intrinsic permeability parameters |
| double krg0,s_crit,ag; //gas relative permeability parameters |
| //parameters for the arificial material: |
| /trunk/SIFEL/TRFEL/SRC/consol_wf1.cpp |
|---|
| 5,7 → 5,10 |
| sources: Lewis and Schrefler's book |
| unknowns: number of unknowns=1, pw = liqud water pressure |
| */ |
| #include <iostream> |
| #include <cmath> |
| #include <cstdlib> |
| #include <string.h> |
| #include <stdio.h> |
| #include <stdlib.h> |
| #include <math.h> |
| 14,6 → 17,7 |
| #include "consol_wf1.h" |
| #include "globalt.h" |
| #include "globmatt.h" |
| #include "errno.h" |
| con_wf1mat::con_wf1mat() |
| { |
| 26,11 → 30,12 |
| mefel_units = 1.0; //basic units for pressures = Pa (Pascals) in MEFEL part |
| vol_strain_effect = 0; //volumetric strain rate influence: 0=no; 1=yes |
| wrc_vol_strain_effect = 0; //volumetric strain rate influence on water retention curve: 0=no; 1=yes |
| alpha = 1.0; //Biot's constant [-] alpha = 1 - kt/ks |
| ks = 2.167e9; //bulk modulus of solid phase (grains) [Pa] |
| ks0 = 2.167e9; //bulk modulus of solid phase (grains) [Pa] |
| kw = 2.0e9; //bulk modulus of water [Pa] |
| phi0 = 0.297; //initial porosity [-] |
| phi0 = 0.4927; //initial porosity [-] |
| kintr0 = 4.5e-13; //intrinsic permeability [m^2] |
| rhow0 = 1000.0; //initial water density [kg/m^3] |
| muw0 = 1.0e-3; //initial water viscosity [Pa.s] |
| 50,6 → 55,17 |
| gamma = 0.0; |
| lambda0 = 0.0; |
| s_entry = 0.0; |
| //parameters for Masin hypopl. eff. stress factor chi - debug version: |
| sairentry0 = -2700.0; |
| aer = 1.0; |
| a_scan = 0.0; |
| eM0 = 0.2; |
| kappam = 0.1; |
| smstar = -1000.0; |
| emstar = 1.0; |
| lambdap0 = 1.0; |
| em = 0.5; //micro porosity //basic setup |
| } |
| con_wf1mat::~con_wf1mat() |
| 190,7 → 206,7 |
| xfscanf (in,"%k%m","waterflowtype",&waterflowtype_kwdset, &model_type); //water flow model type |
| xfscanf (in,"%d", &compress); //compressibility of grains |
| // common material parameters |
| xfscanf (in,"%le %le %le %le %le %le %d %d %d %d %d", &alpha, &ks, &rhos0, &pw_bc, &t0, &kintr0, &por_type, &kintr_type, &krw_type, &sr_type, &xi_type); |
| xfscanf (in,"%le %le %le %le %le %le %d %d %d %d %d", &alpha, &ks0, &rhos0, &pw_bc, &t0, &kintr0, &por_type, &kintr_type, &krw_type, &sr_type, &xi_type); |
| switch (model_type){ |
| case lewis_and_schrefler:{//Lewis and Schrefler's model approach |
| 197,7 → 213,7 |
| break; |
| } |
| case lewis_and_schrefler_mefel:{//Lewis and Schrefler's model approach coupled with MEFEL |
| xfscanf (in,"%le %d", &mefel_units, &vol_strain_effect); |
| xfscanf (in,"%le %d %d", &mefel_units, &vol_strain_effect, &wrc_vol_strain_effect); |
| break; |
| } |
| default:{ |
| 344,6 → 360,10 |
| } |
| case masin_xi: |
| break; |
| case masin_hypopl_xi:{ |
| xfscanf (in,"%le %le %le %le %le %le %le %le %le", &sairentry0,&aer,&a_scan,&eM0,&kappam,&smstar,&emstar,&lambdap0,&em); |
| break; |
| } |
| default:{ |
| print_err ("unknown model type is required", __FILE__, __LINE__, __func__); |
| abort(); |
| 504,22 → 524,43 |
| //xi = gamma*get_sw(pw,ipp); |
| sr = get_sw(pw,ipp); |
| //xi = pow(sr,(gamma/lambda0)); |
| xi = pow(sr,gamma); |
| xi = pow(sr,lambda0); |
| xi = (1-gamma)*xi; |
| //debug??!! |
| if (sr < 1.0) |
| xi = 0; |
| else |
| xi = 1.0; |
| break; |
| case biot_masin_xi:{//according to masin for testing |
| sr = get_sw(pw,ipp); |
| suc = -pw; |
| if (suc>=s_entry) |
| if (suc>=s_entry){ |
| //xi = get_sw(pw,ipp); |
| xi = pow((s_entry/suc),gamma); |
| else |
| //if (xi < 0.1) |
| //xi = 0.0; |
| //else |
| xi = (1-gamma*(1-xi))*xi; |
| //xi = pow(xi,2.0); |
| } |
| else{ |
| xi = 1.0; |
| if (suc>=s_entry) |
| xi = (1-gamma)*xi; |
| else |
| xi = 1.0; |
| } |
| break; |
| } |
| /* case biot_masin_xi:{//according to masin for testing |
| suc = -pw; |
| if (suc>=s_entry){ |
| xi = pow((s_entry/suc),gamma); |
| xi = (1-gamma)*xi; |
| } |
| else{ |
| xi = 1.0; |
| } |
| break; |
| } |
| */ |
| case masin_xi:{ |
| double e=0.0,dpw=0.0,por=0.0; |
| if(Tm->nontransq != NULL){ |
| 532,13 → 573,40 |
| dpw = Tm->ip[ipp].av[0]-Tm->ip[ipp].pv[0]; |
| xi = masin_ret.psi(-pw,-dpw,e,t0);//positive value of suction |
| break; |
| } |
| case masin_hypopl_xi:{ |
| double e=0.0,dpw=0.0,por=0.0; |
| if(Tm->nontransq != NULL){ |
| por = Tm->givenontransq(porosity, ipp);// from mefel |
| e = por/(1-por); |
| } |
| else{ |
| e = phi0/(1-phi0); |
| } |
| if (e == 0.0) |
| e = phi0/(1-phi0); |
| //debug ??!! |
| e = phi0/(1-phi0); |
| dpw = Tm->ip[ipp].av[0]-Tm->ip[ipp].pv[0]; |
| xi = masin_hypopl_psi(pw,dpw,e,ipp); |
| break; |
| } |
| default:{ |
| print_err ("unknown model type is required", __FILE__, __LINE__, __func__); |
| abort(); |
| } |
| } |
| if (ipp == 0){ |
| fprintf(Outt," ipp = %ld, suction = %e,xi = %e\n",ipp,pw,xi); |
| fflush(Outt); |
| } |
| return(xi); |
| } |
| 767,6 → 835,7 |
| } |
| case 7:{ |
| krw = van_genuchten_ret.get_krw(-pw,t0); |
| break; |
| } |
| case 9:{//bazant |
| sw = get_sw(pw,ipp); |
| 863,11 → 932,50 |
| */ |
| double con_wf1mat::get_alpha() |
| { |
| return(alpha); |
| return(alpha); //is set as constant value |
| } |
| /** |
| function computes bulk modulus of solid phase |
| @param ipp - number of integration point |
| @retval ks - bulk modulus of solid phase |
| 16/11/2023 TKr |
| */ |
| double con_wf1mat::get_ks(long ipp) |
| { |
| double ks; |
| ks = ks0; |
| switch (model_type){ |
| case artificial3: |
| case lewis_and_schrefler3:{//Lewis and Schrefler's book |
| break; |
| } |
| case lewis_and_schrefler3_mefel:{//Lewis and Schrefler's book, bulk modulus obtained from mefel; |
| if(Tm->nontransq != NULL){ |
| ks = Tm->givenontransq(bulk_modulus, ipp);// from mefel |
| // if it is calculated from the actual stiffness matrix, it is bulk modulus of solid skeleton, not the solid matrix (grains) |
| ks = ks*mefel_units; |
| } |
| break; |
| } |
| default:{ |
| print_err("unknown model type is required", __FILE__, __LINE__, __func__); |
| } |
| } |
| return(ks); |
| } |
| /** |
| function returns kw bulk modulus of water [Pa] |
| @retval kw - bulk modulus of water [Pa] |
| 922,12 → 1030,14 |
| */ |
| double con_wf1mat::get_capww(double pw, long ipp) |
| { |
| double sw,dsw_dpw,n,capww; |
| double sw,dsw_dpw,n,ks,capww; |
| sw = get_sw(pw,ipp); |
| dsw_dpw = get_dsw_dpw(pw,ipp); |
| n = get_porosity(ipp); |
| ks = get_ks(ipp); |
| if(compress == 1){ |
| //compressible grains: |
| capww = (alpha-n)/ks*sw*(sw + dsw_dpw*pw)+n*sw/kw; |
| 1046,11 → 1156,12 |
| } |
| Tm->ip[ipp].eqother[1] = dsr_depsv; //this is not necessary |
| sw = get_sw(pw,ipp); |
| fwu = -alpha*sw*depsv_r;//volumetric strain effect |
| //if (sr_type == mefel_sr){ |
| fwu = fwu - n*dsr_depsv*depsv_r;//volumetric strain effect on retention curve |
| if(wrc_vol_strain_effect == 1) |
| fwu = fwu - n*dsr_depsv*depsv_r;//volumetric strain effect on retention curve |
| //} |
| } |
| break; |
| 1353,7 → 1464,7 |
| switch (compother){ |
| case 0:{//capillary pressure |
| other = -pw; |
| other = -pw/1000.0; //in kPa //debug??!! |
| break; |
| } |
| case 1:{//saturation |
| 1367,6 → 1478,10 |
| case 3:{//specific moisture content |
| other = get_w(pw,ipp); |
| break; |
| } |
| case 4:{//stress factor chi |
| other = get_xi(pw,ipp); |
| break; |
| } |
| default:{ |
| print_err("unknown type of component is required in function", __FILE__, __LINE__, __func__); |
| 1403,7 → 1518,11 |
| fprintf (out,"Moisture (water) content (kg/kg) "); |
| //fprintf (out,"Moisture content (kg/kg) "); |
| break; |
| } |
| } |
| case 4:{//stress factor chi |
| fprintf (out,"Effective stress chi (-) "); |
| break; |
| } |
| default:{ |
| print_err ("unknown type of component is required in function", __FILE__, __LINE__, __func__); |
| abort(); |
| 1670,11 → 1789,12 |
| } |
| case lewis_and_schrefler_mefel:{//Lewis and Schrefler's book, saturation degree ind its derivative are obtained from mefel; |
| antq[porosity-1] = 1; |
| antq[bulk_modulus-1] = 1; |
| antq[mmean_stress-1] = 1; |
| if (sr_type == mefel_sr){ |
| antq[saturation_deg-1] = 1; |
| antq[der_saturation_deg-1] = 1; |
| antq[der_saturation_deg_depsv-1] = 1; |
| antq[bulk_modulus-1] = 1; |
| } |
| if(vol_strain_effect == 1) |
| antq[strain_vol_rate-1] = 1; |
| 1686,3 → 1806,84 |
| } |
| } |
| } |
| /** |
| function computes effective stress factor xi accordin to generlmod from D. Masin; |
| this was created for the use of effective stress factor chi for bentonite materials |
| - not working properly due to incorrect initialization of effective stresses (initial conditions) in mechanical models |
| @param s - suction |
| @param ds - suction increment |
| @retval xi - factor xi |
| 03/10/2023, TKr |
| */ |
| double con_wf1mat::masin_hypopl_psi(double suction,double dsuction, double e, long ipp) |
| { |
| double xi=0.0,sr=0.0,SrM=0.0,rlambda=0.0,sewM=0.0; |
| //double rlambda_for_se=0.0; |
| double scanpower=3; |
| double SrMlimit=0.2; //was SrMlimit=0.75; |
| double SrMmax=1.0; |
| double WRCpower=1.0; //was WRCpower=1.1; |
| //double Sepower=3.0; |
| double eM = 0.0; //macro porosity |
| double evoid = 0.0; |
| double signetm; |
| gamma = lambdap0; |
| evoid = e; //total porosity |
| signetm = Tm->givenontransq(mmean_stress, ipp);// from mefel - net stress is equal to total stress |
| double pefsat=signetm+suction/1000.0; |
| double em_pefsat=exp(kappam*log(smstar/pefsat)+log(1+emstar))-1; |
| em=exp(log(1+em_pefsat))-1.0; |
| double min_void=0.0001; |
| if(em<min_void) em=min_void; |
| if(em>evoid-min_void) em=evoid-min_void; |
| eM=(evoid-em)/(1.0+em); //macro porosity |
| sewM=aer*sairentry0*eM0/eM; |
| sr = get_sw(suction,ipp); //actual degree saturation |
| //Sr = SrM + em/e*(Srm - SrM); //Srm = saturation degree at microlevel |
| //Sr = (SrM*(e-em)+em)/e; |
| SrM = (sr*e-em)/(e-em); |
| double fact=a_scan; //wetting within main curves |
| if(dsuction>0) fact=1-a_scan; //drying within main curves |
| rlambda=1; |
| //rlambda_for_se=1; |
| if(fact<1.e-10) fact=0; |
| if(aer<1) |
| rlambda=pow(fact, scanpower); |
| if(SrM>1) |
| SrM=1; |
| if(suction>=sewM) {//drying from 0 |
| rlambda=0; |
| //rlambda_for_se=0; |
| } |
| //wetting to 0, rlambda_for_ascan remains 1 here |
| if(dsuction>0 && SrM>SrMlimit) { |
| rlambda=pow((SrMmax-SrM)/(SrMmax-SrMlimit), WRCpower); |
| //rlambda_for_se=pow((SrMmax-SrM)/(SrMmax-SrMlimit), Sepower); |
| } |
| xi = SrM*(1-gamma*rlambda); |
| //xi = SrM; |
| if (ipp == 0){ |
| fprintf(Outt," ipp = %ld, suction = %e, SrM = %e, sr = %e, rlambda = %e, xi = %e, e = %e, signetm = %e\n",ipp,suction,SrM,sr,rlambda,xi,e,signetm); |
| fflush(Outt); |
| } |
| check_math_errel(Tm->elip[ipp]); |
| return xi; |
| } |
| /trunk/SIFEL/TRFEL/SRC/consol_wf1.h |
|---|
| 12,6 → 12,12 |
| #include "masin_reten.h" |
| #include "febex_granit_reten.h" |
| //debug version: |
| #include <iostream> |
| #include <cmath> |
| #include <cstdlib> |
| #include "errno.h" |
| /** |
| This class defines model for water flow in soils. |
| */ |
| 46,6 → 52,7 |
| double get_kintr(long ipp); |
| double get_k(double pw,long ipp); |
| double get_alpha(); |
| double get_ks(long ipp); |
| double get_kw(); |
| double get_muw(); |
| 83,6 → 90,7 |
| double give_suction(long ipp); |
| double give_saturation_degree(long ipp); |
| double get_w(double pw,long ipp); |
| double masin_hypopl_psi(double suction,double dsuction, double e, long ipp); |
| /// marks required non-transport quantities |
| void give_reqntq(long *antq); |
| 110,9 → 118,9 |
| waterflowtype model_type; |
| int compress; |
| int vol_strain_effect,por_type,kintr_type,krw_type,sr_type,xi_type; |
| int vol_strain_effect,wrc_vol_strain_effect,por_type,kintr_type,krw_type,sr_type,xi_type; |
| double mefel_units; //basic units for pressures = Pa (Pascals) |
| double alpha,phi0,ks; |
| double alpha,phi0,ks0; |
| double kw,rhow0,muw0,t0; |
| double kintr0,rhos0; |
| double krw0,sirr,ssat,lambda_krw,beta_krw; |
| 119,6 → 127,9 |
| double b1,phi01; |
| double pw_bc; |
| double gamma,lambda0,s_entry; |
| //debug version: |
| double sairentry0,aer,a_scan,eM0,kappam,smstar,emstar,lambdap0,em; |
| }; |
| #endif |
| /trunk/SIFEL/TRFEL/SRC/globmatt.cpp |
|---|
| 1185,7 → 1185,7 |
| // implementation due to saltmat4 |
| copy_nodval (0); |
| } |
| // loop over elements |
| for (i=0;i<Tt->ne;i++){ |
| if (Gtt->leso[i]==1){ |
| 1199,7 → 1199,6 |
| } |
| Tm->aux_values (0); |
| if (Tp->nvs==1){ |
| // arrays of actual values are moved to arrays of previous values |
| 1296,6 → 1295,8 |
| Tm->aux_values (0); |
| // arrays of actual values are moved to arrays of previous values |
| actual_previous_change (); |
| if (Tp->nvs==1){ |
| // arrays of actual values are moved to arrays of previous values |
| /trunk/SIFEL/TRFEL/SRC/npsolvert.cpp |
|---|
| 619,6 → 619,8 |
| // (C + alpha dt K) v_{n+1} = f_{n+1} - K dd |
| // time derivatives v_{n+1} are solved |
| Tp->ssle->solve_system (Gtt,Cmat,tdlhs,f,Outt); |
| check_math_err();//debug??!! |
| // nodal values computed from nodal derivatives |
| // d_{n+1} = dd + alpha dt v_{n+1} |
| 626,6 → 628,9 |
| lhs[j] = d[j] + alpha*dt*tdlhs[j]; |
| } |
| // arrays of actual values are moved to arrays of previous values |
| actual_previous_change (); |
| // physically corrected solution |
| solution_correction (); |
| /trunk/SIFEL/TRFEL/SRC/transmat.cpp |
|---|
| 2485,7 → 2485,7 |
| switch (Tp->mednam)//names of transported media |
| { |
| case water:{ |
| ncompo=4; |
| ncompo=5; |
| break; |
| } |
| case moisture:{ |
| /trunk/SIFEL/TRFEL/SRC/van_genuchten.cpp |
|---|
| 78,7 → 78,7 |
| xfscanf (in,"%lf %lf %lf %lf %lf %lf", &ssat, &sirr, &p0, &lambda0, &pd, &lambdap); |
| break; |
| } |
| case 3:{//extended for hight suction - other expression |
| case 3:{//extended for hight suction and temperature effect |
| xfscanf (in,"%lf %lf %lf %lf %lf %lf %lf %lf", &ssat, &sirr, &p0, &lambda0, &pd, &lambdap, &sig0, &t0); |
| break; |
| } |
| 108,11 → 108,11 |
| fprintf (out,"\n%lf %lf %lf %lf\n", ssat, sirr, p0, lambda0); |
| break; |
| } |
| case 2:{//extended for expansive clays |
| case 2:{//extended for hight suction |
| fprintf (out,"\n%lf %lf %lf %lf %lf %lf\n", ssat, sirr, p0, lambda0, pd, lambdap); |
| break; |
| } |
| case 3:{//extended for expansive clays and temperature dependent |
| case 3:{//extended for hight suction and temperature effect |
| fprintf (out,"\n%lf %lf %lf %lf %lf %lf %lf %lf\n", ssat, sirr, p0, lambda0, pd, lambdap, sig0, t0); |
| break; |
| } |
| 185,9 → 185,9 |
| a = (374.15-(t-273.15))/647.3;// temperature in deg.C |
| sigt = (1.0 - 0.625*a)*(0.2358*pow(a,1.256)); |
| } |
| else |
| else{ |
| sigt = 0.0019106*exp(0.05*(633.15-t)); |
| } |
| p0 = p0*sigt/sig0; |
| expm = 1.0/(1.0-lambda0); |
| fd = pow((1.0-(-pw)/pd),lambdap); |
| 203,6 → 203,9 |
| } |
| } |
| if(sw > 1.0) |
| sw = 1.0; |
| return(sw); |
| } |
| 218,13 → 221,15 |
| */ |
| double van_genuchten_reten::dsw_dpw(double pw, double t) |
| { |
| double dsw_dpw,hp,fd,theta,dtheta_dpw,dfd_dpw,sigt,a; |
| double sr,dsw_dpw,hp,fd,theta,dtheta_dpw,dfd_dpw,sigt,a; |
| sr = sw(pw,t); |
| switch (vg_ret_type){ |
| case 0:{//classical expression |
| pw = pw/1000.0;//in kPa |
| if(pw >= 0.0){ |
| if(pw >= 0.0 || sr == 1.0){ |
| dsw_dpw = 0.0; |
| } |
| else{ |
| 238,7 → 243,7 |
| break; |
| } |
| case 1:{//other expression (in Pa) |
| if(pw >= 0.0){ |
| if(pw >= 0.0 || sr == 1.0){ |
| dsw_dpw = 0.0; |
| } |
| else{ |
| 250,7 → 255,7 |
| break; |
| } |
| case 2:{//extended for hight suction (in Pa) |
| if(pw >= 0.0){ |
| if(pw >= 0.0 || sr == 1.0){ |
| dsw_dpw = 0.0; |
| } |
| else{ |
| 265,7 → 270,7 |
| break; |
| } |
| case 3:{//extended for hight suction and temperature effect (in Pa) |
| if(pw >= 0.0){ |
| if(pw >= 0.0 || sr == 1.0){ |
| dsw_dpw = 0.0; |
| } |
| else{ |
| 276,9 → 281,9 |
| a = (374.15-(t-273.15))/647.3;// temperature in deg.C |
| sigt = (1.0 - 0.625*a)*(0.2358*pow(a,1.256)); |
| } |
| else |
| else{ |
| sigt = 0.0019106*exp(0.05*(633.15-t)); |
| } |
| p0 = p0*sigt/sig0; |
| theta = pow(1+(pow((-pw/p0),expm)),-lambda0); |
| dtheta_dpw = -lambda0*pow(1+(pow((-pw/p0),expm)),-lambda0-1.0)*expm*pow((-pw/p0),expm-1)*(-1.0/p0); |
| 309,10 → 314,12 |
| */ |
| double van_genuchten_reten::dsw_dt(double pw, double t) |
| { |
| double dsw_dt; |
| double sr,dsw_dt,fd,a,da_dt,dsigt_dt,dp0_dt,dtheta_dt; |
| sr = sw(pw,t); |
| switch (vg_ret_type){ |
| case 0://classical expression |
| case 0://classical expression (in kPa) |
| case 1://other expression (in Pa) |
| case 2:{//extended for hight suction (in Pa) |
| dsw_dt = 0.0; |
| 320,6 → 327,30 |
| } |
| case 3:{//extended for hight suction and temperature effect (in Pa) - not implemented yet |
| dsw_dt = 0.0; |
| if(pw >= 0.0 || sr == 1.0){ |
| dsw_dt = 0.0; |
| } |
| else{ |
| expm = 1.0/(1.0-lambda0); |
| fd = pow((1.0-pw/pd),lambdap); |
| if (t < 633.15){ |
| a = (374.15-(t-273.15))/647.3;// temperature in deg.C |
| da_dt = -1.0/647.3; |
| //sigt = (1.0 - 0.625*a)*(0.2358*pow(a,1.256)); |
| dsigt_dt = -0.332478*pow(a,0.256)*(a-0.89078)*da_dt; |
| } |
| else{ |
| //sigt = 0.0019106*exp(0.05*(633.15-t)); |
| dsigt_dt = -5.335571e9*pow(0.951229,t); |
| } |
| dp0_dt = p0*dsigt_dt/sig0; |
| dtheta_dt = -lambda0*pow(1+(pow((-pw/p0),expm)),-lambda0-1.0)*expm*pow((-pw/p0),expm-1)*(pw/p0/p0)*dp0_dt; |
| dsw_dt = dtheta_dt*fd; |
| dsw_dt = (ssat - sirr)*dsw_dt; |
| } |
| break; |
| } |
| default:{ |