Heat Transfer Project

Mathematical Modeling of Convection on a Parabolic Fin

Gerald Condon

July 2015

Introduction


Heat transfer is an integral part of the engineering and design process. All reactions/movements/changes of state causes the distribution of heat to change. This project will investigate the movement of heat in a heat disipation fin using core mathematical principles and finite differencing approaches. The heat disipation fin for this study described bellow:

Figure 1



Upper Fin Edge:   $y_u(x)=\frac{t}{2}\left[\frac{(L-x)}{L}\right]^2$

Lower Fin Edge:   $y_l(x)=\frac{-t}{2}\left[\frac{(L-x)}{L}\right]^2$

Length:   $L=0.1$

Thickness:   $t=0.02$


Thermal Cond:
$k=20$ (W/mK)
Convection rate:
$h=100$ (W/m^2)

Base Temp:
$T_b=120$ (C)
Flow Temp:
$T_\infty=10$ (C)

Theory



For a concave parabolic fin, the profile function is

$f(x)=\frac{\delta}{2}\left[\frac{(x)}{b}\right]^2$


therefore

$\frac{df(x)}{dx}=\frac{\delta}{b}\frac{x}{b}$


when placed within the standard equation for temperature excess:   $\theta(x)=T(x)-T_s$

$x^2\frac{d^2\theta}{dx^2}+2x\frac{d\theta}{dx}-m^2b^2\theta=0$


Where $m=\left(\frac{2h}{k\delta}\right)^\frac{1}{2}$

This equation is known as the Euler equation. To solve such an equation we substitute $x=e^v$ or $v=ln(x)$ as follows:

$\frac{d\theta}{dx}=\frac{1}{x}\frac{d\theta}{dv}$


and

$\frac{d^2\theta}{dx^2}=-\frac{1}{x^2}\frac{d\theta}{dv}+\frac{1}{x}\frac{d\frac{d\theta}{dv}}{dx}$


which becomes

$\frac{d^2\theta}{dx^2}=-\frac{1}{x^2}\frac{d\theta}{dv}+\frac{1}{x^2}\frac{d^2\theta}{dv^2}$


when substituted into the Euler equation

$x^2\left(\frac{1}{x^2}\frac{d^2\theta}{dv^2}-\frac{1}{x^2}\frac{d\theta}{dv}\right)+2x\left(\frac{1}{x}\frac{d\theta}{dv}\right)-m^2b^2\theta=0$


which simplifies to:

$\frac{d^2\theta}{dv^2}+\frac{d\theta}{dv}-m^2b^2\theta=0$


whose solution is of form:

$\theta=C_1e^{\alpha v}+C_2e^{\beta v}$


Thus:

$q_b=\frac{k\delta L \theta_b}{2b}\left[-1+\sqrt[2]{1+(2mb)^2}\right]$


To acomplish efficiency, equation $q_b$ is devided by the ideal heat flow:

$q_{id}=2hbL\theta_b$


Thus:

$\eta=\frac{2}{1+\sqrt[2]{1+(2mb)^2}}$

Numerical Implementation


In comparision to an analytical solution, which allows for temperature determination at any and all points in a medium. Finite differencing allows a problem to be broken down into distinct and discrete points. Treating each particular point as a small fundamentally calculatable problem. Known in the realm of Computer Science as a "Divide and Conquer" technique. The acuracy of the problem is a function of the amount of nodal points and convergance criteria applied.

The Finite-Differencing method used for this project will be derived from Energy Balance Equation for steady state where   $\dot{E_in}+\dot{E_g}=0$ meaning that the net Energy rate entering a node and being generated in a node is equivalent to zero. In a 2 dymentional case where there exists only 4 possible enterences in a node ( $\uparrow,\downarrow,\rightarrow,\leftarrow$ ) the equation is as follows:

$$\sum_{i=1}^{4}q_{(i)\rightarrow(m,n)}+\dot{q}=0$$


We use 1 dymentional heat rate equations for convection and conduction as shown below:


$q_{(m-1,n)->(m,n)}=k\Delta_{x,y}\frac{T_{m-1,n}-T_{m,n}}{\Delta_{x,y}}$

$q_{(m+1,n)->(m,n)}=k\Delta_{x,y}\frac{T_{m+1,n}-T_{m,n}}{\Delta_{x,y}}$

