Multiphysics Modeling Using
Cellular Automata

Brian Vick

Mechanical Engineering Dept
Virginia Tech
Blacksburg, VA
bvick@vt.edu

Objective

To propose a new method for modeling and solving multi-physics phenomena involving mechanical, thermal, chemical, and electrical processes.

Basic Idea

Combination of fundamental ideas from traditional finite differencing, analytical analysis, and cellular automata.

A computationally explicit yet numerically stable set of rules is developed for elementary processes.

Elementary processes are assembled in a cascading sequence to form complete, complex processes.

Fundamental Idea

U=collection of state variables

General Process

These state variables change with time due a total process consisting of N elementary processes

((∂U)/(∂t))=(∑_(n = 1)^N)(Subscript[P, n])[U]
	(Subscript[P, n])  =  an elementary process

Elementary Processes

If the state of the system evolved according to only one of the elementary processes, the rate equation would be

((∂U)/(∂t))=(Subscript[P, n])[U]

This process is applied over a discrete network of cells and is advanced in time in discrete steps.  The discrete rule is

(U^(k + 1))=(Subscript[ϕ, n])[(U^k)]

Construction of Complex Processes

Once rules for elementary processes have been developed, a complex process consisting of N elementary processes is constructed as a cascade of the elementary processes.

(U^(k + 1))=(Subscript[ϕ, N])[...(Subscript[ϕ, 2])[(Subscript[ϕ, 1])[(U^k)]]...]

Transport Processes

General Transport Processes

The dynamics is governed by transport due to the following fundamental mechanisms: 
1) Advection due to bulk motion of the medium (Convection)
2) Diffusion 
3) Reactions kinetics (chemical or biological processes)
4) Interactions with an external environment

((∂U)/(∂t))=-V∇U+D(∇^2)U+R[U]+S[r, t, U]

Elementary Transport Processes

External Interactions

((∂U)/(∂t))=S[r, t, U]

Advection

((∂U)/(∂t))=-V∇U

Diffusion

((∂U)/(∂t))=D(∇^2)U

Reaction Kinetics

((∂U)/(∂t))=R[U]

Rules for the Elementary Processes

Criteria for the Rules:
  1) Physically realistic
  2) Computationally explicit
  3) Numerically stable

External Interactions

((∂U)/(∂t))=S[r, t, U] = -H U  +  g

The rule for advancing the state variable over a time increment Δt can be immediately deduced from the solution of this linear ordinary differential equation.

(U^(k + 1))=(Subscript[ϕ, ext])[(U^k)]
	(Subscript[ϕ, ext])[U]=U*exp[- H*Δt]+(g/H)(1-exp[-H*Δt ])

Rules for the Elementary Processes

Advection

For motion in only one direction

((∂U)/(∂t))=-V((∂U)/(∂x))

Since this differential equation has a traveling wave solution, a simple and physically correct rule for advancing the state variable is

(U^(k + 1))=(Subscript[ϕ, adv])[(U^k)]
	(Subscript[ϕ, adv])[U]=Rotate[U,V*Δt/Δx]

Rules for the Elementary Processes

Reaction Kinetics

Usually this term is of a nonlinear nature and as a result a convenient rule based on the analytical solution is impossible to find.  However, a first-order Taylor series expansion about the current state variable provides a suitable accurate rule

