Dear all,
I am willing to solve a basic transient heat transfer problem with Freefelm++. While solving, I am facing the issue of a minimum temperature of 100K, which normally should not be the case with my problem. Here I am explaining my problem along with the script for freefem++. There are one edp and another idp file.
Problem: Cube with imposed temperature (300K) on one face (label 3). Imposed flux on the opposite face (label=4). SOlving basic transient heat transfer problem for 1 iteration with dt=0.1s.
IDP File
load “medit”
load “msh3”
func mesh3 Cube (int[int] &NN, real[int, int] &BB, int[int, int] &L)
{
real x0 = BB(0,0), x1 = BB(0,1);
real y0 = BB(1,0), y1 = BB(1,1);
real z0 = BB(2,0), z1 = BB(2,1);
int nx = NN[0], ny = NN[1], nz = NN[2];
// 2D mesh
mesh Thx = square(nx, ny, [x0+(x1-x0)*x, y0+(y1-y0)*y]); //
// 3D mesh
int[int] rup = [0, L(2,1)], rdown=[0, L(2,0)];
int[int] rmid=[1, L(1,0), 2, L(0,1), 3, L(1,1), 4, L(0,0)];
mesh3 Th = buildlayers(Thx, nz, zbound=[z0,z1], labelmid=rmid, labelup = rup, labeldown = rdown);
return Th;
}
EDP File
FreeFem Code
lload “gmsh”;
load “medit”;
load “iovtk”;
include “CubeByFreeFem.idp”
load “msh3”
int[int] NN = [10,10,10]; //the number of step in each direction
real [int, int] BB = [[1,2],[1,2],[1,2]]; //the bounding box
int [int, int] L = [[1,2],[3,4],[5,6]]; //the label of the 6 face left,right, front,back, top down
mesh3 Th = Cube(NN, BB, L);
medit(“Th”, Th);
savemesh(Th, “Steelmesh.mesh”);
real dt=0.1;
int nbstep =100;//0.5/dt;//1.2/dt; //number of time steps
real Tinit = 300;
macro grad(field)[dx(field),dy(field),dz(field)]//
macro normal [N.x, N.y, N.z]//
macro fluxVECTOR(L) [0, L,0 ]//
fespace Vh(Th,P1);
Vh T,Told,TT; //current temperature, previous temperature, test function
Told = 300;
for (int it=0; it < 1 ; it++)
{
/////// thermal properties /////////
fespace ConductTensor(Th,P0);
ConductTensor k11,k12,k13,k21,k22,k23,k31,k32,k33;
k11 = 0.0171*Told + 1.043;
k12 = 0;
k13 = 0;
k21 = 0;
k22 = 0.0171*Told + 1.043;
k23 = 0;
k31 = 0;
k32 = 0;
k33 =0.0171*Told + 1.043;
real rho=4430;
real Cp= 372;
//////////////////////////////////////////////////////////////////
solve transient(T,TT)=
int3d(Th,0)(k11*dx(T)dx(TT)+k12dy(T)dx(TT)+k13dz(T)dx(TT))
+int3d(Th,0)(k21dx(T)dy(TT)+k22dy(T)dy(TT)+k23dz(T)dy(TT))
+int3d(Th,0)(k31dx(T)dz(TT)+k32dy(T)dz(TT)+k33dz(T)*dz(TT))
+ int3d(Th,0) ((rho*Cp/dt) * T * TT ) + int3d(Th,0) ( - (rho*Cp/dt) * Told * TT) //evolution terms
+int2d(Th,4)( fluxVECTOR(-1.2e8) '* normal*TT )
+on(3,T=300)
;
Told = T;
int[int] Order = [1,1];
/* fespace flux(Th,P0);
flux qx,qy,qz, qtest;
qx=-k11*dx(T);
qy=-k22*dy(T);
qz=-k33*dz(T);
string DataName = “Temp”;
savevtk(“Temp_steel_Cube_Freefem”+it+".vtu", Th, T,[qx,qy,qz], dataname=“Temp HeatFlux”, order=Order); */
}
Thanks in advance for your help.
