Tutorial for heat diffusion Simulator & Screenshots

EDNA 2013 (ENERGY DIFFUSSION & NUMERICAL APPLICATIONS) LINUX Release 1.0  is a fiendly user interface program which simulates the field temperature of a 3D solid body by mean of the heat diffusion equation solution. EDNA contains 4 modules:



The “out_XXX.dat” file is a file used to display 3D surfaces in a visualizing software as TECPLOT 360 by default but it can be changed soon.

Download from google code site:

First of all open a linux terminal
Captura de pantalla de 2013-07-05 16:20:12

Go to the directory where you downloaded EDNA.tar

Captura de pantalla de 2013-07-05 16:21:10

Extract EDNA.tar:

Captura de pantalla de 2013-07-05 16:21:55

Outuput screen seems like:

Captura de pantalla de 2013-07-05 16:22:06

Go to the new directory that it was already created “EDNA”

Captura de pantalla de 2013-07-05 16:22:22
Compile with gcc 

Captura de pantalla de 2013-07-05 16:22:55
Output screen from make

Captura de pantalla de 2013-07-05 16:23:03

Clean .o and .mod files

Captura de pantalla de 2013-07-05 16:23:25

output screen


Captura de pantalla de 2013-07-05 16:23:31

run program. Next time you want to run EDNA just start from here and forget steps have been explained before

Captura de pantalla de 2013-07-05 16:23:53

Menu for EDNA. Chose 1 for running

Captura de pantalla de 2013-07-05 16:23:59

like this

Captura de pantalla de 2013-07-05 16:24:08

Fist module: Dimensions and grid. You can change dimensions for the cube and number of points. Be careful for the number of points because if there are too many points, the computational demand would be very high.

Captura de pantalla de 2013-07-05 16:24:12

Second Module: You can change material of the body by mean of a list materials or you can insert your own material information according to the diffusion coefficient. Anysotropic materials are valid.

Captura de pantalla de 2013-07-05 16:24:34

Boundary conditions: You can fix any wall temperature according to the coordinated system.

Captura de pantalla de 2013-07-05 16:24:41

Time options are available and you can change the time in second you want to finish. Also you can provide “timing save” which is the time when a file is going to be written. Be careful with this number

Captura de pantalla de 2013-07-05 16:24:46

You need to agree with the number of files in order to compute

Captura de pantalla de 2013-07-05 16:24:57

Files are written

Captura de pantalla de 2013-07-05 16:25:20

Finish with zero

Captura de pantalla de 2013-07-05 16:25:44

Open Tecplot. This sofware is required in order to see the results from out_XXX.dat files

Captura de pantalla de 2013-07-05 16:25:56

Go to file menu->Load data

Captura de pantalla de 2013-07-05 16:28:08

Chose Tecplot data loader->ok

Captura de pantalla de 2013-07-05 16:28:24

okCaptura de pantalla de 2013-07-05 16:28:30

Select multiple filesCaptura de pantalla de 2013-07-05 16:28:49 Select all files from EDNA directory at the same timeCaptura de pantalla de 2013-07-05 16:29:08 click on “Add”Captura de pantalla de 2013-07-05 16:29:15 like thisCaptura de pantalla de 2013-07-05 16:29:18

Files are loaded and take little time depending of the number of files and number of points chose 3D Cartesian
Captura de pantalla de 2013-07-05 16:31:16
and tick “Show first zone only”

“>Now Tick Slices
Captura de pantalla de 2013-07-05 16:31:27 click on options buttonCaptura de pantalla de 2013-07-05 16:31:39

click on “2”

Captura de pantalla de 2013-07-05 16:31:51

tick “Show group 2” in order to see slices in “y” axis

Captura de pantalla de 2013-07-05 16:31:58 Captura de pantalla de 2013-07-05 16:32:12 Captura de pantalla de 2013-07-05 16:32:16 Captura de pantalla de 2013-07-05 16:32:22

close the options silce window and go to menu toolbar “Animate”

Captura de pantalla de 2013-07-05 16:32:36 Chose zone animation and a window will apear like this

Captura de pantalla de 2013-07-05 16:32:51


Chose the last file you got from EDNA and click on “Animate”. The field temperature will be displayed on the screen. You can optionaly save your simulation on a avi movie file in high definition (resolution:1000 )


Captura de pantalla de 2013-07-05 16:33:11 Captura de pantalla de 2013-07-05 16:33:16 Captura de pantalla de 2013-07-05 16:33:27Support:


Numerical Simulation of Supersonic Flow Over OREX Reentry Capsule

