Computational methods

In this page, you can download some MATLAB and FreeFem scripts I wrote during my studies at Politecnico di Torino.

To download the programs, please click on the highlighted links: they will lead you to a new page in Google Drive, where you will have to click the “Download” button.

All the programs in this page are free of charge: they are under a Creative Commons licence and cannot be used for profit (CC-BY-NC).


This algorithm allows the calculation of the pinch point temperature in a chosen network of heat exchangers. Moreover, it computes the minimum heat needed by the network. The pinch point temperature difference is chosen by the user in the MATLAB Workspace directly and the Grand Composite Curve is computed and printed. The network of heat exchangers is defined in a txt file of this kind:
1 5 200 100
2 11 170 55
3 14 120 165
4 8 50 150
The first column indicates the fluid number, the second column is the product between the specific heat capacity and the mass flow (in kW/K), and the third and fourth columns represent the inlet and outlet temperatures, respectively. The order of the fluids does not influence the result and the txt file needs to be put in the folder where the script is.
I also wrote a JAVA PROGRAM (platform-independent .jar file) for Pinch Point analysis: it can be downloaded HERE.


Below this point you will find some exercises that I solved.

TRANSIENT – 1-D heat transfer (Download)

A transport current of 50 A flows in an electrical conductor (laterally insulated, with length L = 1 m, resistivity rho = 1.75e-8 Ωm; diffusivity α = 11e-5 m2/s, conductivity k=350 W/m/K), with diameter D = 4mm. Compute the temperature evolution in the middle of the conductor, which is initially at  T(x,0)=T_0= 300K , with the boundaries kept at the given constant temperature T_0.

STEADY STATE- 1-D heat conduction with radial geometry (Download)

Pressurized water at ~ 400 K flows in a cylindrical pipe (inner diameter D_in = 20 mm, outer diameter D_out = 24 mm, conductivity k = 0.035 W/m/K, heat transfer coefficient h_in = 100 W/m2/K). On the outside the pipe is cooled by air at 300 K (heat transfer coefficient h_out = 10 W/m2/K). Compute numerically the temperature profile in the pipe wall cross section, using the FD method.

TRANSIENT – 1-D heat conduction with radial geometry (Download)

A transport current of 10 A flows in a hollow cylinder (inner diameter D_in = 2 mm, outer diameter D_out = 4 mm, resistivity ρ = 2.5e-7 Ωm; diffusivity α = 11e-5 m2/s, conductivity k = 3.5 W/m/K), having a length much larger than the transverse dimension. The conductor is insulated on the outer side, while the inner side contains a fluid at a constant temperature T_0 = 300 K, which is also the initial temperature of the cylinder (heat transfer coefficient with fluid h = 100 W/m2/K). Solve numerically the transient problem and compute the steady-state temperature profile on the conductor cross section, using the FD method.

TRANSIENT – 1-D heat conduction (Download)

A superconducting strand (diameter D = 0.81 mm, length L = 1 m, diffusivity α = 11e-5 m2/s, conductivity k = 350 W/m/K), initially at 4.5 K, carries a current of 50 A. The electric field generated by the current while the strand is superconducting is E = E_0 (I/I_C)^n, where E_0 = 1e-5 V/m and n = 1.2. I_C is the critical current and it decreases with temperature: I_C = 50 – T_2. The boundaries are kept at 4.5 K. Compute numerically the maximum temperature increase in the middle of the strand.

ADVECTION with constant speed (Download)

A fluid (ρ = 130 kg/m3, c_p = 5500 J/kg/K) flows in a pipe (diameter D = 4 mm, length L = 2 m). Consider the 1D flow along the conduit length, in view of the fact that L>>D. Write a program which uses the upwind finite difference scheme to compute the numerical solution of the 1D pure advection problem: dT/dt+c*dT/dx=0, with c = 0.005 m/s, and: T(x,0)=4.5K; T(0, t=0) =15K.
Note that we are specifying boundary values for the solution on the left boundary x = 0, but not on the right boundary x = 2 m: why? Plot the temperature profile along the pipe at t = 50, 100, 150 and 200 s.

ADVECTION: Burger’s Equation with shock (Download)

The temperature profile in a pipe (-10<=x<=10) evolves according to the equation:
dT/dt+T*dT/dx=0 with T(x,0)=f(x), and:
f(x)= 5 (x<-5); 2*cos(pi*(x-5)/2)+3 (-5<=x-3); 1 (x>-3)
Using MATLAB, compute and plot the characteristic lines in the (x,t) piane, exiting from x0 = -6, -5, -4,-3, … , 3, 4, 5. Compute analytically and numerically the temperature evolution at times t = 0, 1, 2, 3, 4, 5 in the portion of the domain [-5, 10]. Compare the analytical and numerical solutions in a single plot.


The Finite Elements method is studied with the FREE software FreeFem++. On the FreeFem++ website, you will find the official documentation for the program. The solution of some exercises in Italian language is available here.

FreeFem++ example (Download)

A hollow cylinder (r_in=0.02 m, r_out=0.1 m) has a starting temperature of 25 degrees C. It has a density of 1500 kg/m3, a thermal conductivity of 1 W/mK, and a specific heat capacity of 500 J/kgK. A fluid with a temperature of 30 degrees °C and a convective heat transfer coefficient of 500 W/m2K flows around the cylinder. Compute the temperature evolution of a section of the cylinder with a tolerance of 10E-6.