Added: Thermal expansion coefficient as material property

git-svn-id: http://svn.sintef.no/trondheim/IFEM/trunk@2882 e10b68d5-8a6e-419e-a041-bce267b0401d

LinIsotropic.C

@@ -12,9 +12,12 @@
 //==============================================================================
 
 #include "LinIsotropic.h"
+#include "Utilities.h"
+#include "Functions.h"
 #include "Field.h"
 #include "Tensor.h"
 #include "Vec3.h"
+#include "tinyxml.h"
 
 
 LinIsotropic::LinIsotropic (bool ps, bool ax) : planeStress(ps), axiSymmetry(ax)
@@ -25,31 +28,63 @@ LinIsotropic::LinIsotropic (bool ps, bool ax) : planeStress(ps), axiSymmetry(ax)
   Emod = 2.05e11;
   nu = 0.29;
   rho = 7.85e3;
+  Afunc = NULL;
+  alpha = 1.2e-7;
 }
 
 
 LinIsotropic::LinIsotropic (RealFunc* E, double v, double den, bool ps, bool ax)
-  : Efunc(E), Efield(NULL), nu(v), rho(den), planeStress(ps), axiSymmetry(ax)
+  : Efunc(E), Efield(NULL), nu(v), rho(den), Afunc(NULL), alpha(1.2e-7),
+    planeStress(ps), axiSymmetry(ax)
 {
   Emod = -1.0; // Should not be referenced
 }
 
 
 LinIsotropic::LinIsotropic (Field* E, double v, double den, bool ps, bool ax)
-  : Efunc(NULL), Efield(E), nu(v), rho(den), planeStress(ps), axiSymmetry(ax)
+  : Efunc(NULL), Efield(E), nu(v), rho(den), Afunc(NULL), alpha(1.2e-7),
+    planeStress(ps), axiSymmetry(ax)
 {
   Emod = -1.0; // Should not be referenced
 }
 
 
+void LinIsotropic::parse (const TiXmlElement* elem)
+{
+  if (utl::getAttribute(elem,"E",Emod))
+    std::cout <<" "<< Emod;
+  if (utl::getAttribute(elem,"nu",nu))
+    std::cout <<" "<< nu;
+  if (utl::getAttribute(elem,"rho",rho))
+    std::cout <<" "<< rho;
+  if (utl::getAttribute(elem,"alpha",alpha))
+    std::cout <<" "<< alpha;
+
+  const TiXmlNode* aval = NULL;
+  const TiXmlElement* child = elem->FirstChildElement();
+  for (; child; child = child->NextSiblingElement())
+    if (!strcasecmp(child->Value(),"thermalexpansion"))
+    {
+      std::cout <<" ";
+      std::string type;
+      utl::getAttribute(child,"type",type,true);
+      if ((aval = child->FirstChild()))
+        Afunc = utl::parseTimeFunc(aval->Value(),type);
+    }
+
+  if (!aval) std::cout << std::endl;
+}
+
+
 void LinIsotropic::print (std::ostream& os) const
 {
   std::cout <<"LinIsotropic: ";
   if (axiSymmetry)
     std::cout <<"axial-symmetric, ";
   else if (planeStress)
     std::cout <<"plane stress, ";
-  std::cout <<"E = "<< Emod <<", nu = "<< nu <<", rho = "<< rho << std::endl;
+  std::cout <<"E = "<< Emod <<", nu = "<< nu <<", rho = "<< rho
+            <<"alpha = "<< alpha << std::endl;
 }
 
 
@@ -201,3 +236,9 @@ bool LinIsotropic::evaluate (Matrix& C, SymmTensor& sigma, double& U,
 
   return true;
 }
+
+
+double LinIsotropic::getThermalExpansion (double T) const
+{
+  return Afunc ? (*Afunc)(T) : alpha;
+}

LinIsotropic.h

@@ -39,7 +39,7 @@ class LinIsotropic : public Material
   LinIsotropic(double E, double v = 0.0, double densty = 0.0,
                bool ps = false, bool ax = false)
     : Efunc(NULL), Efield(NULL), Emod(E), nu(v), rho(densty),
-      planeStress(ps), axiSymmetry(ax) {}
+      Afunc(NULL), alpha(0.0), planeStress(ps), axiSymmetry(ax) {}
   //! \brief Constructor initializing the material parameters.
   //! \param[in] E Young's modulus (spatial function)
   //! \param[in] v Poisson's ratio
@@ -57,16 +57,21 @@ class LinIsotropic : public Material
   LinIsotropic(Field* E, double v = 0.0, double density = 0.0,
                bool ps = false, bool ax = false);
   //! \brief The destructor deletes the stiffness function, if defined.
-  virtual ~LinIsotropic() { delete Efunc; delete Efield; }
+  virtual ~LinIsotropic() { delete Efunc; delete Efield; delete Afunc; }
 
-  //! \brief Returns \e false if plane stress in 2D.
-  virtual bool isPlaneStrain() const { return !planeStress; }
+  //! \brief Parses material parementers from an XML element.
+  virtual void parse(const TiXmlElement* elem);
 
   //! \brief Prints out material parameters to the given output stream.
   virtual void print(std::ostream& os) const;
 