The Orbital Re-Entry Flight Experiment (OREX) was the Japanise experiment implemented in 1994 to investigate flight conditions of capsules. A better performance of thermal protection systems and aerodynamic heating were the goals. Experiments were similar to those used by NASA.

A high-speed flow over reentry capsule generates a bow shock wave. This shock wave causes high pressure and temperature over the surface, generating damages for the reentry vehicle.


Figure 1

In front of the body, (zone a) figure 1, the bow shock is detached and between them there is a supersonic zone flow and subsonic flow. Each zone is divided by the sonic line (zone b) whose position is of important interest for reserachers since it defines the magnitude of  global drag coefficient.

When the flow comes at the edge of the vehicle,  it turns and expands very fast (zone c), and the boundary layer detaches from the body. This phenomenon creates a free shear layer between two different velocities (zone e) thus the free shear layer divides the recirculating region behind reentry capsule (zone d) and the free stream. Eventually, the flow is recompressed downstream by the tail shock (zone f), which is formed because the flow becomes supersonic when it passes through the expansion fan on the corner.

The Numerical simulation was carried out using sixth order schemes to solve derivatives of Navier-Stokes equations. Mach number is 1.5 and the fluid is air. Boundary condions were implemented according to [1]. LES technique was used to solve turbulent structures.

Thanks to PhD. Martin Salinas of the Engineering Intitute of UNAM.


[1] Poinsot and Lele “Boundary Conditions for Direct Simulations of Compressible Viscous Flow “. Stanford 1991

More info related with OREX http://www.rocket.jaxa.jp/fstrc/0c01.html

How to Solve the Heat Equation?

Transient problems in heat transfer are restricted to simple geometries and simple analytical solutions. Nevertheless, in many cases the geometry and boundary conditions make difficult the use of analytical techniques, therefore we can make use of the finite-difference methods.

Numerical method for the heat equation

For the two-dimensional system of Figure 1 under transient conditions with constant properties, the heat equation is:



Figure 1. Discretized computational domain with finite-differences.

We may use the central-difference approximations to the spatial derivates. The subscripts (i, j) are used to designate the locations of the  nodal points,in x and y respectively. The heat equation equation must be discretized in space and in time to obtain a numerical solution which depends on x, y and t. To achieve this purpose, the integer p is used to calculate the time:


Finite-difference approximation to the time derivate is represented as:


p denotes the time dependence of T,  the time derivative is expressed in terms of the difference in temperatures associated with the new temporal step (p+1)  and previous temporal step (p). Finite-difference approximation to the spatial derivate is represented as:

figura 2

The explicit form of the finite-difference equation for the interior node i, j is:


The last expresion is found just subtituting the explicit form of the spatial derivative in the heat equation. Then, solving for the nodal temperature T(i,j) at the new temporal step (p+1) and assuming that Δx= Δy, it follows that:


The last equation is explicit because unknown nodal temperatures for the new time are determined by known nodal temperatures at the previous time, and F0  is the Fourier number defined as:


The temperature of each interior node is known at t=0 (p=0) from prescribed initial conditions, the calculations begin at t= Δt , (p=1). With temperatures known for t=Δt, the appropriate finite-difference equation is the applied at each node to determinate its temperature at t=2 Δt (p=2). The transient temperature distribution is obtained by marching out in time, using intervals of Δt.

In order to prevent errors due to instability of the method, the prescribed value of Δt must be maintained below a certain limit, which depends on Δx and other parameters of the system. The stability criterion is determined by requiring that the coefficient associated with the node of interest at the previous time is greater than or equal to zero. The stability criterion depends on Fourier number, and for a two dimensional interior node is:



An easy example is a square plate. At first we are going to consider a material with Fourier number of 0.2, and boundary conditions with constant temperature. Temperature on the top is 20°C, and the others sides are keep constant at 100°C. The temperature at the beginning of the time is also 100°C.

The computational domain is shown below and, it represents those points in which the temperature can be calculated.


According to the domain, we are going to implement the heat equation in each node, and results are shown in the next table. We must notice that nodes (1,1) (2,1) (3,1) (4,1) (1,2) (1,3) (1,4) (4,1) (4,2) (4,3) and (4,4) can not be calculated because they are constant boundary conditions.


The result for three dimensional  heat equation with p=1500 steps is shown below.

Download EDNA

Download f90

EDNA is a powerful program which simulates the field temperature of a 3D solid body by mean of the heat diffusion equation solution.


“Fundamentals of Heat and Mass Transfer”, Incropera, sixth edition.

Peristaltic Flow of a Newtonian Fluid

