OilfieldToolsNet/IFEM
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Added axisymmetric linear elasticity integrand

Browse Source 
git-svn-id: http://svn.sintef.no/trondheim/IFEM/trunk@921 e10b68d5-8a6e-419e-a041-bce267b0401d
This commit is contained in:
kmo
2011-04-30 16:06:02 +00:00
committed by Knut Morten Okstad
parent c00aba37a7
commit 5a532e2c61
13 changed files with 493 additions and 39 deletions
Download Patch File Download Diff File Expand all files Collapse all files

15
Apps/LinearElasticity/SIMLinEl2D.C
View File

@@ -14,6 +14,7 @@
#include "SIMLinEl2D.h"
#include "LinIsotropic.h"
#include "LinearElasticity.h"
#include "AxSymmElasticity.h"
#include "AnalyticSolutions.h"
#include "Functions.h"
#include "Utilities.h"
@@ -25,19 +26,25 @@
bool SIMLinEl2D::planeStrain = false;
bool SIMLinEl2D::axiSymmetry = false;
SIMLinEl2D::SIMLinEl2D (int)
{
myProblem = new LinearElasticity(2);
if (axiSymmetry)
myProblem = new AxSymmElasticity();
else
myProblem = new LinearElasticity(2);
}
SIMLinEl2D::~SIMLinEl2D ()
{
for (size_t i = 0; i < mVec.size(); i++)
delete mVec[i];
}
bool SIMLinEl2D::parse (char* keyWord, std::istream& is)
{
char* cline = 0;
@@ -69,7 +76,7 @@ bool SIMLinEl2D::parse (char* keyWord, std::istream& is)
double E = atof(strtok(NULL," "));
double nu = atof(strtok(NULL," "));
double rho = atof(strtok(NULL," "));
mVec.push_back(new LinIsotropic(E,nu,rho,!planeStrain));
mVec.push_back(new LinIsotropic(E,nu,rho,!planeStrain,axiSymmetry));
if (myPid == 0)
std::cout <<"\tMaterial code "<< code <<": "
<< E <<" "<< nu <<" "<< rho << std::endl;
@@ -233,7 +240,7 @@ bool SIMLinEl2D::parse (char* keyWord, std::istream& is)
{
std::cout <<"\tMaterial for all patches: "
<< E <<" "<< nu <<" "<< rho << std::endl;
mVec.push_back(new LinIsotropic(E,nu,rho,!planeStrain));
mVec.push_back(new LinIsotropic(E,nu,rho,!planeStrain,axiSymmetry));
}
else
{
@@ -245,7 +252,7 @@ bool SIMLinEl2D::parse (char* keyWord, std::istream& is)
std::cout <<"\tMaterial for P"<< patch
<<": "<< E <<" "<< nu <<" "<< rho << std::endl;
myProps.push_back(Property(Property::MATERIAL,mVec.size(),pid,3));
mVec.push_back(new LinIsotropic(E,nu,rho,!planeStrain));
mVec.push_back(new LinIsotropic(E,nu,rho,!planeStrain,axiSymmetry));
}
}

3
Apps/LinearElasticity/SIMLinEl2D.h
View File

@@ -51,7 +51,8 @@ protected:
virtual bool initNeumann(size_t propInd);
public:
static bool planeStrain; //!< Plane strain/stress option
static bool planeStrain; //!< Plane strain/stress option
static bool axiSymmetry; //!< Axisymmtry option
protected:
std::vector<Material*> mVec; //!< Material data

24
Apps/LinearElasticity/Test/Cylinder-Axisymm.inp Normal file
View File

@@ -0,0 +1,24 @@
# $Id$
PATCHFILE rectangle.g2
RAISEORDER 1
# patch ru rv
1 2 2
REFINE 1
# patch ru rv
1 3 1
CONSTRAINTS 2
# patch edge code
1 3 2
1 4 2
ISOTROPIC 1
# code E nu rho
0 2.05e11 0.27 7850.0
PRESSURE 1
# patch edge direction pressure
1 1 0 -1.0e6

