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
|
#include "ViscousStress.h"
#include "ATC_Error.h"
#include "LammpsInterface.h"
#include <iostream>
using ATC_Utility::command_line;
using ATC_Utility::str2dbl;
using std::string;
using std::vector;
using std::fstream;
namespace ATC {
//=============================================================================
// isotropic constant viscosity
//=============================================================================
ViscousStressConstant::ViscousStressConstant(fstream &fileId)
: ViscousStress(), viscosity_(0)
{
if (!fileId.is_open()) throw ATC_Error("cannot open material file");
vector<string> line;
while(fileId.good()) {
command_line(fileId, line);
if (line[0] == "end") {
if (viscosity_ < 0.0)
throw ATC_Error("ViscousStressConstant:: bad constant viscosity");
return;
}
else if (line[0]=="viscosity") viscosity_ = str2dbl(line[1]);
else throw ATC_Error( "viscosity constant - unrecognized material function");
}
}
//=============================================================================
// compute the stress at N integration points from the velocity gradients
// T_{ij} = viscosity * du_i/dx_j
//=============================================================================
void ViscousStressConstant::viscous_stress(const FIELD_MATS & /* fields */,
const GRAD_FIELD_MATS &gradFields,
DENS_MAT_VEC &sigma)
{
GRAD_FIELD_MATS::const_iterator du_itr = gradFields.find(VELOCITY);
const DENS_MAT_VEC &du = du_itr->second;
CLON_VEC DuxDx(du[0],CLONE_COL,0);
CLON_VEC DuxDy(du[1],CLONE_COL,0);
CLON_VEC DuxDz(du[2],CLONE_COL,0);
CLON_VEC DuyDx(du[0],CLONE_COL,1);
CLON_VEC DuyDy(du[1],CLONE_COL,1);
CLON_VEC DuyDz(du[2],CLONE_COL,1);
CLON_VEC DuzDx(du[0],CLONE_COL,2);
CLON_VEC DuzDy(du[1],CLONE_COL,2);
CLON_VEC DuzDz(du[2],CLONE_COL,2);
const INDEX N = DuxDx.size(); // # of integration pts
sigma.assign(3, DENS_MAT(N,3));
// precompute the pressure and copy to the diagonal
column(sigma[0],0) = DuxDx;
column(sigma[0],1) = DuyDx;
column(sigma[0],2) = DuzDx;
column(sigma[1],0) = DuxDy;
column(sigma[1],1) = DuyDy;
column(sigma[1],2) = DuzDy;
column(sigma[2],0) = DuxDz;
column(sigma[2],1) = DuyDz;
column(sigma[2],2) = DuzDz;
sigma[0] *= -1.*viscosity_;
sigma[1] *= -1.*viscosity_;
sigma[2] *= -1.*viscosity_;
}
void ViscousStressConstant::viscosity(const FIELD_MATS &fields,
DENS_MAT &coefs) const
{
const DENS_MAT & v = (fields.find(VELOCITY))->second;
coefs.resize(v.nRows(),v.nCols());
coefs = -1.*viscosity_;
}
}// end atc namespace
|
