This is an update to my project work. The previous post can be found here, that provides some speed-up stats in matrix multiplication using CUDA.

The next problem to tackle before targeting Fluid Flow problems would be related to heat transfer. As we know, heat conduction is governed by a single partial differential equation:

$\frac{\partial T}{\partial t} = \alpha\nabla^2T.$

Where, T is Temperature; t is Time; $\alpha$ is Thermal Diffusivity

And for a 2-Dimensional plate, this equation would look like:

$\frac{\partial T}{\partial t} = \alpha(\frac{\partial^2 T}{\partial^2 x} + \frac{\partial^2 T}{\partial^2 y}).$

In order to solve this equation numerically, you need a grid. I chose a square plate and divided it into a grid consisting of 512 or 1024 nodes. It would look something similar to this:

And since we intend to solve the above PDE numerically, we must find out a way to represent it by a Finite Difference Scheme. A possible approximation could be:

$\frac{\Delta T}{\Delta t} = \alpha(\frac{T_{i+1,j} - 2T_{i,j} + T_{i-1,j}}{\Delta x^2} + \frac{T_{i,j+1} - 2T_{i,j} + T_{i,j-1}}{\Delta y^2}).$

Since we choose a square plate and a square grid, we divide the $x$ and $y$ directions in the same number of small divisions, i.e. $\Delta x = \Delta y.$

In order to convert some of the parameters to non-dimensional form, we rewrite the equation in the form:

$\Delta T = \frac{\alpha.\Delta t}{\Delta x^2}(T_{i+1,j} - 2T_{i,j} + T_{i-1,j} + T_{i,j+1} - 2T_{i,j} + T_{i,j-1}).$

and we denote the expression $\frac{\alpha.\Delta t}{\Delta x^2}$   as   $\alpha_{factor}.$

I looked up in the bible for Heat Transfer by Prof. Frank Incropera and found a numerical example as a reference to my problem. According to the literature, a decent assumption would be:   $\alpha_{factor} < 0.25.$

So my inputs to the problem were:

1. 512×512 square grid.
2. Initialized temperature for the 4 corners of the plate, i.e. indirectly setting up the boundary edge temperatures.
3. Time steps till approximate steady-state arrival $t_{steps} = 30000.$
4. $\alpha_{factor} = 0.2.$

We can now iterate over each grid node, updating it in each time step and running this operation over a period of time defined by 30000 time steps.

This can be done easily by a C code:

void cpu_compute(int ni, int nj, float alpha, float *temp_in, float *temp_out){

int i, j, i00, im10, ip10, i0m1, i0p1;
float d2tdx2, dt2dy2;

for(j=1; j<nj-1; j++){
for(i=1; i<ni-1; i++){
i00 = I2D(ni, i, j);
im10 = I2D(ni, i-1, j);
ip10 = I2D(ni, i+1, j);
i0m1 = I2D(ni, i, j-1);
i0p1 = I2D(ni, i, j+1);

d2tdx2 = temp_in[im10] - 2*temp_in[i00] + temp_in[ip10];
d2tdy2 = temp_in[i0m1] - 2*temp_in[i00] + temp_in[i0p1];

temp_out[i00] = temp_in[i00] + alpha*(d2tdx2 + d2tdy2);
}
}
}


We use the Gauss-Seidel iterative method here for updating the node values. After the function is run, we swap the input and output temperature arrays. This continues for the initialized number of time steps.

The GPU version of this can be found here in another post. I have primarily explained the functioning of the GPU kernel.

The initial temperatures at the 4 corners A, B, C and D were taken as 200, 200, 300 and 300 K respectively. Since the temperatures are only added or subtracted, they need not compulsorily be in Kelvin scale, we may even use the Celsius scale. After the corner temperatures are set, we calculate the border or edge temperatures. We could argue that the steady state temperature at the edge AB would be constant at 200 K. Similarly edge CD would have a constant temperature of 300 K. Now, edge AD and BC could be approximated to have temperatures varying linearly with distance (from T = 200 K to T = 300 K). The Gauss-Seidel iterative algorithm is used to find the temperature of the inner grid nodes, i.e. the 510×510 center square.

Satisfactory results were obtained after 30000 iterations. We may export the obtained data as a text file or write it to a data file (.dat). I used MATLAB to plot 2-D and 3-D contours for these outputs. The following contours were obtained.

A detailed explanation for the dimensional analysis can be found here. I am currently trying to simulate the same problem on ANSYS FLUENT and compare results. Any new developments would be added here.

From → GPU & GPGPU, Project

Hí! Nice post!!
Can you help me?
If I have all 3d steady state heat equations developed for a lot of boundary conditions (interiror node, external node with convection, face node with heat flux etc), finite element method, how can I solve my mxmxm thermal matrix, using Gauss Seidel & Matlab?