Numerical Method


The diffusion simulation is base on the rules :

 1. The temperature changes follow the formula of diffusion
            Ut = D * Uxx
 2. The boundary condition is isolated
            Uxn+1 = Uxn

Here is the core part of this program.   The formula of diffusion will be
convert into program code.  

The Formula of Diffusion : Ut = D * Uxx

U is a function that record the temperature values alone x-axis.    Meanwhile , it
changes by time t.
That means :
        U = U(x,t)

Ux means the slope of U at x point.   Or ,
how temperature changes when it moves alone x-axis at x point.

Let's consider about temperature value function U at point x , U(x).
Assume h is very small distance alone x-axis , or h -> 0
        2h  =  U(x+1)  -  U(x)  =  U(x)  -  U(x-1)
        U(x+2h) = U(x+1)
Then by definition
        Ux(x) = [ U(x+h) - U(x-h) ] / 2h

    
    Uxx(x) = [ Ux(x+h) - Ux(x-h) ] / 2h

        Ux(x+h) = [ U(x+h+h) - U(x+h-h) ] / 2h = [ U(x+1) - U(x) ] / 2h
        Ux(x-h) = [ U(x-h+h) - U(x-h-h) ] / 2h = [ U(x) - U(x-1) ] / 2h

 Then we have ,
  

        Uxx(x) = [ Ux(x+h) - Ux(x-h) ] / 2h = [ U(x+1) - 2U(x) + U(x-1) ] / (4h^2)


Apply this to y-axis

        Uyy(y) = [ U(y+1) - 2U(y) + U(y-1) ] / (4h^2)

Apply this to Ut (For 2D , dt -> 0 ),

        Ut = D * ( Uxx + Uyy) =  D *
{ [ U(x+1) - 2U(x) + U(x-1) ] / (4h^2) + [ U(y+1) - 2U(y) + U(y-1) ] / (4h^2) }
             = D * [ U(x+1,y) + U(x-1,y) + U(x,y+1) + U(x,y-1) - 4U(x,y) ] / (2h^2)
 
      U(t+1) = U(t) + dt * D1 / h^2
* [ U(x+1,y) + U(x-1,y) + U(x,y+1) + U(x,y-1) - 4U(x,y) ]


This is numerical form of Diffusion Equation ( Ut = D * Uxx ) !
Each calculation will get the temperature value after time interval dt.

This method is called Finite Difference Method.

[The Meaning of dt * D1 / h^2 ]

Assume c =
dt * D1 / h^2
                   = dt * D2 / h1^2
                   = dt * Dcof / dx^2


Here h1 = 2h = U(x+1) - U(x) = dx = The width of each calculation nod
        D2 = Dcof = The diffusion coefficient of a material

The unit of :
        U(x) : Celcius
        dt : sec
        dx : cm
        Dcof : cm^2 / sec

The dt , dx should be very small.    Then we can set
        dt = 0.01 sec
        dx = 0.01 cm

So , we can use 500 nods to present a material that has 5cm width.
And each calculation will get the result of the time passed by 0.01 sec

In this program , the c = 0.200 , it is a constant in here.

        
c =
dt * Dcof / dx^2
        Dcof = 0.002 (cm^2/sec)

This says :
IF , we set dt = 0.01 sec , dx = 0.01 cm
Then , the program shows the simulation result of a material that has
Dcof = 0.002
(cm^2/sec)

[Remark]

You'll get wrong result if writing code like this :

    
    U(x) = U(x) + c * [ U(x-1) - 2U(x) + U(x+1) ]

This is the correct form :

    
     U2(x) = U(x) + c * [ U(x-1) - 2U(x) + U(x+1) ]


And , the value of c should not exceed 0.5.



[Hint About Programming]

class Material {
  bool bActive;
  float value;
} u[n];

In this 1D case , the user can pick up any point from 1 to (n-1) for diffusion.
Say , n == 100 , but the material only occupy u[5] to u[10]
So ,
// initialize
for(int i = 5 ; i<=10 ; i++ ) {
  u[i].bActive = true;
  u[i].value = ini_value ;
}

// for calculation
for(int i = 1 ; i<=(n-1) ; i++) {
  if( u[i].bActive == true ) {
     // deal with boundary condition
     if( u[i-1].bActive == false ) u[i-1].value = u[i].value;
     if( u[i+1].bActive == false ) u[i+1].value = u[i].value;
     normal calculation
  }
}

// draw on screen
show( u[i] );

The example on the web page only calculate 50 times per sec ,
update show( ) once after more calculations can make it much faster.





Copyright © 2010 bdragon All rights reserved.
bdragong@gmail.com

WebRep
Overall rating
 
Comments