diff --git a/Apps/LinearElasticity/SIMLinEl2D.C b/Apps/LinearElasticity/SIMLinEl2D.C index 862bbf00..3997ed94 100644 --- a/Apps/LinearElasticity/SIMLinEl2D.C +++ b/Apps/LinearElasticity/SIMLinEl2D.C @@ -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)); } } diff --git a/Apps/LinearElasticity/SIMLinEl2D.h b/Apps/LinearElasticity/SIMLinEl2D.h index 879a62d7..19dad203 100644 --- a/Apps/LinearElasticity/SIMLinEl2D.h +++ b/Apps/LinearElasticity/SIMLinEl2D.h @@ -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 mVec; //!< Material data diff --git a/Apps/LinearElasticity/Test/Cylinder-Axisymm.inp b/Apps/LinearElasticity/Test/Cylinder-Axisymm.inp new file mode 100644 index 00000000..e52cf716 --- /dev/null +++ b/Apps/LinearElasticity/Test/Cylinder-Axisymm.inp @@ -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 diff --git a/Apps/LinearElasticity/Test/rectangle.g2 b/Apps/LinearElasticity/Test/rectangle.g2 new file mode 100644 index 00000000..2f53a30b --- /dev/null +++ b/Apps/LinearElasticity/Test/rectangle.g2 @@ -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 diff --git a/Apps/LinearElasticity/Test/run.sh b/Apps/LinearElasticity/Test/run.sh index f9afcc90..6658004d 100755 --- a/Apps/LinearElasticity/Test/run.sh +++ b/Apps/LinearElasticity/Test/run.sh @@ -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 diff --git a/Apps/LinearElasticity/main_LinEl3D.C b/Apps/LinearElasticity/main_LinEl3D.C index 86f2c62f..6b1ae5a5 100644 --- a/Apps/LinearElasticity/main_LinEl3D.C +++ b/Apps/LinearElasticity/main_LinEl3D.C @@ -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] - <<" [-dense|-spr|-superlu|-samg|-petsc]\n" - <<" [-free] [-lag] [-spec] [-2D[pstrain]] [-nGauss ]\n" - <<" [-vtf ] [-nviz ]" - <<" [-nu ] [-nv ] [-nw ]\n" - <<" [-eig ] [-nev ] [-ncv ]\n" + <<" [-dense|-spr|-superlu[]|-samg|-petsc]\n " + <<" [-free] [-lag] [-spec] [-2D[pstrain|axisymm]] [-nGauss ]\n" + <<" [-vtf [-nviz ]" + <<" [-nu ] [-nv ] [-nw ]]\n" + <<" [-eig [-nev ] [-ncv ]]\n" <<" [-ignore ...] [-fixDup]" <<" [-checkRHS] [-check] [-dumpASC]\n"; return 0; diff --git a/src/Integrands/AnaSol.h b/src/Integrands/AnaSol.h index 9aeef942..dc84ad42 100644 --- a/src/Integrands/AnaSol.h +++ b/src/Integrands/AnaSol.h @@ -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) {} diff --git a/src/Integrands/AxSymmElasticity.C b/src/Integrands/AxSymmElasticity.C new file mode 100644 index 00000000..639bc380 --- /dev/null +++ b/src/Integrands/AxSymmElasticity.C @@ -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& 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(); +} diff --git a/src/Integrands/AxSymmElasticity.h b/src/Integrands/AxSymmElasticity.h new file mode 100644 index 00000000..51457662 --- /dev/null +++ b/src/Integrands/AxSymmElasticity.h @@ -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& 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 diff --git a/src/Integrands/Elasticity.C b/src/Integrands/Elasticity.C index 6123af15..110fb45c 100644 --- a/src/Integrands/Elasticity.C +++ b/src/Integrands/Elasticity.C @@ -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; } diff --git a/src/Integrands/Elasticity.h b/src/Integrands/Elasticity.h index 55affb51..ccad0bad 100644 --- a/src/Integrands/Elasticity.h +++ b/src/Integrands/Elasticity.h @@ -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& 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 diff --git a/src/Integrands/LinIsotropic.C b/src/Integrands/LinIsotropic.C index 1b81dd72..a7a3fdd3 100644 --- a/src/Integrands/LinIsotropic.C +++ b/src/Integrands/LinIsotropic.C @@ -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)); } diff --git a/src/Integrands/LinIsotropic.h b/src/Integrands/LinIsotropic.h index 5620fce8..021ce9b4 100644 --- a/src/Integrands/LinIsotropic.h +++ b/src/Integrands/LinIsotropic.h @@ -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