[ff3d-users] a question about improving accuracy in ff3d

[Wei Huang]

Hi!

I just started to use ff3d for my research on thermal 
characterization of silicon ICs. I am using Gmesh for mesh 
generation and am modifying a sample code which I obtained 
from another researcher.

I found some discrepancies between ff3d and ANSYS for the 
same thermal FE analysis. I am wondering if there is 
anything I am doing wrong when using ff3d. Here is the case 
I am working on: I have a silicon bulk that is 20mm x 20mm 
x 0.5mm in dimensions. On the bottom surface, there is an 
area of 10mm x 10mm with a 2E+6 W/m^2 heat flux, and on the 
top surface is an isotherm of 0 degrees. I assume silicon's 
thermal conductivity to be a constant value of 100W/m-K. 
All other surfaces are adiabatic.

> From ANSYS, I got 10 degrees for the center of the 10mm x 
10mm heat source, ff3d also reports the same results 
(9.99997), everything is perfect so far. But at the middle 
points of four edges of the 10mm x 10mm heat source area, I 
got 5.42122, 5.44859, 5.30292, 5.00623 from ff3d, whereas 
in ANSYS, I got 5.0 (+/-0.1 at most). Since the correct 
temperature on the edges should be (10+0)/2=5, I am 
wondering what should I do to reduce the error in ff3d? I 
have tried to generate a finer mesh, more iterations, but 
without success. Any suggestions will be appreciated! Below 
are my ff3d and gmesh files:

(PS: I used "gmsh -o test.msh -3 -optimize test.geo" to 
generate the mesh, also tried two values--0.0005 and 
0.001--for "characteristic length" of target element sizes 
in Gmesh. I am also running gmesh with 0.00005 for the 
fourth element in the "Point" definition, but it seems to 
take forever to finish on my PC...)


FF3D:
////////////////////////////
double q1 = 0 ;
double q2 = 0 ;
double q3 = 2000000 ;
double q4 = 0 ;
double q5 = 0 ;

double k1 = 100.000000;
double z1 = 0.000500;

function  uTop  = 
one((x>=-0.01)&&(x<=0.01)&&(y>=0.005)&&(y<=0.01));
function  uLeft  = 
one((x>=-0.01)&&(x<=-0.005)&&(y>=-0.005)&&(y<=0.005));
function  uCenter  = 
one((x>=-0.005)&&(x<=0.005)&&(y>=-0.005)&&(y<=0.005));
function  uRight  = 
one((x>=0.005)&&(x<=0.01)&&(y>=-0.005)&&(y<=0.005));
function  uBottom  = 
one((x>=-0.01)&&(x<=0.01)&&(y>=-0.01)&&(y<=-0.005));

function qTop   =    q1 * uTop  ;
function qLeft   =    q2 * uLeft  ;
function qCenter   =    q3 * uCenter  ;
function qRight   =    q4 * uRight  ;
function qBottom   =    q5 * uBottom  ;

function q = qTop + qLeft + qCenter + qRight + qBottom ;


function f1 = one((z>=0)&&(z<z1));
function fk = k1*f1 ;

// Load a external mesh build with GMSH
mesh M = read(gmsh,"test.msh");

solve(u) in M 
krylov(type=cg,precond=ichol),cg(maxiter=100000000,epsilon=1e-8),memory(matrix=sparse)
{
 test(w)
// Variational form of the Heat equation
// int[2000075]*(w*q) is the boundary condition in the 
bottom face   (Plane surface number 2000075 in test.geo 
file)

                int(fk*grad(u)*grad(w))   =   int[2000075](w*q);
                u = 0 on 2000077;

};

cout << u(0,-0.005,0)<<"\n";
cout <<"Top  "<< "Left "<<"Center   "<< "Right   "<<"Bottom 
"<<"\n";
cout << u(0,0.005,0)<< " ";
cout << u(-0.005,0,0)<< " ";
cout << u(0,0,0)<< " ";
cout << u(0.005,0,0)<< " ";
cout << u(0,-0.005,0)<< " ";

;
///////////

GMESH geometry file (test.geo):
-------------

Point(1) = {-0.01, -0.01, 0, 0.0005};
Point(2) = {0.01, -0.01, 0, 0.0005};
Point(3) = {0.01, 0.01, 0, 0.0005};
Point(4) = {-0.01, 0.01, 0, 0.0005};

Point(5) = {-0.01, -0.01, 0.0005, 0.0005};
Point(6) = {0.01, -0.01, 0.0005, 0.0005};
Point(7) = {0.01, 0.01, 0.0005, 0.0005};
Point(8) = {-0.01, 0.01, 0.0005, 0.0005};

Line (1) = {4, 3};
Line (2) = {3, 2};
Line (3) = {2, 1};
Line (4) = {1, 4};

Line (5) = {8, 7};
Line (6) = {7, 6};
Line (7) = {6, 5};
Line (8) = {5, 8};

Line (9) = {1, 5};
Line (10) = {2, 6};
Line (11) = {3, 7};
Line (12) = {4, 8};

Line Loop (3000075) = {1, 2, 3, 4};
Plane Surface (2000075) = {3000075};
Line Loop (3000077) = {5, 6, 7, 8};
Plane Surface (2000077) = {3000077};
Line Loop (3000079) = {6, -10, -2, 11};
Plane Surface (2000079) = {3000079};
Line Loop (3000081) = {3, 9, -7, -10};
Plane Surface (2000081) = {3000081};
Line Loop (3000083) = {9, 8, -12, -4};
Plane Surface (2000083) = {3000083};
Line Loop (3000085) = {12, 5, -11, -1};
Plane Surface (2000085) = {3000085};

Surface Loop (3000087) = {2000075, 2000081, 2000079, 
2000077, 2000085, 2000083};
Volume (2000087) = {3000087};

------------------------