# Transient heat transfer problem with FreeFem++

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

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, ny = NN, nz = NN;

// 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

include “CubeByFreeFem.idp”

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 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)(k21
dx(T)dy(TT)+k22dy(T)dy(TT)+k23dz(T)dy(TT))
+int3d(Th,0)(k31
dx(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); */

}