$q_{(m,n-1)->(m,n)}=k\Delta_{x,y}\frac{T_{m,n-1}-T_{m,n}}{\Delta_{x,y}}$

$q_{(m,n+1)->(m,n)}=k\Delta_{x,y}\frac{T_{m,n+1}-T_{m,n}}{\Delta_{x,y}}$

Note that $k$ may be substituted with $h\Delta_{x,y}$ and $T$ with $T_\infty$ respectively.

In an iterative timestep finite-differencing system, Node temperatures are iterated over in a physical nodal matrix

$$T_{(m,n)}(s)=T_{(m,n)}(s-1)+\sum_{i=1}^{4}q_{(i)\rightarrow(m,n)}\Delta_{timestep}$$



Results

Case 1


Case 2


Case 3


Combined


Totals

$q_{0.01}=$
$q_{0.005}=$
$q_{0.001}=$
$q_{par}=$
868.8
$q_{id}=$
2400
$\eta_{0.01}=$
$\eta_{0.005}=$
$\eta_{0.001}=$
$\eta_{nom}=$
0.362831

Conclusion

In all numerical modeling methods, sacrefices are made in order to descritize a complex problem into tangible and computeable parts. One may only approach the theoretical values upon investing more computational effort. In this project, it is shown that while some mesh resolution may seem to be insufficent in computing a correct heat disipation value (e.g. 0.01 resoution is near rectangular fin shape), it is the binding mathematical principles that allow the solution to converge. What is important to take from this simulation is that it is indeed possible to solve the cases using Javascript and by object-oriented matrix interation. By iterating over the descrete nodes and using fundamental Energy Balance over time, convergance is possible in the time/timestep dymention, as in reality.

Apendix

Finite-Differencing Javascript Engine

