Wafer Heating
 Adapted mesh Temperature 3d view of Temperature

// simulation of wafer heating by a laser.
/**************************************************************
From:  Fabian Dortu

Materials and Components Analysis group (MCA)
SPDT division

IMEC, Kapeldreef 75, B-3001 Leuven, Belgium
tel: +32/16 28 8774
e-mail: Fabian.Dortu@imec.be <mailto:Fabian.Dortu@imec.be>
**************************************************************/
// The poisson equation (here the heat equation) is solved in cylindrical coordinates:
//   $\Delta u = 1/r \partial_r(r \partial_r u) + \partial_{zz} u + 1/r^2 \partial_{\theta^2} u$
// so the variationnal formulation is after integer on the 3d cylinder
// rotation invariant given $\partial_{\theta} u =0$
//  $- \Delta u = f$
//  on the 1/2 plan $\theta =0$
//  $\int r ( \partial_r u \partial_r v + \partial_z u \partial_z v ) = \int r f v - \int_{\Gamma} r \partial_n u v$
//  the radial coordinate 'r' is called 'x'.
//  the depth coordinate 'z' is called 'y'.
// The wafer stuck is a rectangle with upper left corner at (0,0)
//  and lower right corner at (radius,-thickness)
// The laser beam it the surface from top to bottom at the center (0,0).
// The units assumed for the distance is the micrometer (um)

//***********************
//*** User parameters ***
//***********************

real radius=1000;   // the wafer radius (um)
real thick=600;     // the wafer thickness (um)
real thC=1.5e-4;    // the thermal conductivity (W/K/um)
real beamr=2;     // the laser beam radius (in the sens of a gaussian)
real refl=0.55;     // reflection coefficient Air/Silicon
real Eg=1.12;       // Band gap of Silicon (eV)
real En=1.49;       // Energy of laser (eV) - must be greater than Eg.
real I0=1e-2;       // Laser beam intensity (W/um^2)
real alpha=0.2;     // absorption coefficient (um^-1)
func f=I0 * (1-refl) * (En-Eg)/En * alpha * exp(alpha*y) * exp(-x*x/beamr/beamr);
// the heat generation function: a laser beam of radius 'beamr'
// lateral shape is gaussian.
// in depth shape is decreasing exponential because of absorption.
real g=300;          // temperature at right and bottom surfaces (in Kelvin)

//********************************
//*** Geometry/Mesh definition ***
//********************************

//--- less basic mesh ---
border bottom1(t=0,radius/5) {x=t; y=-thick; label=1;}
border right1(t=-thick,0) {x=radius; y=t; label=3;};  // 'right' actually means external surface
border top3(t=radius/5,0) {x=t; y=0; label=6;};
border left1(t=0,-thick/5) {x=0; y=t; label=7;}; // 'left' actually means center of cylinder
border left2(t=-thick/5,-thick) {x=0; y=t; label=8;}; // 'left' actually means center of cylinder
mesh Th = buildmesh ( bottom1(20) + bottom2(10) + right1(20) + top1(10) + top2(20) + top3(50) + left1(50) + left2(25) );

plot(Th,wait=1,fill=1,ps="wafer_mesh2.eps");
//savemesh(Th,"wafer_stuck.msh");

//**************************
//*** Problem definition ***
//**************************

fespace Vh(Th,P1);
Vh u,v,zero;

u=0;     // The temperature variable.
v=0;     // The weight function.
zero=0;  // used to set initial condition
int i=0; // variable used for mesh reconstruction
real error=0.1, coef=0.1^(1./5.);

//--- Variational Form ----
problem Problem1(u,v,solver=CG,init=i,eps=1.0e-8) =
int2d(Th) ( x*dx(u)*dx(v) + x*dy(u)*dy(v) )
+ int2d(Th) ( -v*x*f/thC )    // the source term = laser heating
+ on( 1,2,3,u=g ) ;             // fixed temperature at bottom and right surface.

//*************
//*** Solve ***
//*************

real cpu=clock();

cout << "***USER*** " << "Begin solve/adapt iterations" << endl;

for (i=0;i<10;i++)
{
cout << "***USER*** " << "Iteration number: " << i << endl;
real d = clock();
Problem1;
plot(u,wait=1);
cout << " CPU = " << clock()-d << endl;

error = error * coef;

} ;

plot(Th,wait=1,bb=[[0,0],[-10,-10]],nbiso=20,ps="adaptedmesh.eps"); // plot the last adapted mesh
plot(u,wait=1,bb=[[0,0],[-10,-10]],nbiso=20,ps="temperature.eps"); // plot the solution

//Plot a cut section at r=0
real[int] xx(10),yy(10);
for (i=0;i<10;i++)
{
x=0.; y=i/10.; // this line is used by the next one
xx[i]=i;
yy[i]=u; // value of u at point (0. , i/10.)
}

plot([xx,yy],ps="likegnu.eps",wait=true);

cout << " CPU = " << clock()-cpu << endl;


Page last modified on April 03, 2014, at 12:55 PM