+  //! \brief Returns \e false if plane stress in 2D.
+  virtual bool isPlaneStrain() const { return !planeStress; }
+
   //! \brief Evaluates the mass density at current point.
   virtual double getMassDensity(const Vec3&) const { return rho; }
+  //! \brief Evaluates the thermal expansion coefficient for given temperature.
+  virtual double getThermalExpansion(double T) const;
 
   //! \brief Evaluates the constitutive relation at an integration point.
   //! \param[out] C Constitutive matrix at current point
@@ -98,6 +103,8 @@ class LinIsotropic : public Material
   double Emod;        //!< Young's modulus (constant)
   double nu;          //!< Poisson's ratio
   double rho;         //!< Mass density
+  ScalarFunc* Afunc;  //!< Thermal expansion coefficient function
+  double alpha;       //!< Thermal expansion coefficient (constant)
   bool   planeStress; //!< Plane stress/strain option for 2D problems
   bool   axiSymmetry; //!< Axi-symmetric option
 };

MaterialBase.h

@@ -21,6 +21,7 @@ class Tensor;
 class SymmTensor;
 class FiniteElement;
 class Field;
+class TiXmlElement;
 struct TimeDomain;
 
 
@@ -38,12 +39,15 @@ class Material
   //! \brief Empty destructor.
   virtual ~Material() {}
 
-  //! \brief Returns \e false if plane stress in 2D.
-  virtual bool isPlaneStrain() const { return true; }
+  //! \brief Parses material parementers from an XML element.
+  virtual void parse(const TiXmlElement*) {}
 
   //! \brief Prints out material parameters to the given output stream.
   virtual void print(std::ostream&) const {}
 
+  //! \brief Returns \e false if plane stress in 2D.
+  virtual bool isPlaneStrain() const { return true; }
+
   //! \brief Initializes the material with the number of integration points.
   virtual void initIntegration(size_t) {}
   //! \brief Initializes the material model for a new integration loop.
@@ -57,6 +61,8 @@ class Material
 
   //! \brief Evaluates the mass density at current point.
   virtual double getMassDensity(const Vec3&) const { return 0.0; }
+  //! \brief Evaluates the thermal expansion coefficient for given temperature.
+  virtual double getThermalExpansion(double) const { return 0.0; }
 
   //! \brief Evaluates the constitutive relation at an integration point.
   //! \param[out] C Constitutive matrix at current point
@@ -84,9 +90,8 @@ class Material
   virtual int getNoIntVariables() const { return 0; }
   //! \brief Returns an internal variable associated with the material model.
   virtual double getInternalVariable(int, char*, size_t=0) const { return 0.0; }
-
   //! \brief Returns whether the material model has diverged.
-  virtual bool diverged(size_t=0) const { return false; }
+  virtual bool diverged(size_t = 0) const { return false; }
 };
 
 #endif

SIMElasticity.h

@@ -18,11 +18,12 @@
 #include "SIM3D.h"
 #include "Elasticity.h"
 #include "LinIsotropic.h"
+#include "TimeStep.h"
 #include "AnaSol.h"
 #include "Functions.h"
 #include "Utilities.h"
 #include "tinyxml.h"
-#include "TimeStep.h"
+
 
 /*!
   \brief Driver class for isogeometric FEM analysis of elasticity problems.
@@ -58,7 +59,7 @@ template<class Dim> class SIMElasticity : public Dim
     Elasticity* elp = dynamic_cast<Elasticity*>(Dim::myProblem);
     if (elp)
       elp->advanceStep(tp.time.dt,tp.time.dtn);
-    
+
     return true;
   }
 
@@ -344,25 +345,18 @@ template<class Dim> class SIMElasticity : public Dim
       return this->Dim::parse(elem);
 
     const TiXmlElement* child = elem->FirstChildElement();
-    for (; child; child = child->NextSiblingElement()) {
+    for (; child; child = child->NextSiblingElement())
       if (this->parseDimSpecific(child))
         continue;
 
-      if (!strcasecmp(child->Value(),"isotropic")) {
+      else if (!strcasecmp(child->Value(),"isotropic")) {
         int code = this->parseMaterialSet(child,mVec.size());
-
-        double E = 1000.0, nu = 0.3, rho = 1.0;
-        utl::getAttribute(child,"E",E);
-        utl::getAttribute(child,"nu",nu);
-        utl::getAttribute(child,"rho",rho);
-
+        std::cout <<"\tMaterial code "<< code <<":";
         if (Dim::dimension == 2)
-          mVec.push_back(new LinIsotropic(E,nu,rho,!planeStrain,axiSymmetry));
+          mVec.push_back(new LinIsotropic(!planeStrain,axiSymmetry));
         else
-          mVec.push_back(new LinIsotropic(E,nu,rho));
-
-        std::cout <<"\tMaterial code "<< code <<": "
-                  << E <<" "<< nu <<" "<< rho << std::endl;
+          mVec.push_back(new LinIsotropic(false,false));
+        mVec.back()->parse(child);
       }
 
       else if (!strcasecmp(child->Value(),"gravity")) {
@@ -394,8 +388,7 @@ template<class Dim> class SIMElasticity : public Dim
       }
 
       else
-        return this->Dim::parse(elem);
-    }
+        this->Dim::parse(child);
 
     if (!mVec.empty())
       this->getIntegrand()->setMaterial(mVec.front());