The peristaltic flows of Newtonian fluids in complex geometries are interesting to medical and industrial applications. In human body this plays a important role in the transport of food, urine, blood and egg to the uterus.

In the present investigation was analysed the 3D peristaltic flow in a complex channel with progressive waves of wall contraction. The peristaltic flow is investigated in the wave frame with c relative velocity with the frame of laboratory. This peristaltic flow was investigated varying wave amplitude, wave velocity, wave length and phase difference.


To describe this complex geometry BFC (Body Fitted Coordinates) method has been employed in this study.


Boundary conditions were imposed over walls to simulate peristaltic wave movement, and this peristaltic wave had a constant speed.


Amplitude wave velocity  results.


Length velocity  results.


Phase difference velocity results.


Independence of velocity


Numerical Simulation of Taylor Vortices

The annular flow between concentric cylinders is known as Taylor-Couette flow, and when dynamic equilibrium between the radial forces and th
e radial pressure disappears, a secondary flow results from this imbalance. In this secondary flow, some  toroidal structures are formed along axial direction.

In this analysis, 3D Taylor-Couette flow has been studied solving Navier-Stokes equations using finite volume technique and a commercial CFD software Package (PHOENICS).

A cylindrical mesh with 440,000 tetrahedral elements was used for system description.


Boundary conditions play an important role for formation of Taylor-Couette flow and rotating boundary conditions were imposed on inner cylinder, the highest speed was found on the surface of the inner cylinder and a speed equal to zero on the outer cylinder.perfil vel_1

The toroidal structures appear in couple along the axial direction and spin in the same axis, but also these structures are rotating among themselves, as it is shown in the next figures.


Taylor number was varied to study the appearance of this structure known as Taylor vortex, and how it changes the vortex shape.

taylor vortex_1

Velocity profile has a sinusoidal behavior for all Taylor numbers employed in this study as it is shown below.


Numerical Simulation of Supersonic Jet

Turbulent under-expanded jet is an important fundamental flow, involving interactions between supersonic shock-wave and turbulent compressible flow. Jets have practical applications, like exhaust plumes of propulsion vehicles, investigations on supersonic combustion, reduction of noise emissions like screech tones which are a particular interest because acoustic equipment and materials properties can be affected and gaseous fuel injection system. Bearing in main this considerations, a under-expanded supersonic jet of air was simulated during the present work using LES method. The main reason to do this work is because measurements of the structure of under-expanded jets are limited, due to difficulties caused for variation of static pressures, velocities, and shock waves at the time of introducing equipment.

There is a complicated structure inside of the jet which has some remarkable features. First, the jet boundary oscillates as the supersonic jet periodically expands and converges in zones where expansion fans and shock-waves take place. These shocks waves (red color) alternate with expansions fans (blue color). The supersonic jet expands and cools as it flows through the expansion fans and is compressed and heats as it flows through the shock wave. If the supersonic jet is circular, then the axisymmetric shock wave and fans have conical shape.

Shock-wave cells of Supersonic Jet


Number of nodes nx=155, ny=109, nz=109

Solution for Cartesian Navier-Stokes equations.

Boundary condition for compressible Navier-Stokes equations [1]

Sixth-order finite differences scheme (spatial discretization).

Third-order Runge-Kutta method (time discretization)

Wave radation Ma=1.4 of supersonic Jet


Transient velocity vectors indicating the pressure field.

Isosurfaces of criterion Q which represent physically speaking, the development of turbulence.


[1] Poinsot and Lele “Boundary Conditions for Direct Simulations of Compressible Viscous Flow “. Stanford 1991.

Periodic Boundary Conditions


Real Boundary conditions for Navier-Stokes are extremely complex. It is possible to find solutions for Navier-Stokes equations introducing periodic boundary conditions. Flows simulated appear in a computational domain which is moving at a constant velocity.

Periodic boundary conditions for Navier-Stokes are applied using the same value of the rate of physical properties of the flow in the first and the last node of the computational domain. Some high-order schemes [1] uses more nodes, and it is necessary apply the same procedure, for example the second node has the same value of the penultimate node, and so on. Periodic boundary conditions have the advantage of short time of computation.

We have simulated a compressible jet which has periodic boundary conditions in II-UNAM.

  • Fluid: Air
  • Periodic Boundary Conditions
  • Free code
  • By Abraham Aguilar II_UNAM
[1] “On the Use of Higher-Order Finite-Difference Schemes on Curvilinear and Deforming Meshes”,  Miguel R.Visbal and Datta V. Gaitonde 2002.