Simulation of Dynamic Systems with Engineering Applications Using Mathematica

Dr. John (Kohoutek^( *†))

* Nanofabrication Research Group, Center for Nanoscale Science and Technology, National Institute of Standards and Technology, Gaithersburg, MD
† Institute for Research in Electronics and Applied Physics, University of Maryland, College Park, MD

Abstract

A finite difference method has been developed from scratch using Mathematica.  Here, this method is compared to both a commercial software package and a theoretical calculation; convergence testing has also been performed.  Computational differences in the scratch method are noted, which may explain the deviation between the commercial software package result and the theoretical result for a given problem.  Which computational method is correct?  The audience can decide.

A few words

-This method was developed from scratch with minimal help from the literature or other sources.

Finite Difference Method

The finite difference method uses finite differences to approximate derivatives in different directions.  In one dimension, a first derivative looks like:

f'(x)=((f(x2)-f(x1))/(Δx))

A second derivative looks like:

f''(x)=(((f(x3)-f(x2))-(f(x2)-f(x1)))/((Δx)^2))=((f(x1)-2 f(x2)+f(x3))/((Δx)^2))

In the time domain, you might have something like:

(df/dt)=(((d^2) f)/(dx^2))

(Subscript[f, new])=((dt (f(x1)-2 f(x2)+f(x3)))/((Δx)^2))+(Subscript[f, old])

Finite Element Method

The finite element method is an alternative to the finite difference method for discretization of a real space.  The finite element method is used to solve the following equation for u[x] and f[x]:

u''(x)=f(x)

∫f[x]v[x]=∫u''[x]v[x]

∫f[x]v[x]=-∫u'[x]v'[x]=-ϕ[u,v]

∫f[x]v[x]=u[x]ϕ[(Subscript[v, i]),(Subscript[v, j])]

Matrix form:

M.f=L.u

v[x] can be any function, but it is chosen to be a basis, which represents the mesh [1].

Problem to Solve

A 1 mW laser beam with a 1.5 μm waist illuminates a semi-infinite silicon half-space.  What is the heating profile and temperature increase caused by the laser beam at the surface of the piece?  What is the frequency dependence of the heating at the surface (time constant)?

Heat Equation

The problem that I was interested in solving was the heat equation, which is a form of Poisson’s equation [2].

(dT/dt)=α(∇^2)T

Where T is temperature, α = k/cρ, k = thermal conductivity, c = specific heat capacity, ρ = density of material, and Q represents a heat source in terms of power.  As you can see, this problem should be solvable by both the finite difference and finite element methods.

FEM:

u''[x]=f[x]

u[x]=T[x]
f[x]=((dT[x])/dt)
((dT[x])/dt)=T''[x]

FDTD (as above)

(Subscript[f, new])=((dt (f(x1)-2 f(x2)+f(x3)))/((Δx)^2))+(Subscript[f, old])
(((Subscript[f, new])-(Subscript[f, old]))/dt)=(( (f(x1)-2 f(x2)+f(x3)))/((Δx)^2))
((dT[x])/dt)=T''[x]

Features of this Finite Difference Implementation

-Allows for asymmetric mesh/mesh refinement.
-Allows the use of arbitrary materials at arbitrary locations in terms of k,c, ρ for a given mesh cell.
-Allows the use of arbitrary time stepping.
-Can feature 1D or 2D time-varying sources.
-Features real-time monitoring.
-The entire code takes about 500kb.

Features of this Finite Difference Implementation

Monitor[For[i=1,i<=totali,i++,
buffer[[1]]=buffer[[1]]+tablefunc[i];
AppendTo[dt,Last[dt]+Δt[[i]]];
AppendTo[values,buffer[[1,Length[xmat]/2,1]]];
AppendTo[diff,Last[Differences[values]]/Last[Differences[dt]]];
If[i>5000&&Last[diff]<1,Break[],0];
],GraphicsGrid[{{ListLinePlot[Transpose[{dt[[2;;i]]-Δt[[1]],values[[2;;i]]-beginning}],PlotRange->All],
ListLogPlot[Transpose[{dt[[2;;i]]-Δt[[1]],diff[[2;;i]]}],PlotRange->All,Joined->True]},{ListLinePlot[Transpose[{ymat[[1;;plotr]]-ymat[[1]],((buffer[[1,Length[xmat]/2,1;;plotr]]-beginning))}],PlotRange->All,Mesh->All],ListLinePlot[Transpose[{xmat[[Length[xmat]/2]]-(xmat[[1;;Length[xmat]/2]]-xmat[[1]]),((buffer[[1,1;;Length[xmat]/2,1]]-beginning))}],PlotRange->All,Mesh->All]}},ImageSize->900]];

Mesh

-Largest mesh shown on left
-Smallest mesh shown on right

Border Source Implementation Method 1

In the first method, Dirichlet boundaries are placed on three edges, clamping them at RT (shown in yellow).  The source is placed *just outside* the simulation region (gradient at bottom), and the simulation develops in all but the bottom row.  The bottom boundary is set up so that no heat can escape through the source edge.

Border Source Implementation Method 2

In the second method, the three edges are still clamped at RT as before, but the source is placed inside the simulation region, and the simulation develops with the source inside of it, still including all but the bottom row.   The bottom boundary is again set up so that no heat can escape through the source edge.

