Commit 1e33437f authored by Martin Beck's avatar Martin Beck
Browse files

[feature/ShearStressSink] Most recent version used for Promo

parent ae758e28
......@@ -259,19 +259,19 @@ public:
// calculate the effective porosity
if(problem.coupled() == true)
{
if ((*this)[scvIdx].divU < -(*this)[scvIdx].porosity())
{
if ((*this)[scvIdx].divU < -(*this)[scvIdx].porosity())
{
(*this)[scvIdx].effPorosity = (*this)[scvIdx].porosity();
std::cout<<"volume change too large"<<std::endl;
}
else
// this equation would be correct if the bulk volume could change (Vol_new = Vol_init *(1+div u)), however, we
// have a constant bulk volume therefore we should apply phi_eff = phi_init + div u
// but this causes convergence problems. Since div u is very small here the chosen relation is
// assumed to be a good approximation
(*this)[scvIdx].effPorosity = ((*this)[scvIdx].porosity() + (*this)[scvIdx].divU)/(1.0 + (*this)[scvIdx].divU);
(*this)[scvIdx].effPorosity = (*this)[scvIdx].porosity();
std::cout<<"volume change too large"<<std::endl;
}
else
// this equation would be correct if the bulk volume could change (Vol_new = Vol_init *(1+div u)), however, we
// have a constant bulk volume therefore we should apply phi_eff = phi_init + div u
// but this causes convergence problems. Since div u is very small here the chosen relation is
// assumed to be a good approximation
(*this)[scvIdx].effPorosity = ((*this)[scvIdx].porosity() + (*this)[scvIdx].divU)/(1.0 + (*this)[scvIdx].divU);
}
else
(*this)[scvIdx].effPorosity = (*this)[scvIdx].porosity();
}
......
......@@ -167,16 +167,16 @@ public:
// this evaluation should be moved to another location
DimVector tmpVec;
DimMatrix Keff, Keff_i, Keff_j;
// Keff_i = EffectivePermeabilityModel::effectivePermeability(this->curVolVars_()[fluxVars.face().i],
// this->problem_().spatialParams(),
// this->element_(),
// this->fvGeometry_(),
// fluxVars.face().i);
// Keff_j = EffectivePermeabilityModel::effectivePermeability(this->curVolVars_()[fluxVars.face().j],
// this->problem_().spatialParams(),
// this->element_(),
// this->fvGeometry_(),
// fluxVars.face().j);
Keff_i = EffectivePermeabilityModel::effectivePermeability(this->curVolVars_()[fluxVars.face().i],
this->problem_().spatialParams(),
this->element_(),
this->fvGeometry_(),
fluxVars.face().i);
Keff_j = EffectivePermeabilityModel::effectivePermeability(this->curVolVars_()[fluxVars.face().j],
this->problem_().spatialParams(),
this->element_(),
this->fvGeometry_(),
fluxVars.face().j);
this->problem_().spatialParams().meanK(Keff, Keff_i, Keff_j);
// loop over all phases
......@@ -246,7 +246,7 @@ public:
stabilizationTerm /= dt;
// flux[eqIdx] -= stabilizationTerm;
flux[eqIdx] -= stabilizationTerm;
}
}
......
......@@ -490,7 +490,7 @@ public:
const Scalar lambda = lameParams[0];
const Scalar mu = lameParams[1];
K[eIdx] = lambda + 2/3*mu;
K[eIdx] = lambda + 2.0/3.0*mu;
stressSink[eIdx] = stressSink_[eIdx];
// calculate strain tensor
......@@ -512,8 +512,10 @@ public:
}
if(hasElementFailed_[eIdx])
sigma = this->calculateReducedStress(eIdx, sigma, true);
if (eIdx == 1)
sigma = this->calculateReducedStress(eIdx, sigma, true);
else
sigma = this->calculateReducedStress(eIdx, sigma, false);
// in case of rock mechanics sign convention compressive stresses
// are defined to be positive
if(rockMechanicsSignConvention_){
......@@ -610,10 +612,10 @@ public:
Dune::FieldMatrix<RF, dim, dim> deltaEffStressRotated(0.0);
deltaEffStressRotated = this->calculateRotatedStress(deltaEffStress, faultAngle_[eIdx]);
// if( (eIdx == 1019) /*|| (eIdx == 1025)*/)
// if( (eIdx == 1) /*|| (eIdx == 1025)*/)
// // if(this->problem().getHasElementFailed(eIdx))
// {
// std::cout << "element " << eIdx << " is inclined by " << angle << std::endl;
// std::cout << "element " << eIdx << " is inclined by " << faultAngle_[eIdx] << std::endl;
//
// std::cout << "deltaEffStress:" << std::endl;
// std::cout << deltaEffStress[0][0] << " " << deltaEffStress[0][1] << std::endl;
......
......@@ -821,33 +821,34 @@ public:
// if( this->timeManager().episodeIndex() == 7)
// episodeLength_ = 1.0;
if( this->timeManager().episodeIndex() == 1)
episodeLength_ = 1800.0;
if( this->timeManager().episodeIndex() == 2)
episodeLength_ = 100.0;
// if( this->timeManager().episodeIndex() == 1)
// episodeLength_ = 1800.0;
// if( this->timeManager().episodeIndex() == 2)
// episodeLength_ = 100.0;
if( this->timeManager().episodeIndex() == 15)
episodeLength_ = 0.98;
if( (this->timeManager().episodeIndex() == 24) ||
(this->timeManager().episodeIndex() == 34))
episodeLength_ = 0.99;
if( (this->timeManager().episodeIndex() == 16) ||
(this->timeManager().episodeIndex() == 25) ||
(this->timeManager().episodeIndex() == 35))
episodeLength_ = 4.00;
if( (this->timeManager().episodeIndex() == 17) ||
(this->timeManager().episodeIndex() == 26) ||
(this->timeManager().episodeIndex() == 36))
episodeLength_ = 15.00;
if( (this->timeManager().episodeIndex() == 18) ||
(this->timeManager().episodeIndex() == 27) ||
(this->timeManager().episodeIndex() == 37))
episodeLength_ = 20.00;
if( (this->timeManager().episodeIndex() == 22) ||
(this->timeManager().episodeIndex() == 31) ||
(this->timeManager().episodeIndex() == 41))
episodeLength_ = 100.00;
// if( this->timeManager().episodeIndex() == 15)
// episodeLength_ = 0.98;
// if( (this->timeManager().episodeIndex() == 24) ||
// (this->timeManager().episodeIndex() == 34))
// episodeLength_ = 0.99;
// if( (this->timeManager().episodeIndex() == 16) ||
// (this->timeManager().episodeIndex() == 25) ||
// (this->timeManager().episodeIndex() == 35))
// episodeLength_ = 4.00;
// if( (this->timeManager().episodeIndex() == 17) ||
// (this->timeManager().episodeIndex() == 26) ||
// (this->timeManager().episodeIndex() == 36))
// episodeLength_ = 15.00;
// if( (this->timeManager().episodeIndex() == 18) ||
// (this->timeManager().episodeIndex() == 27) ||
// (this->timeManager().episodeIndex() == 37))
// episodeLength_ = 20.00;
// if( (this->timeManager().episodeIndex() == 22) ||
// (this->timeManager().episodeIndex() == 31) ||
// (this->timeManager().episodeIndex() == 41))
// episodeLength_ = 100.00;
// High K
// if( this->timeManager().episodeIndex() == 29)
// episodeLength_ = 0.96;
// if( (this->timeManager().episodeIndex() == 30) ||
......@@ -866,15 +867,33 @@ public:
// (this->timeManager().episodeIndex() == 36) ||
// (this->timeManager().episodeIndex() == 36))
// episodeLength_ = 100.00;
//
// if( this->timeManager().episodeIndex() == 1)
// episodeLength_ = 1800.0;
// if( this->timeManager().episodeIndex() == 2)
// episodeLength_ = 100.0;
// if( (this->timeManager().episodeIndex() == 9) ||
// (this->timeManager().episodeIndex() == 23))
// episodeLength_ = 0.98;
// Low K
if( this->timeManager().episodeIndex() == 1)
episodeLength_ = 900.0;
if( this->timeManager().episodeIndex() == 2)
episodeLength_ = 100.0;
if( this->timeManager().episodeIndex() == 10)
episodeLength_ = 0.99;
if( (this->timeManager().episodeIndex() == 11) ||
(this->timeManager().episodeIndex() == 11) ||
(this->timeManager().episodeIndex() == 11))
episodeLength_ = 4.00;
if( (this->timeManager().episodeIndex() == 12) ||
(this->timeManager().episodeIndex() == 12) ||
(this->timeManager().episodeIndex() == 12))
episodeLength_ = 15.00;
if( (this->timeManager().episodeIndex() == 13) ||
(this->timeManager().episodeIndex() == 13) ||
(this->timeManager().episodeIndex() == 13))
episodeLength_ = 20.00;
if( (this->timeManager().episodeIndex() == 17) ||
(this->timeManager().episodeIndex() == 17) ||
(this->timeManager().episodeIndex() == 17))
episodeLength_ = 100.00;
// if( this->timeManager().episodeIndex() == 12)
// episodeLength_ = 0.99;
// if( (this->timeManager().episodeIndex() == 13) ||
......
......@@ -640,7 +640,7 @@ public:
{
if(globalPos[1] > this->bBoxMax()[1]-eps_)
{
Scalar time = this->timeManager().time();
Scalar time = this->timeManager().time() - tInitEnd_;
Scalar factor= GET_RUNTIME_PARAM(TypeTag, Scalar,FailureParameters.DisplacementFactor);
values[uyIdx] = factor*time;
......
......@@ -374,8 +374,8 @@ private:
Scalar xCoord = globalPos[0];
Scalar yCoord = globalPos[1];
Scalar xCenterFault = 50.0;
Scalar yCenterFault = 50.0;
Scalar xCenterFault = 0.5;
Scalar yCenterFault = 0.5;
// Scalar xFaultDistanceCenterToBothTips = 87.5; // x-distance in m to both tips of the fault
Scalar xFaultDistanceCenterToBothTips = 100; // x-distance in m to both tips of the fault
......@@ -385,18 +385,19 @@ private:
Scalar deltaX = xCenterFault - xCoord;
if (std::abs(deltaX) < xFaultDistanceCenterToBothTips + eps_ &&
if(
(std::abs(deltaX) < xFaultDistanceCenterToBothTips + eps_ &&
yCoord < yCenterFault - tan(55 * M_PI / 180) * deltaX + spatialTolarance_ &&
yCoord > yCenterFault - tan(55 * M_PI / 180) * deltaX - spatialTolarance_ &&
// yCoord < 50 &&
// yCoord > 40 &&
yCoord > yFaultBottom &&
yCoord < yFaultTop)
return true;
else if (xCoord < 480 &&
xCoord > 470 &&
yCoord > 1020 &&
yCoord < 1030)
// for Heterogeneous
&&
(!(xCoord < 0.21 && xCoord > 0.09 && yCoord > 0.00 && yCoord < 0.11))
&&
(!(xCoord < 0.91 && xCoord > 0.79 && yCoord > 0.89 && yCoord < 1.01)) )
return true;
// else if (xCoord < 50 &&
// xCoord > 45 &&
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment