/usr/include/ASL/num/aslFDPoroElasticity.h is in libasl-dev 0.1.6-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 | /*
* Advanced Simulation Library <http://asl.org.il>
*
* Copyright 2015 Avtech Scientific <http://avtechscientific.com>
*
*
* This file is part of Advanced Simulation Library (ASL).
*
* ASL is free software: you can redistribute it and/or modify it
* under the terms of the GNU Affero General Public License as
* published by the Free Software Foundation, version 3 of the License.
*
* ASL is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Affero General Public License for more details.
*
* You should have received a copy of the GNU Affero General Public License
* along with ASL. If not, see <http://www.gnu.org/licenses/>.
*
*/
#ifndef ASLFDPOROELASTICITY_H
#define ASLFDPOROELASTICITY_H
#include "aslNumMethod.h"
#include "acl/aclMath/aclVectorOfElementsDef.h"
#include "aslFDElasticity.h"
#include "utilities/aslUValue.h"
namespace asl
{
/// Numerical method which computes homogenious isotropic poro-elasticity equation
/**
\ingroup Elasticity
\ingroup NumMethods
The classic poroelastic equations were originally developed by Biot (1941)
to represent biphasic soil consolidation. The equations can be written
in the following form:
\f[ (K+\mu/3)\nabla_j \nabla_k u_k+ \mu \Delta u_j - a \nabla_j p = - (\rho_s-\rho_f)g_j \f]
\f[ a\nabla_i\dot u_i+ \frac{1}{S} \dot p = \nabla_i k \nabla_i p, \f]
where
- \f$\vec u\f$ is the displacement vector,
- \f$ p \f$ the interstitial pressure,
- \f$ a \f$ the ratio of fluid volume extracted to volume change of the tissue under compression,
- \f$\rho_s\f$ the solid fraction density,
- \f$\rho_f\f$ the fluid density,
- \f$\vec g \f$ the gravitational vector,
- \f$ 1/S \f$ the amount of fluid which can be forced into the tissue under constant volume,
- \f$ k \f$ the hydraulic conductivity,
- \f$K\f$ is the bulk modulus,
- \f$\mu\f$ is the shear modulus,.
In order to solve the first equation we will introduce a fictisionus time \f$ \tau \f$:
\f[ \partial_\tau u = (K+\mu/3)\nabla_j \nabla_k u_k+ \mu \Delta u_j - a \nabla_j p + (\rho_s-\rho_f)g_j \f]
\f[ a\partial_t\nabla_i u_i+ \frac{1}{S} \partial_t p = \nabla_i k \nabla_i p, \f]
This implementation is aided to improove stability for incompressible materials.
Let us reformulate the first eqation by the following way:
\f[ \partial_\tau u = \mu \Delta u_j - \nabla_j \tilde p - a\nabla_j p + (\rho_s-\rho_f)g_j, \f]
\f[ \tilde p = -(K+\mu/3) \nabla_k u_k. \f]
then the second eqation will turn to:
\f[ \partial_t p = \frac{aS}{K+\mu/3}\partial_t\tilde p + S\nabla_i k \nabla_i p, \f]
The equation for \f$\tilde p\f$ leads to an instability for low values of the Poisson ration (for nearly incompressible materials).
Therefore the last equation we would like to replace by a relaxation one:
\f[ \partial_\tau \tilde p = - w\left( (K+\mu/3)\nabla_k u_k + \tilde p\right), \f]
where \f$ w \f$ is a relaxation parameter.
Let's return to single time variable. Two time steps are related as follows
\f$ \delta_t = N \delta_\tau \f$. Than \f$ \partial_\tau = N\partial_t \f$.
The resulting set of equations is:
\f[ \partial_\tau u = \mu \Delta u_j - \nabla_j \tilde p - a\nabla_j p + (\rho_s-\rho_f)g_j, \f]
\f[ \partial_\tau \tilde p = - w\left( (K+\mu/3)\nabla_k u_k + \tilde p\right), \f]
\f[ \partial_\tau p = \frac{aS}{K+\mu/3}\partial_\tau\tilde p + \frac{Sk}{N}\Delta p, \f]
*/
class FDPoroElasticity: public ElasticityCommonA
{
private:
Data pressureData;
Data pressureInternalData;
Data pressureLiquidData;
Data pressureLiquidInternalData;
Param hydraulicCondactivity;
Param compresibility;
Param nSubsteps;
public:
FDPoroElasticity();
/**
\param d is a displacement field
\param pl is a pressure of liquid field
\param bM is the bulk modulus
\param sM is the shear modulus
\param k is hydraulic conductivity
\param vT is a vector template
*/
FDPoroElasticity(Data d, Data pl, Param bM,
Param sM, Param k,
const VectorTemplate* vT);
~FDPoroElasticity();
virtual void init();
virtual void execute();
inline Data getPressureData() const;
inline Data getLiquidPressureData() const;
/// defaul value 10
void setNSubsteps(unsigned int n);
};
typedef std::shared_ptr<FDPoroElasticity> SPFDPoroElasticity;
/**
\param d is a displacement field
\param pl is a pressure of liquid field
\param bM is the bulk modulus
\param sM is the shear modulus
\param k is hydraulic conductivity
\param vT is a vector template
*/
SPFDPoroElasticity generateFDPoroElasticity(SPDataWithGhostNodesACLData d,
SPDataWithGhostNodesACLData pl,
double bM,
double sM,
double k,
const VectorTemplate* vT);
/**
\param d is a displacement field
\param pl is a pressure of liquid field
\param bM is the bulk modulus
\param sM is the shear modulus
\param k is hydraulic conductivity
\param vT is a vector template
*/
template <typename T>
SPFDPoroElasticity generateFDPoroElasticity(SPDataWithGhostNodesACLData d,
SPDataWithGhostNodesACLData pl,
UValue<T> bM,
UValue<T> sM,
UValue<T> k,
const VectorTemplate* vT);
//-------------------------IMPLEMENTATION------------------------
inline ElasticityCommonA::Data FDPoroElasticity::getPressureData() const
{
return pressureData;
}
inline ElasticityCommonA::Data FDPoroElasticity::getLiquidPressureData() const
{
return pressureLiquidData;
}
} // asl
#endif // ASLFDELASTICITY_H
|