Theoretical Solution

An analytical method is used to solve the problem in cylindrical coordinates [2].

(dT/dt)=(((d^2) T)/(dz^2))+((α d (dT r))/(dr dt r))

with an interface condition:

(dT/dz)=((Q[r])/k)

at

z=0

A Gaussian source is used:

Q[r]=2(P0/(π (R^2)))Exp[-((2(r^2))/(R^2))]

The equation above transforms into:

(((d^2)(SuperStar[T]))/(dz^2))-((γ^2)-(σ^2))(SuperStar[T])=0

The above ODE is solved and solution inverse Hankel transformed using Mathematica:

(SuperStar[f])[γ]=∫r f[r]BesselJ[0,γ r]\[DifferentialD]r,          f[r]=∫γ (SuperStar[f])[γ]BesselJ[0,γ r]\[DifferentialD]γ

Where R is the beam radius and P0 is the power of the beam in Watts.  This is the same source used in the simulation.  For an arbitrary layer structure, a numerical inverse Hankel transform method is solved with Mathematica. The homogenous case that I described previously has a steady-state solution with a magnitude at the surface of:

ΔT=(P0/((Sqrt[2π])R k))= 1.785 K

Convergence to Theory

When border source implementation method 1 is used and the mesh is refined further and further, the result converges to the analytical result:

Convergence to Commercial FEM Package

When border source implementation method 2 is used and the mesh is refined further and further, the result converges to that of a commercial software package:

Results - 2D Method 1

Theory result shown on left, FDTD result shown on right.

Results - 2D Method 2

Commercial software package result shown on left, FDTD result shown on right.

Results - Step Response

The step response at the point of highest temperature is shown below for the FDTD versus the commercial software package method.  These results are comparable, with time constants on the same order of magnitude.

Future Work

-Simulation of arbitrary devices/inhomogeneous media.
-Better utilization of nonuniform time stepping - perhaps using some form of Newton’s method to find the next “best” time step.
-Integration with System Modeler.
-Multiphysics modeling (Maxwell’s equations):

eqn2=(Eyx3-Eyx1)/dx-(Exy3-Exy1)/dy==-(Bznew-Bzold)/dt;
eqn4=(Bzy3-Bzy1)/dy==u0e0 (Exnew-Exold)/dt;
eqn5=-(Bzx3-Bzx1)/dx==u0e0 (Eynew-Eyold)/dt;
Expand[Solve[{eqn2,eqn4,eqn5},{Exnew,Eynew,Bznew}]]

layermatrixex[x_,y_,i_]:=Which[True,{{0,-Δt[[i]]/(Δy[[y]] mu eps[[x,y]]),0},{0,1,0},{0,Δt[[i]]/(Δy[[y]] mu eps[[x,y]]),0}}];
layermatrixey[x_,y_,i_]:=Which[True,{{0,0,0},{Δt[[i]]/(Δx[[x]] mu eps[[x,y]]),1,-Δt[[i]]/(Δx[[x]] mu eps[[x,y]])},{0,0,0}}];
layermatrixbz[x_,y_,i_]:=Which[True,{{0,-Δt[[i]]/Δy[[y]],0},{Δt[[i]]/Δx[[x]],1,-Δt[[i]]/Δx[[x]]},{0,Δt[[i]]/Δy[[y]],0}}];

extest={{0,exy1,0},{exx1,exold,exx3},{0,exy3,0}};
MatrixForm[extest];
eytest={{0,eyy1,0},{eyx1,eyold,eyx3},{0,eyy3,0}};
bztest={{0,bzy1,0},{bzx1,bzold,bzx3},{0,bzy3,0}};
funccey[ex_,ey_,bz_,x_,y_]:=Which[True,{{0,bz[[x,y-1]],0},{0,ey[[x,y]],0},{0,bz[[x,y+1]],0}}];
funccex[ex_,ey_,bz_,x_,y_]:=Which[True,{{0,0,0},{bz[[x-1,y]],ex[[x,y]],bz[[x+1,y]]},{0,0,0}}];
funccbz[ex_,ey_,bz_,x_,y_]:=Which[True,{{0,ey[[x,y-1]],0},{ex[[x-1,y]],bz[[x,y]],ex[[x+1,y]]},{0,ey[[x,y+1]],0}}];

eps={{eps1,eps1,eps1},{eps1,eps1,eps1},{eps1,eps1,eps1}};
Δx={Δx1,Δx1,Δx1};
Δy={Δy1,Δy1,Δy1};
Δt={Δt1};

Tr[layermatrixey[2,2,1].funccey[extest,eytest,bztest,2,2]]
Tr[layermatrixex[2,2,1].funccex[extest,eytest,bztest,2,2]]
Tr[layermatrixbz[2,2,1].funccbz[extest,eytest,bztest,2,2]]

Conclusion

-Demonstrated some of the capabilities of a commercial software package that may cost multiple thousands of dollars in a 400kb package.
-Can feature things such as mesh refinement, arbitrary time stepping, different-dimensional time varying sources, real time monitoring, arbitrary structures, etc, some of which I have used here.
-I have shown a possible cause for the difference between a commercial FEM package result and analytical theory result, which could be based on the placement of the heat source inside or outside of the simulation region.
-Multiphysics FDTD should be easy to implement with this method, given time away from other obligations.