IFEM-HM/HeatTransfer.C at master · akva2/IFEM-HM · GitHub

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
// $Id$
//==============================================================================
//!
//! \file HeatTransfer.C
//!
//! \date Mar 25 2020
//!
//! \author Arne Morten Kvarving / SINTEF
//!
//! \brief Integrand implementations for heat transfer problems.
//!
//==============================================================================

#include "HeatTransfer.h"
#include "ElmMats.h"
#include "FiniteElement.h"
#include "LocalIntegral.h"
#include "TimeDomain.h"
#include "Utilities.h"
#include "Vec3.h"
#include "Vec3Oper.h"

#include <ext/alloc_traits.h>
#include <iostream>
#include <memory>


HeatTransfer::HeatTransfer (unsigned short int n,
                            bool mass,
                            TimeIntegration::Method method) :
  IntegrandBase(n),
  withMass(mass),
  bdf(TimeIntegration::Order(method))
{
  primsol.resize(TimeIntegration::Steps(method) + 1);
  registerVector("concentration1",&concentration);
}


bool HeatTransfer::initElement (const std::vector<int>& MNPC,
                                LocalIntegral& A)
{
  if (!this->IntegrandBase::initElement(MNPC, A))
      return false;

  int ierr = 0;
  if (withMass) {
    A.vec.resize(A.vec.size()+1);
    if ((ierr = utl::gather(MNPC,1,concentration,A.vec.back())))
      std::cerr <<" *** " << __PRETTY_FUNCTION__ << ":  Detected "
                << ierr <<" node numbers out of range."<< std::endl;
  }

  return ierr == 0;
}


bool HeatTransfer::evalInt (LocalIntegral& elmInt,
                            const FiniteElement& fe,
                            const TimeDomain& time,
                            const Vec3& X) const
{
  ElmMats& elMat = static_cast<ElmMats&>(elmInt);

  double T = elMat.vec[1].dot(fe.N);
  double Tbdf = -T * bdf[1] / time.dt;
  for (int i = 2; i <= bdf.getOrder(); ++i)
    Tbdf -= fe.N.dot(elMat.vec[i]) * bdf[i] / time.dt;

  WeakOps::Mass(elMat.A[0], fe, props.meatDensity() * props.meatHeatCapacity() * bdf[0] / time.dt);
  WeakOps::Laplacian(elMat.A[0], fe, props.meatThermalConductivity());

  if (withMass) {
    Vector dT;
    fe.dNdX.multiply(elMat.vec[1], dT, true);
    Vector dC;
    fe.dNdX.multiply(elMat.vec.back(), dC, true);
    Vec3 waterVelocity = props.waterVelocity(T, dC, dT);
    WeakOps::Advection(elMat.A[0], fe, waterVelocity,
                       props.waterDensity() * props.waterHeatCapacity(),
                       WeakOperators::CONSERVATIVE);
  }

  WeakOps::Source(elMat.b[0], fe, props.meatDensity() * props.meatHeatCapacity() * Tbdf);

  return true;
}


bool HeatTransfer::evalBou(LocalIntegral& elmInt, const FiniteElement& fe,
                           const Vec3& X, const Vec3& normal) const
{
  ElmMats& elMat = static_cast<ElmMats&>(elmInt);

  double T = elMat.vec[1].dot(fe.N);
  WeakOps::Source(elMat.b[0], fe, props.heatTransfer() * (props.ovenTemperature() - T));

  return true;
}


std::string HeatTransfer::getField1Name (size_t, const char* prefix) const
{
  if (!prefix) return "T";

  return prefix + std::string(" T");
}


HeatTransfer::Robin::Robin(unsigned short int n, const HeatTransfer& itg) :
  IntegrandBase(n),
  integrand(itg)
{
}


bool HeatTransfer::Robin::initElementBou(const std::vector<int>& MNPC,
                                                 LocalIntegral& elmInt)
{
  return const_cast<HeatTransfer&>(integrand).initElement(MNPC, elmInt);
}


bool HeatTransfer::Robin::evalBou(LocalIntegral& elmInt, const FiniteElement& fe,
                                          const Vec3& X, const Vec3& normal) const
{
  ElmMats& elMat = static_cast<ElmMats&>(elmInt);
  const HMProperties& props = integrand.getProps();

  if (integrand.enableMass()) {
    double T = elMat.vec[1].dot(fe.N);
    Vector dT;
    fe.dNdX.multiply(elMat.vec[1], dT, true);
    Vector dC;
    fe.dNdX.multiply(elMat.vec.back(), dC, true);
    Vec3 uw = props.waterVelocity(T, dC, dT);
    double val = (1.0 - props.dissipationCoefficient()) * props.heatTransfer();
    WeakOps::Mass(elMat.A[0], fe, (uw * normal) * props.waterHeatCapacity() * props.waterDensity());
    WeakOps::Source(elMat.b[0], fe, val * (props.ovenTemperature() - T));
  } else {
    std::cerr << "No Robin conditions without mass contributions, use Neumann conditions." << std::endl;
    return false;
  }

  return true;
}