[Getdp] Heat equation
Matt Koch
mattkoch at alum.mit.edu
Sat Dec 21 01:48:14 CET 2002
Fabio:
To begin with, isn't the heat conduction equation as follows?
rho d(cp T)/dt = +k div(grad(T)) + f(T,t)
In other words, if you set u = cp T, then your version should read as
follows!
rho d(u)/dt = +k/cp div(grad(u)) + f(T,t)
That could be the reason why you are not able to reach steady state
within the time frame you expect.
Anyway, apart from that, I received very good help from the GetDP group
on the topic of instationary heat conduction and was able to write a 3D
problem based on that input. I have attached the .geo file and the .pro
file, and it should be fairly simple to bring these from 3D to 2D.
Regards,
Matt Koch
-------------- next part --------------
// === Matt Koch ==================================== Fri-10-31-2002 ===
// Example from GetDP Mailing List!
// =====================================================================
//
// --- Group -----------------------------------------------------------
//
Group {
SfcInr = Region[54];
SfcOtr = Region[55];
SfcLwr = Region[56];
SfcUpr = Region[57];
Vol = Region[58];
Sfc = Region[{SfcInr,SfcOtr,SfcLwr,SfcUpr}];
Dmn = Region[{Sfc,Vol}];
}
//
// --- Function --------------------------------------------------------
//
Function {
TheCon[Vol] = 50.0;
TheMas[Vol] = 50.0 * 1000.0;
TmeFct[] = ($Time < 50.0e-03) ? 1.0 : 0.0;
PwrSrc[Vol] = 1.0e+06 * TmeFct[];
HTCInr[SfcInr] = 5.0;
HTCOtr[SfcOtr] = 15.0;
TmpInr[SfcInr] = 20.0+273.15;
TmpOtr[SfcOtr] = 30.0+273.15;
TmpLwr[Dmn] = 10.0+273.15;
TmpUpr[Dmn] = 12.0+273.15;
TmeStt = 0.0;
TmeStp = 1.0;
TmeDel = 0.1;
SveFct[] = !($TimeStep % 10);
}
//
// --- Constraint ------------------------------------------------------
//
Constraint {
{Name TmpCst;
Case {
{Region SfcLwr; Value TmpLwr[];}
{Region SfcUpr; Value TmpUpr[];}
}
}
}
//
// --- Jacobian --------------------------------------------------------
//
Jacobian {
{Name TmpJacVol;
Case {
{Region All; Jacobian Vol;}
}
}
{Name TmpJacSfc;
Case {
{Region All; Jacobian Sur;}
}
}
}
//
// --- Integration -----------------------------------------------------
//
Integration {
{Name TmpInt;
Case {
{Type Gauss;
Case {
{GeoElement Point ; NumberOfPoints 1;}
{GeoElement Line ; NumberOfPoints 3;}
{GeoElement Triangle ; NumberOfPoints 4;}
{GeoElement Quadrangle ; NumberOfPoints 4;}
{GeoElement Tetrahedron; NumberOfPoints 4;}
{GeoElement Hexahedron ; NumberOfPoints 6;}
{GeoElement Prism ; NumberOfPoints 6;}
}
}
}
}
}
//
// --- FunctionSpace ---------------------------------------------------
//
FunctionSpace {
{Name TmpSpc; Type Form0;
BasisFunction {
{Name sn; NameOfCoef Tmpn; Function BF_Node; Support Dmn; Entity NodesOf[All];}
}
Constraint {
{NameOfCoef Tmpn; EntityType NodesOf; NameOfConstraint TmpCst;}
}
}
}
//
// --- Formulation -----------------------------------------------------
//
Formulation {
{Name TmpFrm; Type FemEquation;
Quantity {
{Name Tmp; Type Local; NameOfSpace TmpSpc;}
}
Equation {
Galerkin { [ +TheCon[] * Dof{d Tmp}, {d Tmp} ];
In Vol ; Integration TmpInt; Jacobian TmpJacVol;}
Galerkin { Dt [ +TheMas[] * Dof{ Tmp}, { Tmp} ];
In Vol ; Integration TmpInt; Jacobian TmpJacVol;}
Galerkin { [ -PwrSrc[], {Tmp} ];
In Vol ; Integration TmpInt; Jacobian TmpJacVol;}
Galerkin { [ +HTCInr[] * Dof{Tmp}, {Tmp} ];
In SfcInr; Integration TmpInt; Jacobian TmpJacSfc;}
Galerkin { [ -HTCInr[] * TmpInr[], {Tmp} ];
In SfcInr; Integration TmpInt; Jacobian TmpJacSfc;}
Galerkin { [ +HTCOtr[] * Dof{Tmp}, {Tmp} ];
In SfcOtr; Integration TmpInt; Jacobian TmpJacSfc;}
Galerkin { [ -HTCOtr[] * TmpOtr[], {Tmp} ];
In SfcOtr; Integration TmpInt; Jacobian TmpJacSfc;}
}
}
}
//
// --- Resolution ------------------------------------------------------
//
Resolution {
{Name TmpRes;
System {
{Name Tmp; NameOfFormulation TmpFrm;}
}
Operation {
InitSolution Tmp ; SaveSolution Tmp;
TimeLoopTheta {
Time0 TmeStt; TimeMax TmeStp; DTime TmeDel; Theta 1.;
Operation {
Generate Tmp ; Solve Tmp;
Test[SveFct[]] {SaveSolution Tmp;}
}
}
}
}
}
//
// --- PostProcessing --------------------------------------------------
//
PostProcessing {
{Name TmpPst; NameOfFormulation TmpFrm;
Quantity {
{Name Tmp; Value{ Local{[ {Tmp} ]; In Vol;} }}
}
}
}
//
// --- PostOperation ---------------------------------------------------
//
PostOperation {
{Name TmpOps; NameOfPostProcessing TmpPst;
Operation {
Print[Tmp, OnElementsOf Vol, File "Tmp_3D.pos"];
}
}
}
//
// =====================================================================
-------------- next part --------------
// === Matt Koch ==================================== Fri-10-31-2002 ===
// Hollow Cylinder!
// =====================================================================
//
// --- Dimensions ------------------------------------------------------
//
dRdsInr = 26.0e-03;
dRdsOtr = 30.0e-03;
dHgt = 100.0e-03;
dMsh = 2.0e-03;
//
// --- Points ----------------------------------------------------------
//
Point(1) = { 0.0, 0.0, 0.0, dMsh};
Point(2) = {+dRdsInr, 0.0, 0.0, dMsh};
Point(3) = { 0.0, +dRdsInr, 0.0, dMsh};
Point(4) = {-dRdsInr, 0.0, 0.0, dMsh};
Point(5) = { 0.0, -dRdsInr, 0.0, dMsh};
Point(6) = {+dRdsOtr, 0.0, 0.0, dMsh};
Point(7) = { 0.0, +dRdsOtr, 0.0, dMsh};
Point(8) = {-dRdsOtr, 0.0, 0.0, dMsh};
Point(9) = { 0.0, -dRdsOtr, 0.0, dMsh};
//
// --- Circles ---------------------------------------------------------
//
Circle(1) = {2,1,3};
Circle(2) = {3,1,4};
Circle(3) = {4,1,5};
Circle(4) = {5,1,2};
Circle(5) = {6,1,7};
Circle(6) = {7,1,8};
Circle(7) = {8,1,9};
Circle(8) = {9,1,6};
//
// --- Geometry Surfaces -----------------------------------------------
//
LopInr = newreg;
Line Loop(LopInr) = {1,2,3,4};
LopOtr = newreg;
Line Loop(LopOtr) = {5,6,7,8};
SfcGeoBtm = newreg;
Plane Surface(SfcGeoBtm) = {LopOtr,LopInr};
//
// --- Geometry Volumes ------------------------------------------------
//
VolGeo = newreg;
Extrude Surface {SfcGeoBtm, {0.0,0.0,dHgt} } {Layers{20,VolGeo,1.0}; Recombine;};
//
// --- Physics Surfaces ------------------------------------------------
//
SfcPhyInr = newreg; // 54
Physical Surface(SfcPhyInr) = {40,44,48,52};
SfcPhyOtr = newreg; // 55
Physical Surface(SfcPhyOtr) = {24,28,32,36};
SfcPhyLwr = newreg; // 56
Physical Surface(SfcPhyLwr) = {11};
SfcPhyUpr = newreg; // 57
Physical Surface(SfcPhyUpr) = {53};
//
// --- Physics Volumes -------------------------------------------------
//
VolPhy = newreg; // 58
Physical Volume(VolPhy) = {12};
//
// =====================================================================