((∂U)/(∂t))=R[U]=R[(U^k)]+(((R'))^k)(U-(U^k))

This linearized equation provides for a rule for reaction kinetics as

(U^(k + 1))=(Subscript[ϕ, react])[(U^k)]
	(Subscript[ϕ, react])[U]=U*exp[ R'*Δt]-((R[U] - R'*U)/(R'))(1-exp[R'*Δt ])

An example is the logistic growth model used in population dynamics

R[U]=A*U (1 - U/B)

Rules for the Elementary Processes

Diffusion

2 Cell Interaction

Using a finite difference type approach, we get

{(((Subscript[u, i]))^(k + 1)),(((Subscript[u, i + 1]))^(k + 1))}=(Subscript[ϕ, diff, cells])[{(((Subscript[u, i]))^k),(((Subscript[u, i]))^k)}]
	(Subscript[ϕ, diff, cells])[{(Subscript[u, i]), (Subscript[u, i + 1])}]={(Subscript[u, i]) + Δu, 
              (Subscript[u, i + 1]) - Δu}
	Δu=(1/2)((Subscript[u, i + 1]) - (Subscript[u, i]))(1-exp[-2(D/(Δ(x^2)))*Δt ])

Rules for the Elementary Processes

Multi-dimensional Systems

All these rules are easily extended to multi-dimensional grids of cells and irregular networks of cells .

Rules for Elementary Processes: Mathematica Code

External Interaction

Clear[ϕext, ϕadv, ϕdiffCells, ϕdiff, ϕreact]

ϕext[u_List] := 
      u*Exp[-h*Δt] + 
        g/h*(1 - Exp[-h*Δt])
ϕext[u_List] := u + g*Δt /; h == 0

Advection

ϕadv[u_List] := 
      RotateRight[
        u, ((V*Δt)/Δ)x]
ϕadv[u_List] := u /; V == 0

Diffusion

ϕdiffCells[{u1_, u2_}] := {u1, 
          u2} + {1, -1}*(1/2) (u2 - u1) (1.  - 
              Exp[-2 (α/Δ)(x^2) Δt])

ϕdiff[
        u_List] := 
Module[{i, 
          unew}, 
unew = 
          u; 
Do[
unew[[{i, 
                  i + 1}]] = ϕdiffCells[
              unew[[{i, i + 1}]]], 
{i, 1, Ncells - 1, 
            2}
]; 
Do[
unew[[{i, 
                  i + 1}]] = ϕdiffCells[
              unew[[{i, i + 1}]]], 
{i, 2, Ncells - 1, 
            2}
]; 
unew[[{Ncells, 
                1}]] = ϕdiffCells[
            unew[[{Ncells, 
                  1}]]]; 
unew
]

Reaction Kinetics

ϕreact[u_List] := 
      u*Exp[Rp[u]*Δt] - (((R[u] - Rp[u]*u)/Rp)[u])*(1 - 
              Exp[Rp[u]*Δt])

R[u_] := a*u*(1 - u/b);
Rp[u_] := a*(1 - 2*u/b);

Initial Conditions

RandomList[n_] := Table[Random[Integer], {n}]

RandomList[n_, k_] := 
      Table[Random[Integer, {0, k - 1}], {n}]

RandomListReal[n_, x_] := Table[Random[Real, {0, x}], {n}]

CenterList[n_Integer, a_List] := 
      With[{len = ((n-Length[a] )/2)}, 
        Join[Table[0, {⌊len⌋}], a, 
          Table[0, {⌈len⌉}]]]

Case Studies: Elementary Processes

Diffusion

Square Wave Initial Distribution

Random Initial Distribution

Effect of Diffusivity: Square Wave Initial Distribution

\(\(D\^+\) = \(dimensionless\ diffusivity = \(D*Δt\)\/Δx\^2\)\)

Case Studies: Collection of Elementary Processes

Square Wave Initial Distribution

Random Initial Distribution

Case Studies: Complex Processes

Advection/Diffusion/External Interactions

Effect of Velocity and External Interaction (D=0.25)

Case Studies: Complex Processes

Reaction Kinetics/Diffusion/Advection

Accuracy Studies

Summary and Conclusions

> A new methodology, based on a combination of differencing, analytical analysis, and cellular automata concepts, has been developed.  It involves a cascading sequence of simple, explicit rules of evaluation, rather than complicated partial differential equations.   The resulting scheme is computationally explicit yet numerically stable.  

> The power and flexibility of the method is demonstrated with a general transport process consisting of elementary processes due to diffusion, advection, reaction kinetics, and external interactions.  Numerical experiments show that the method is highly accurate if the time step is sufficiently small.  
 
> The value of the method is that a library of rules for elementary processes can be developed and these processes can be combined to form any application of interest. The method offers a relatively simple means to model complex, multi-physics processes involving coupled mechanical, thermal, chemical, and electrical effects.

> The proposed methodology frees the user from the shackles of differential equations.