10
Apps/LinearElasticity/Test/rectangle.g2 Normal file
View File

@@ -0,0 +1,10 @@
200 1 0 0
2 0
2 2
0 0 1 1
2 2
0 0 1 1
1.0 0.0
2.0 0.0
1.0 5.0
2.0 5.0

2
Apps/LinearElasticity/Test/run.sh
View File

@@ -45,6 +45,8 @@ run Cylinder-NURBS.inp -checkRHS -vtf 1 -nviz 5
run Cylinder-Lagrange.inp -checkRHS -vtf 1 -lagrange
run Cylinder-Spectral.inp -checkRHS -vtf 1 -spectral
run Cylinder-Axisymm.inp -2Daxisymm -vtf 1 -nviz 5
eigOpt="-eig 5 -nev 20 -ncv 40 -shift 2.0e-4"
run SquarePlate-Splines.inp $eigOpt -nGauss 3 -vtf 1 -nviz 3
run SquarePlate-Lagrange.inp $eigOpt -nGauss 3 -vtf 1 -lagrange

15
Apps/LinearElasticity/main_LinEl3D.C
View File

@@ -52,6 +52,7 @@
\arg -fixDup : Resolve co-located nodes by merging them into a single node
\arg -2D : Use two-parametric simulation driver (plane stress)
\arg -2Dpstrain : Use two-parametric simulation driver (plane strain)
\arg -2Daxisymm : Use two-parametric simulation driver (axi-symmetric solid)
\arg -lag : Use Lagrangian basis functions instead of splines/NURBS
\arg -spec : Use Spectral basis functions instead of splines/NURBS
*/
@@ -137,8 +138,10 @@ int main (int argc, char** argv)
vizRHS = true;
else if (!strcmp(argv[i],"-fixDup"))
fixDup = true;
else if (!strcmp(argv[i],"-2Dpstrain"))
else if (!strncmp(argv[i],"-2Dpstra",8))
twoD = SIMLinEl2D::planeStrain = true;
else if (!strncmp(argv[i],"-2Daxi",6))
twoD = SIMLinEl2D::axiSymmetry = true;
else if (!strncmp(argv[i],"-2D",3))
twoD = true;
else if (!strncmp(argv[i],"-lag",4))
@@ -153,11 +156,11 @@ int main (int argc, char** argv)
if (!infile)
{
std::cout <<"usage: "<< argv[0]
<<" <inputfile> [-dense|-spr|-superlu<nt>|-samg|-petsc]\n"
<<" [-free] [-lag] [-spec] [-2D[pstrain]] [-nGauss <n>]\n"
<<" [-vtf <format>] [-nviz <nviz>]"
<<" [-nu <nu>] [-nv <nv>] [-nw <nw>]\n"
<<" [-eig <iop>] [-nev <nev>] [-ncv <ncv] [-shift <shf>]\n"
<<" <inputfile> [-dense|-spr|-superlu[<nt>]|-samg|-petsc]\n "
<<" [-free] [-lag] [-spec] [-2D[pstrain|axisymm]] [-nGauss <n>]\n"
<<" [-vtf <format> [-nviz <nviz>]"
<<" [-nu <nu>] [-nv <nv>] [-nw <nw>]]\n"
<<" [-eig <iop> [-nev <nev>] [-ncv <ncv] [-shift <shf>]]\n"
<<" [-ignore <p1> <p2> ...] [-fixDup]"
<<" [-checkRHS] [-check] [-dumpASC]\n";
return 0;

2
src/Integrands/AnaSol.h
View File

@@ -40,7 +40,7 @@ public:
//! \param[in] v3 Secondary solution field, symmetric stress tensor
//!
//! \details It is assumed that all the arguments are pointers to dynamically
//! allocation objects, as the class destructor will attempt to delete them.
//! allocated objects, as the class destructor will attempt to delete them.
AnaSol(RealFunc* s1 = 0, VecFunc* s2 = 0,
VecFunc* v1 = 0, TensorFunc* v2 = 0, STensorFunc* v3 = 0)
: scalSol(s1), scalSecSol(s2), vecSol(v1), vecSecSol(v2), stressSol(v3) {}

268
src/Integrands/AxSymmElasticity.C Normal file
View File

@@ -0,0 +1,268 @@
// $Id$
//==============================================================================
//!
//! \file AxSymmElasticity.C
//!
//! \date Apr 28 2011
//!
//! \author Knut Morten Okstad / SINTEF
//!
//! \brief Integrand implementations for axial-symmetric elasticity problems.
//!
//==============================================================================
#include "AxSymmElasticity.h"
#include "FiniteElement.h"
#include "MaterialBase.h"
#include "Utilities.h"
#include "Tensor.h"
#include "Vec3Oper.h"
#ifndef epsR
//! \brief Zero tolerance for the radial coordinate.
#define epsR 1.0e-16
#endif
AxSymmElasticity::AxSymmElasticity () : Elasticity(2)
{
// Only the current solution is needed
primsol.resize(1);
}
void AxSymmElasticity::print (std::ostream& os) const
{
std::cout <<"Axial-symmetric Elasticity problem\n";
this->Elasticity::print(os);
}
/*!
The strain-displacement matrix for an axiallly symmetric 3D continuum element
is formally defined as:
\f[
[B] = \left[\begin{array}{cc}
\frac{\partial}{\partial r} & 0 \\
0 & \frac{\partial}{\partial z} \\
\frac{1}{r} & 0 \\
\frac{\partial}{\partial z} & \frac{\partial}{\partial r}
\end{array}\right] [N]
\f]
where
[\a N ] is the element basis functions arranged in a [2][2*NENOD] matrix.
*/
bool AxSymmElasticity::kinematics (const Vector& N, const Matrix& dNdX,
double r, SymmTensor& eps) const
{
const size_t nenod = N.size();
Bmat.resize(8,nenod,true);
if (dNdX.cols() < 2)
{
std::cerr <<" *** AxSymmElasticity::kinematics: Invalid dimension on dNdX, "
<< dNdX.rows() <<"x"<< dNdX.cols() <<"."<< std::endl;
return false;
}
else if (r < -epsR)
{
std::cerr <<" *** AxSymmElasticity::kinematics: Invalid point r < 0, "
<< r << std::endl;
return false;
}
#define INDEX(i,j) i+4*(j-1)
// Strain-displacement matrix for 3D axisymmetric elements:
//
// | d/dr 0 |
// [B] = | 0 d/dz | * [N]
// | 1/r 0 |
// | d/dz d/dr |
for (size_t i = 1; i <= nenod; i++)
{
// Normal strain part
Bmat(INDEX(1,1),i) = dNdX(i,1);
Bmat(INDEX(2,2),i) = dNdX(i,2);
// Hoop strain part
Bmat(INDEX(3,1),i) = r <= epsR ? dNdX(i,1) : N(i)/r;
// Shear strain part
Bmat(INDEX(4,1),i) = dNdX(i,2);
Bmat(INDEX(4,2),i) = dNdX(i,1);
}
Bmat.resize(4,2*nenod);
// Evaluate the strains
if (eV && !eV->empty() && eps.dim() > 0)
return Bmat.multiply(*eV,eps); // eps = B*eV
return true;
}
bool AxSymmElasticity::evalInt (LocalIntegral*& elmInt, const FiniteElement& fe,
const Vec3& X) const
{
bool lHaveStrains = false;
SymmTensor eps(2,true), sigma(2,true);
if (eKm || eKg || iS)
{
// Compute the strain-displacement matrix B from N, dNdX and r = X.x,
// and evaluate the symmetric strain tensor if displacements are available
if (!this->kinematics(fe.N,fe.dNdX,X.x,eps)) return false;
if (!eps.isZero(1.0e-16)) lHaveStrains = true;
// Evaluate the constitutive matrix and the stress tensor at this point
double U;
if (!material->evaluate(Cmat,sigma,U,X,eps,eps))
return false;
}
// Axi-symmetric integration point volume; 2*pi*r*|J|*w
const double detJW = 2.0*M_PI*X.x*fe.detJxW;
if (eKm)
{
// Integrate the material stiffness matrix
CB.multiply(Cmat,Bmat).multiply(detJW); // CB = C*B*|J|*w
eKm->multiply(Bmat,CB,true,false,true); // EK += B^T * CB
}
if (eKg && lHaveStrains)
// Integrate the geometric stiffness matrix
this->formKG(*eKg,fe.dNdX,sigma,detJW);
if (eM)
// Integrate the mass matrix
this->formMassMatrix(*eM,fe.N,X,detJW);
if (iS && lHaveStrains)
{
// Integrate the internal forces
sigma *= -detJW;
if (!Bmat.multiply(sigma,*iS,true,true)) // ES -= B^T*sigma
return false;
}
if (eS)
// Integrate the load vector due to gravitation and other body forces
this->formBodyForce(*eS,fe.N,X,detJW);
return this->getIntegralResult(elmInt);
}
bool AxSymmElasticity::evalBou (LocalIntegral*& elmInt, const FiniteElement& fe,
const Vec3& X, const Vec3& normal) const
{
if (!tracFld)
{
std::cerr <<" *** AxSymmElasticity::evalBou: No tractions."<< std::endl;
return false;
}
else if (!eS)
{
std::cerr <<" *** AxSymmElasticity::evalBou: No load vector."<< std::endl;
return false;
}
// Axi-symmetric integration point volume; 2*pi*r*|J|*w
const double detJW = 2.0*M_PI*X.x*fe.detJxW;
// Evaluate the surface traction
Vec3 T = (*tracFld)(X,normal);
// Store the traction value for vizualization
if (!T.isZero()) tracVal[X] = T;
// Integrate the force vector
Vector& ES = *eS;
for (size_t a = 1; a <= fe.N.size(); a++)
{
ES(2*a-1) += T.x*fe.N(a)*detJW;
ES(2*a ) += T.y*fe.N(a)*detJW;
}
return this->getIntegralResult(elmInt);
}
bool AxSymmElasticity::evalSol (Vector& s, const Vector& N,
const Matrix& dNdX, const Vec3& X,
const std::vector<int>& MNPC) const
{
int ierr = 0;
if (eV && !primsol.front().empty())
if ((ierr = utl::gather(MNPC,2,primsol.front(),*eV)))
{
std::cerr <<" *** AxSymmElasticity::evalSol: Detected "
<< ierr <<" node numbers out of range."<< std::endl;
return false;
}
// Evaluate the strain tensor, eps
SymmTensor eps(2,true);
if (!this->kinematics(N,dNdX,X.x,eps))
return false;
// Calculate the stress tensor through the constitutive relation
SymmTensor sigma(2,true); double U;
if (!material->evaluate(Cmat,sigma,U,X,eps,eps))
return false;
// Congruence transformation to local coordinate system at current point
if (locSys) sigma.transform(locSys->getTmat(X));
s = sigma;
s.push_back(sigma.vonMises());
return true;
}
bool AxSymmElasticity::evalSol (Vector& s, const Vector& N,
const Matrix& dNdX, const Vec3& X) const
{
if (!eV || eV->empty())
{
std::cerr <<" *** AxSymmElasticity::evalSol: No displacement vector."
<< std::endl;
return false;
}
else if (eV->size() != dNdX.rows()*2)
{
std::cerr <<" *** AxSymmElasticity::evalSol: Invalid displacement vector."
<<"\n size(eV) = "<< eV->size() <<" size(dNdX) = "
<< dNdX.rows() <<","<< dNdX.cols() << std::endl;
return false;
}
// Evaluate the strain tensor, eps
SymmTensor eps(2,true);
if (!this->kinematics(N,dNdX,X.x,eps))
return false;
// Calculate the stress tensor through the constitutive relation
SymmTensor sigma(2,true); double U;
if (!material->evaluate(Cmat,sigma,U,X,eps,eps))
return false;
s = sigma;
return true;
}
const char* AxSymmElasticity::getFieldLabel (size_t i, const char* prefix) const
{
static const char* s[5] = { "s_r","s_z","s_t","s_zr", "von Mises stress" };
static std::string label;
if (prefix)
label = std::string(prefix) + " " + std::string(s[i]);
else
label = s[i];
return label.c_str();
}

99
src/Integrands/AxSymmElasticity.h Normal file
View File

@@ -0,0 +1,99 @@
// $Id$
//==============================================================================
//!
//! \file AxSymmElasticity.h
//!
//! \date Apr 28 2011
//!
//! \author Knut Morten Okstad / SINTEF
//!
//! \brief Integrand implementations for axial-symmetric elasticity problems.
//!
//==============================================================================
#ifndef _AX_SYMM_ELASTICITY_H
#define _AX_SYMM_ELASTICITY_H
#include "Elasticity.h"
/*!
\brief Class representing the integrand of axial-symmetric elasticity problem.
\details Most methods of this class are inherited form the base class.
Only the \a evalInt, \a evalBou and \a evalSol methods, which are specific
for axial-symmetric linear elasticity problems are implemented here.
*/
class AxSymmElasticity : public Elasticity
{
public:
//! \brief Default constructor.
AxSymmElasticity();
//! \brief Empty destructor.
virtual ~AxSymmElasticity() {}
//! \brief Prints out the problem definition to the given output stream.
virtual void print(std::ostream& os) const;
//! \brief Evaluates the integrand at an interior point.
//! \param elmInt The local integral object to receive the contributions
//! \param[in] fe Finite element data of current integration point
//! \param[in] X Cylindric coordinates of current integration point
virtual bool evalInt(LocalIntegral*& elmInt, const FiniteElement& fe,
const Vec3& X) const;
//! \brief Evaluates the integrand at a boundary point.
//! \param elmInt The local integral object to receive the contributions
//! \param[in] fe Finite element data of current integration point
//! \param[in] X Cylindric coordinates of current integration point
//! \param[in] normal Boundary normal vector at current integration point
virtual bool evalBou(LocalIntegral*& elmInt, const FiniteElement& fe,
const Vec3& X, const Vec3& normal) const;
//! \brief Evaluates the secondary solution at a result point.
//! \param[out] s Array of solution field values at current point
//! \param[in] N Basis function values at current point
//! \param[in] dNdX Basis function gradients at current point
//! \param[in] X Cylindric coordinates of current point
//! \param[in] MNPC Nodal point correspondance for the basis function values
virtual bool evalSol(Vector& s, const Vector& N, const Matrix& dNdX,
const Vec3& X, const std::vector<int>& MNPC) const;
//! \brief Evaluates the finite element (FE) solution at an integration point.
//! \param[out] s The FE stress values at current point
//! \param[in] N Basis function values at current point
//! \param[in] dNdX Basis function gradients at current point
//! \param[in] X Cylindric coordinates of current point
virtual bool evalSol(Vector& s, const Vector& N, const Matrix& dNdX,
const Vec3& X) const;
//! \brief Returns the number of secondary solution field components.
virtual size_t getNoFields() const { return 5; }
//! \brief Returns the name of a secondary solution field component.
//! \param[in] i Field component index
//! \param[in] prefix Name prefix for all components
virtual const char* getFieldLabel(size_t i, const char* prefix = 0) const;
//! \brief Returns \e true if this is an axial-symmetric problem.
virtual bool isAxiSymmetric() const { return true; }
protected:
//! \brief Calculates some kinematic quantities at current point.
//! \param[in] N Basis function values at current point
//! \param[in] dNdX Basis function gradients at current point
//! \param[in] r Radial coordinate of current point
//! \param[out] eps Strain tensor at current point
//!
//! \details The strain displacement matrix \b B is established
//! and stored in the mutable class member \a Bmat.
bool kinematics(const Vector& N, const Matrix& dNdX,
double r, SymmTensor& eps) const;
private:
// Work array declared as member to avoid frequent re-allocation
// within the numerical integration loop (for reduced overhead)
mutable Matrix CB; //!< Result of the matrix-matrix product C*B
};
#endif

27
src/Integrands/Elasticity.C
View File

@@ -470,7 +470,8 @@ bool Elasticity::evalSol (Vector& s, const Vector&,
}
bool Elasticity::evalSol (Vector& s, const Matrix& dNdX, const Vec3& X) const
bool Elasticity::evalSol (Vector& s, const Vector&,
const Matrix& dNdX, const Vec3& X) const
{
if (!eV || eV->empty())
{
@@ -601,13 +602,17 @@ bool ElasticityNorm::evalInt (LocalIntegral*& elmInt, const FiniteElement& fe,
// Evaluate the finite element stress field
Vector sigmah, sigma;
if (!problem.evalSol(sigmah,fe.dNdX,X))
if (!problem.evalSol(sigmah,fe.N,fe.dNdX,X))
return false;
else if (sigmah.size() == 4)
else if (sigmah.size() == 4 && Cinv.rows() == 3)
sigmah.erase(sigmah.begin()+2); // Remove the sigma_zz if plane strain
double detJW = fe.detJxW;
if (problem.isAxiSymmetric())
detJW *= 2.0*M_PI*X.x;
// Integrate the energy norm a(u^h,u^h)
pnorm[0] += sigmah.dot(Cinv*sigmah)*fe.detJxW;
pnorm[0] += sigmah.dot(Cinv*sigmah)*detJW;
if (problem.haveLoads())
{
@@ -616,21 +621,21 @@ bool ElasticityNorm::evalInt (LocalIntegral*& elmInt, const FiniteElement& fe,
// Evaluate the displacement field
Vec3 u = problem.evalSol(fe.N);
// Integrate the external energy (f,u^h)
pnorm[1] += f*u*fe.detJxW;
pnorm[1] += f*u*detJW;
}
if (anasol)
{
// Evaluate the analytical stress field
sigma = (*anasol)(X);
if (sigma.size() == 4)
if (sigma.size() == 4 && Cinv.rows() == 3)
sigma.erase(sigma.begin()+2); // Remove the sigma_zz if plane strain
// Integrate the energy norm a(u,u)
pnorm[2] += sigma.dot(Cinv*sigma)*fe.detJxW;
pnorm[2] += sigma.dot(Cinv*sigma)*detJW;
// Integrate the error in energy norm a(u-u^h,u-u^h)
sigma -= sigmah;
pnorm[3] += sigma.dot(Cinv*sigma)*fe.detJxW;
pnorm[3] += sigma.dot(Cinv*sigma)*detJW;
}
return true;
@@ -649,7 +654,11 @@ bool ElasticityNorm::evalBou (LocalIntegral*& elmInt, const FiniteElement& fe,
// Evaluate the displacement field
Vec3 u = problem.evalSol(fe.N);
double detJW = fe.detJxW;
if (problem.isAxiSymmetric())
detJW *= 2.0*M_PI*X.x;
// Integrate the external energy
pnorm[1] += T*u*fe.detJxW;
pnorm[1] += T*u*detJW;
return true;
}

10
src/Integrands/Elasticity.h
View File

@@ -81,15 +81,16 @@ public:
//! \param[in] dNdX Basis function gradients at current point
//! \param[in] X Cartesian coordinates of current point
//! \param[in] MNPC Nodal point correspondance for the basis function values
virtual bool evalSol(Vector& s,
const Vector& N, const Matrix& dNdX,
virtual bool evalSol(Vector& s, const Vector& N, const Matrix& dNdX,
const Vec3& X, const std::vector<int>& MNPC) const;
//! \brief Evaluates the finite element (FE) solution at an integration point.
//! \param[out] s The FE stress values at current point
//! \param[in] N Basis function values at current point
//! \param[in] dNdX Basis function gradients at current point
//! \param[in] X Cartesian coordinates of current point
virtual bool evalSol(Vector& s, const Matrix& dNdX, const Vec3& X) const;
virtual bool evalSol(Vector& s, const Vector& N, const Matrix& dNdX,
const Vec3& X) const;
//! \brief Evaluates the analytical solution at an integration point.
//! \param[out] s The analytical stress values at current point
@@ -133,6 +134,9 @@ public:
//! \param[in] prefix Name prefix for all components
virtual const char* getFieldLabel(size_t i, const char* prefix = 0) const;
//! \brief Returns \e true if this is an axial-symmetric problem.
virtual bool isAxiSymmetric() const { return false; }
protected:
//! \brief Calculates some kinematic quantities at current point.
//! \param[in] dNdX Basis function gradients at current point

39
src/Integrands/LinIsotropic.C
View File

@@ -15,7 +15,7 @@
#include "Tensor.h"
LinIsotropic::LinIsotropic (bool ps) : planeStress(ps)
LinIsotropic::LinIsotropic (bool ps, bool ax) : planeStress(ps), axiSymmetry(ax)
{
// Default material properties - typical values for steel (SI units)
Emod = 2.05e11;
@@ -27,7 +27,9 @@ LinIsotropic::LinIsotropic (bool ps) : planeStress(ps)
void LinIsotropic::print (std::ostream& os) const
{
std::cout <<"LinIsotropic: ";
if (planeStress)
if (axiSymmetry)
std::cout <<"axial-symmetric, ";
else if (planeStress)
std::cout <<"plane stress, ";
std::cout <<"E = "<< Emod <<", nu = "<< nu <<", rho = "<< rho << std::endl;
}
@@ -51,6 +53,14 @@ void LinIsotropic::print (std::ostream& os) const
0 & 0 & \frac{1}{2}-\nu
\end{array}\right] \f]
For 3D axisymmetric solids: \f[
[C] = \frac{E}{(1+\nu)(1-2\nu)} \left[\begin{array}{cccc}
1-\nu & \nu & \nu & 0 \\
\nu & 1-\nu & \nu & 0 \\
\nu & \nu & 1-\nu & 0 \\
0 & 0 & 0 & \frac{1}{2}-\nu
\end{array}\right] \f]
For 3D: \f[
[C] = \frac{E}{(1+\nu)(1-2\nu)} \left[\begin{array}{cccccc}
1-\nu & \nu & \nu & 0 & 0 & 0 \\
@@ -67,7 +77,7 @@ bool LinIsotropic::evaluate (Matrix& C, SymmTensor& sigma, double& U,
char iop, const TimeDomain*) const
{
const size_t nsd = sigma.dim();
const size_t nst = nsd*(nsd+1)/2;
const size_t nst = nsd == 2 && axiSymmetry ? 4 : nsd*(nsd+1)/2;
C.resize(nst,nst,true);
if (nsd == 1)
@@ -90,7 +100,7 @@ bool LinIsotropic::evaluate (Matrix& C, SymmTensor& sigma, double& U,
}
if (iop < 0) // The inverse C-matrix is wanted
if (nsd == 3 || (nsd == 2 && planeStress))
if (nsd == 3 || (nsd == 2 && (planeStress || axiSymmetry)))
{
C(1,1) = 1.0 / Emod;
C(2,1) = -nu / Emod;
@@ -102,12 +112,12 @@ bool LinIsotropic::evaluate (Matrix& C, SymmTensor& sigma, double& U,
}
else
if (nsd == 2 && planeStress)
if (nsd == 2 && planeStress && !axiSymmetry)
{
C(1,1) = Emod / (1.0 - nu*nu);
C(2,1) = C(1,1) * nu;
}
else // 2D plain strain or 3D
else // 2D plain strain, axisymmetric or 3D
{
double fact = Emod / ((1.0 + nu) * (1.0 - nu - nu));
C(1,1) = fact * (1.0 - nu);
@@ -120,7 +130,16 @@ bool LinIsotropic::evaluate (Matrix& C, SymmTensor& sigma, double& U,
const double G = Emod / (2.0 + nu + nu);
C(nsd+1,nsd+1) = iop < 0 ? 1.0 / G : G;
if (nsd > 2)
if (nsd == 2 && axiSymmetry)
{
C(4,4) = C(3,3);
C(3,1) = C(2,1);
C(3,2) = C(2,1);
C(1,3) = C(2,1);
C(2,3) = C(2,1);
C(3,3) = C(1,1);
}
else if (nsd > 2)
{
C(3,1) = C(2,1);
C(3,2) = C(2,1);
@@ -134,11 +153,11 @@ bool LinIsotropic::evaluate (Matrix& C, SymmTensor& sigma, double& U,
if (iop > 0)
{
// Calculate the stress tensor, sigma = C*eps
Vector sig; // Use a local variable to avoid a redim of sigma
Vector sig; // Use a local variable to avoid redimensioning of sigma
if (eps.dim() != sigma.dim())
{
// Account for non-matching tensor dimensions
SymmTensor epsil(sigma.dim());
SymmTensor epsil(sigma.dim(), nsd == 2 && axiSymmetry);
if (!C.multiply(epsil=eps,sig))
return false;
}
@@ -147,7 +166,7 @@ bool LinIsotropic::evaluate (Matrix& C, SymmTensor& sigma, double& U,
return false;
sigma = sig; // Add sigma_zz in case of plane strain
if (!planeStress && nsd == 2 && sigma.size() == 4)
if (!planeStress && ! axiSymmetry && nsd == 2 && sigma.size() == 4)
sigma(3,3) = nu * (sigma(1,1)+sigma(2,2));
}

18
src/Integrands/LinIsotropic.h
View File

@@ -18,7 +18,7 @@
/*!
\brief Class representing a isotropic linear elastic material model.
\brief Class representing an isotropic linear elastic material model.
*/
class LinIsotropic : public Material
@@ -26,10 +26,17 @@ class LinIsotropic : public Material
public:
//! \brief Default constructor.
//! \param[in] ps If \e true, assume plane stress in 2D
LinIsotropic(bool ps = false);
//! \param[in] ax If \e true, assume 3D axi-symmetric material
LinIsotropic(bool ps = false, bool ax = false);
//! \brief Constructor initializing the material parameters.
LinIsotropic(double E, double v = 0.0, double density = 0.0, bool ps = false)
: Emod(E), nu(v), rho(density), planeStress(ps) {}
//! \param[in] E Young's modulus
//! \param[in] v Poisson's ratio
//! \param[in] density Mass density
//! \param[in] ps If \e true, assume plane stress in 2D
//! \param[in] ax If \e true, assume 3D axi-symmetric material
LinIsotropic(double E, double v = 0.0, double density = 0.0,
bool ps = false, bool ax = false)
: Emod(E), nu(v), rho(density), planeStress(ps), axiSymmetry(ax) {}
//! \brief Empty destructor.
virtual ~LinIsotropic() {}
@@ -37,7 +44,7 @@ public:
virtual bool isPlaneStrain() const { return !planeStress; }
//! \brief Prints out material parameters to the given output stream.
virtual void print(std::ostream&) const;
virtual void print(std::ostream& os) const;
//! \brief Evaluates the mass density at current point.
virtual double getMassDensity(const Vec3&) const { return rho; }
@@ -62,6 +69,7 @@ protected:
double nu; //!< Poisson's ratio
double rho; //!< Mass density
bool planeStress; //!< Plane stress/strain option for 2D problems
bool axiSymmetry; //!< Axi-symmetric option
};
#endif