function upperedge(x){ aval=t/2*(Math.pow((L-x)/L,2)); return aval; } function loweredge(x){ aval=-t/2*(Math.pow((L-x)/L,2)); return aval; } plot21=[];////internal nodes; plot22=[];//constant nodes; plot23=[];//convection nodes; plot31=[]; plot32=[]; plot33=[]; plot41=[]; plot42=[]; plot43=[]; endmatrix1=[]; endmatrix2=[]; endmatrix3=[]; function Runtest(steps){ stepx=steps; stepy=steps; xo=0; yo=-0.01; initialT=10; tplate=120; rows=(L)/stepx+2; cols=(t)/stepy+3; timestep=0.01; endtime=30;//10 seconds kval=20; h=100; //INITIALIZE MATRIX console.log("Matrix init"); var matrix=[];//rows for(var i=0;i < rows;i++){ matrix[i]=[]; for(var j=0;j < cols;j++){ var anode={}; anode['id']="("+i+","+j+")"; anode['T']=[]; anode['flux']=[]; anode['T'].push(initialT); anode['x']=i*stepx+xo; anode['y']=j*stepy+yo-stepy; anode['condition']='inactive'; anode['nid']=''; anode['sid']=''; anode['eid']=''; anode['wid']=''; matrix[i][j]=anode; }; }; console.log(matrix); //Prune Matrix for case console.log("Matrix map") for(var i=0;i < rows;i++){ for(var j=0;j < cols;j++){ //check nodes north south east west of node; //check east eexists=1; try{ if(typeof matrix[i-1][j] == 'undefined'){ } } catch(err){//area does not exist eexists=0; } if(eexists){ matrix[i][j]['wid']="("+(i-1)+","+j+")"; } //check west wexists=1; try{ if(typeof matrix[i+1][j] == 'undefined'){ } } catch(err){//area does not exist wexists=0; } if(wexists){ matrix[i][j]['eid']="("+(i+1)+","+j+")"; } //check north nexists=1; try{ if(typeof matrix[i][j+1] == 'undefined'){ } } catch(err){//area does not exist nexists=0; } if(nexists){ matrix[i][j]['nid']="("+i+","+(j+1)+")"; } //check south sexists=1; try{ if(typeof matrix[i][j-1] == 'undefined'){ } } catch(err){//area does not exist sexists=0; } if(sexists){ matrix[i][j]['sid']="("+i+","+(j-1)+")"; } } } console.log(matrix); //check matrix for internal and external nodes; console.log("calculating boundary conditions") for(var i=0;i < rows;i++){ for(var j=0;j < cols;j++){ currentx=i*stepx+xo; currenty=j*stepy+yo-stepy; if(currenty>=upperedge(currentx)+stepy-0.0000001){ //then the nodal point is outside of the fin if(currenty-stepy < upperedge(currentx)+stepy-0.0000001){ //then this point is a convection point matrix[i][j]['condition']='convection'; } else{ matrix[i][j]['condition']='inactive'; } } else if(currenty<=loweredge(currentx)-stepy+0.0000001){ //then the nodal point is outside of pin if(currenty+stepy>loweredge(currentx)-stepy+0.0000001){ matrix[i][j]['condition']='convection'; } else{ matrix[i][j]['condition']='inactive'; } } else{ //nodal point must be part of the fin; matrix[i][j]['condition']='internal'; matrix[i][j]['T'][0]=initialT; } //if at base: if(currentx==0){ //then it is part of the base of the fin and held constant matrix[i][j]['condition']='inactive'; matrix[i][j]['T'][0]=initialT; } if(currentx==0){ //then it is part of the base of the fin and held constant if(j>0&&j < cols-1){ matrix[i][j]['condition']='constant'; matrix[i][j]['T'][0]=tplate; } else{ matrix[i][j]['condition']='convection'; matrix[i][j]['T'][0]=initialT; } } if(currentx>0.1001){ if(Math.abs(currenty) < 0.001){ matrix[i][j]['condition']='convection'; matrix[i][j]['T'][0]=initialT; } else{ matrix[i][j]['condition']='inactive'; matrix[i][j]['T'][0]=initialT; } } } } console.log(matrix); console.log("configure plots"); for(var i=0;i < rows;i++){ for(var j=0;j < cols;j++){ currentx=i*stepx+xo; currenty=j*stepy+yo-stepy; if(matrix[i][j]['condition']=='internal'){ if(stepx==0.01){ plot21.push([currentx,currenty]); } if(stepx==0.005){ plot31.push([currentx,currenty]); } if(stepx==0.001){ plot41.push([currentx,currenty]); }; } if(matrix[i][j]['condition']=='constant'){ if(stepx==0.01){ plot22.push([currentx,currenty]); } if(stepx==0.005){ plot32.push([currentx,currenty]); } if(stepx==0.001){ plot42.push([currentx,currenty]); }; } if(matrix[i][j]['condition']=='convection'){ if(stepx==0.01){ plot23.push([currentx,currenty]); } if(stepx==0.005){ plot33.push([currentx,currenty]); } if(stepx==0.001){ plot43.push([currentx,currenty]); }; } } }; console.log("runtests") for(var k=0;k<(endtime/timestep)+1;k++){ //calculate heat flux for(var i=0;i < rows;i++){ for(var j=0;j < cols;j++){ if(matrix[i][j]['condition']=='internal'){ //calculate heat flux by neighbors var aflux=0; if(matrix[i][j]['wid']!=''){ if(matrix[i-1][j]['condition']=='constant'){ aflux=aflux+(kval*((matrix[i-1][j]['T'][k]-matrix[i][j]['T'][k]))); } if(matrix[i-1][j]['condition']=='internal'){ aflux=aflux+(kval*((matrix[i-1][j]['T'][k]-matrix[i][j]['T'][k]))); } } if(matrix[i][j]['eid']!=''){ if(matrix[i+1][j]['condition']=='internal'){ aflux=aflux+(kval*((matrix[i+1][j]['T'][k]-matrix[i][j]['T'][k]))); } if(matrix[i+1][j]['condition']=='convection'){ aflux=aflux+(h*stepx*(matrix[i+1][j]['T'][k]-matrix[i][j]['T'][k])); } } if(matrix[i][j]['nid']!=''){ if(matrix[i][j+1]['condition']=='internal'){ aflux=aflux+(kval*((matrix[i][j+1]['T'][k]-matrix[i][j]['T'][k]))); } if(matrix[i][j+1]['condition']=='convection'){ aflux=aflux+(h*stepx*(matrix[i][j+1]['T'][k]-matrix[i][j]['T'][k])); } } if(matrix[i][j]['sid']!=''){ if(matrix[i][j-1]['condition']=='internal'){ aflux=aflux+(kval*((matrix[i][j-1]['T'][k]-matrix[i][j]['T'][k]))); } if(matrix[i][j-1]['condition']=='convection'){ aflux=aflux+(h*stepx*(matrix[i][j-1]['T'][k]-matrix[i][j]['T'][k])); } } matrix[i][j]['flux'].push(aflux); } if(matrix[i][j]['condition']=='convection'){ //record heat loss from neighbors var aflux=0; if(matrix[i][j]['wid']!=''){ if(matrix[i-1][j]['condition']=='internal'){ aflux=aflux+(h*stepx*((matrix[i-1][j]['T'][k]-matrix[i][j]['T'][k]))); } } if(matrix[i][j]['eid']!=''){ if(matrix[i+1][j]['condition']=='internal'){ aflux=aflux+(h*stepx*((matrix[i+1][j]['T'][k]-matrix[i][j]['T'][k]))); } } if(matrix[i][j]['nid']!=''){ try{ if(matrix[i][j+1]['condition']=='internal'){ aflux=aflux+(h*stepx*((matrix[i][j+1]['T'][k]-matrix[i][j]['T'][k]))); } if(matrix[i][j+1]['condition']=='constant'){ aflux=aflux+(h*stepx*((matrix[i][j+1]['T'][k]-matrix[i][j]['T'][k]))); } } catch(err){}; } if(matrix[i][j]['sid']!=''){ try{ if(matrix[i][j-1]['condition']=='internal'){ aflux=aflux+(h*stepx*((matrix[i][j-1]['T'][k]-matrix[i][j]['T'][k]))); } if(matrix[i][j-1]['condition']=='constant'){ aflux=aflux+(h*stepx*((matrix[i][j-1]['T'][k]-matrix[i][j]['T'][k]))); } } catch(err){}; } matrix[i][j]['flux'].push(aflux); } if(matrix[i][j]['condition']=='constant'){ //record heat absorbtion by neighbors if(matrix[i][j]['eid']!=''){ if(matrix[i+1][j]['condition']=='internal'){ aflux=aflux+(kval*(stepx)*((matrix[i+1][j]['T'][k]-matrix[i][j]['T'][k])/stepy)); } } matrix[i][j]['flux'].push(aflux); } } }; //distribute heat flux as new T for next timestep; for(var i=0;i < rows;i++){ for(var j=0;j < cols;j++){ //if(Math.abs(matrix[i][j]['flux'][k]) < 1){ // matrix[i][j]['flux'][k]=0; //} //if(Math.abs(matrix[i][j]['flux'][k])>10000){ // matrix[i][j]['flux'][k]=10000*matrix[i][j]['flux'][k]/Math.abs(matrix[i][j]['flux'][k]); //} //if(matrix[i][j]['flux'][k]>0.1){ // matrix[i][j]['flux'][k]=0.1; //} if(matrix[i][j]['condition']=='internal'){ matrix[i][j]['T'][k+1]= matrix[i][j]['T'][k]+(matrix[i][j]['flux'][k])*timestep; } else{ matrix[i][j]['T'][k+1]= matrix[i][j]['T'][k]; } if(matrix[i][j]['T'][k+1]<0){ matrix[i][j]['T'][k+1]=0; } if(matrix[i][j]['T'][k+1]>200){ matrix[i][j]['T'][k+1]=120; } } } }; console.log(matrix); for(var i=0;i < rows;i++){ if(stepx==0.01){ endmatrix1[i]=[]; } if(stepx==0.005){ endmatrix2[i]=[]; } if(stepx==0.001){ endmatrix3[i]=[]; } }; for(var i=0;i < rows;i++){ for(var j=0;j < cols;j++){ if(stepx==0.01){ endmatrix1[i][j]=matrix[i][j]; } if(stepx==0.005){ endmatrix2[i][j]=matrix[i][j]; } if(stepx==0.001){ endmatrix3[i][j]=matrix[i][j]; } } }; } chart20data=[]; chart21data=[]; $(document).ready(function(){ plotchart1(); setTimeout(function(){ alert("Calculations will now start please wait ~30s for results."); Runtest(0.01); Runtest(0.005); Runtest(0.001); console.log("Endmatrix1"); console.log(endmatrix1); plotchart2(); plotchart7(); plotchart8(); plotchart9(); plotchart10(); plotchart11(); console.log("Endmatrix2"); console.log(endmatrix2); plotchart3(); plotchart6(); plotchart12(); plotchart13(); plotchart14(); plotchart15(); console.log("Endmatrix3"); console.log(endmatrix3); plotchart4(); plotchart5(); plotchart16(); plotchart17(); plotchart18(); plotchart19(); plotchart20(); plotchart21(); }, 10000); });