...

Mak-FinalReport.pdf

by user

on
Category: Documents
4

views

Report

Comments

Transcript

Mak-FinalReport.pdf
Modeling Transient Liquid Phase Bonding Process
in the Ni-P System
by
Billy Mak
An Engineering Project Submitted to the Graduate
Faculty of Rensselaer Polytechnic Institute
in Partial Fulfillment of the
Requirements for the degree of
MASTER OF ENGINEERING
Major Subject: MECHANICAL ENGINEERING
Approved:
_________________________________________
Professor Sudhangshu Bose, Primary Project Adviser
_________________________________________
Professor Ernesto Gutierrez-Miravete, Secondary Project Adviser
Rensselaer Polytechnic Institute
Hartford, CT
December 2010
(For Graduation January 2011)
i
© Copyright 2010
by
Billy Mak
All Rights Reserved
ii
CONTENTS
Modeling Transient Liquid Phase Bonding Process ........................................................... i
in the Ni-P System .............................................................................................................. i
NOMENCLATURE ......................................................................................................... vi
LIST OF TABLES ........................................................................................................... vii
LIST OF FIGURES ........................................................................................................ viii
ACKNOWLEDGMENT .................................................................................................. ix
ABSTRACT ...................................................................................................................... x
1. INTRODUCTION ....................................................................................................... 1
1.1
TLP Bonding Process ......................................................................................... 1
1.2
Types of Modeling ............................................................................................. 2
1.3
Problem Statement ............................................................................................. 3
1.3.1
Methods .................................................................................................. 4
2. THEORY AND METHODOLOGY ........................................................................... 6
2.1
Diffusion ............................................................................................................ 6
2.2
Assumptions and Initial Conditions ................................................................... 6
2.3
Fixed-Grid Numerical Modeling........................................................................ 7
2.4
Ikawa’s Method .................................................................................................. 8
2.5
2.6
2.4.1
Equations ................................................................................................ 8
2.4.2
Initial Design Parameters ....................................................................... 9
Forward Finite Difference ................................................................................ 10
2.5.1
Equations .............................................................................................. 10
2.5.2
Initial Design Parameters ..................................................................... 11
COMSOL Theory and Design ......................................................................... 12
2.6.1
Governing Equations ............................................................................ 12
2.6.2
Initial Design Parameters ..................................................................... 13
3. RESULTS AND DISCUSSION ................................................................................ 17
iii
3.1
Li’s Method ...................................................................................................... 17
3.2
Ikawa’s Method ................................................................................................ 18
3.3
3.4
3.5
3.2.1
Program 1 ............................................................................................. 18
3.2.2
Program 2 ............................................................................................. 20
Forward Finite Difference ................................................................................ 23
3.3.1
Case Study 1 – Divergence of Cs ......................................................... 25
3.3.2
Behavior of the Liquid Layer Half Width ............................................ 25
3.3.3
Effect of Arbitrary Reference Value of the Activity [a*] and
Reciprocals of Slopes on Binary Phase Diagram [ks] .......................... 25
COMSOL ......................................................................................................... 26
3.4.1
Effect of Infinitesimal Difference [  ] .............................................. 27
3.4.2
Behavior of the Liquid Layer Half Width ............................................ 28
3.4.3
Behavior of the Effective Concentration Values ................................. 28
3.4.4
COMSOL Limitations.......................................................................... 28
Overall Discussion ........................................................................................... 29
4. CONCLUSIONS ....................................................................................................... 30
4.1
Future Work ..................................................................................................... 32
5. REFERENCES .......................................................................................................... 33
6. APPENDIX A – COMSOL Reference ...................................................................... 35
6.1
Developing the Model: Initial Concentration [cini(x)], Effective Reciprocal
Activity Coefficient [  eff (Ceff)] ....................................................................... 35
6.1.1
6.2
Sample Calculation for Effective Concentration [Ceff,b] and Effective
Reciprocal Activity Coefficient [  eff ] ................................................. 36
COMSOL Data ................................................................................................ 37
6.2.1
Models .................................................................................................. 37
6.2.2
Output Data .......................................................................................... 37
7. APPENDIX B – FORTRAN: Ikawa Method ............................................................ 38
7.1
First Program .................................................................................................... 38
iv
7.2
7.1.1
Code ..................................................................................................... 38
7.1.2
Input Data ............................................................................................. 40
7.1.3
Output Data .......................................................................................... 40
Second Program ............................................................................................... 40
7.2.1
Code ..................................................................................................... 40
7.2.2
Input Data ............................................................................................. 41
7.2.3
Output Data .......................................................................................... 41
8. APPENDIX C – FORTRAN: Forward Finite Difference Method ............................ 42
8.1
Code ................................................................................................................. 42
8.2
Input Data ......................................................................................................... 44
8.3
Output Data ...................................................................................................... 44
9. APPENDIX D – Ikawa Case Studies ........................................................................ 46
9.1
Case Study 1: Effect of Diffusivity .................................................................. 46
9.1.1
9.2
Case Study 2: Liquid Layer Half Width Change [h=2.5x10-5 m] .................... 47
9.2.1
9.3
Data Files ............................................................................................. 46
Data Files ............................................................................................. 48
Case Study 3: Molar Volume Change .............................................................. 48
9.3.1
Data Files ............................................................................................. 49
10. APPENDIX E – Forward Finite Difference Case Studies......................................... 50
10.1 Case Study 3: Concentration of Solid Input Parameter [Cs=0.16%] ............... 50
10.1.1 Data Files ............................................................................................. 50
v
NOMENCLATURE
a
chemical activities
a*
arbitrary reference value of the activity
C
solute atom concentration: COMSOL: m3/mol, % at of P
D
diffusivity, m2/s
D0
diffusion coefficient of the material, m2/s
G
reaction rate, mol/(m3-s)
h
initial liquid layer half width of the crack, m
J
flux, mol/(m2-s)
k
reciprocals of slopes on binary phase diagram
L
length of system, m
Q
activation energy, J
R
gas constant, J/K-mol
t
time, s
T
temperature, K
V
molar volume, m3/mol
x
distance along the system (position coordinate), m
dt
time step, s
dx
element position size, m
erf
error function
hmax
maximum theoretical liquid layer half width with constant molar volume, m
Greek Symbols

activity coefficient
 ts
time scaling coefficient

infinitesimal difference

arbitrary value to calculate moving boundary
Superscripts/Subscripts
*
arbitrary reference number
’
parameterized value
l
liquid phase at the solid-liquid interface
s
solid phase at the solid-liquid interface
b
base metal
i
initial condition at liquid
0
constant for material
o
value of [base metal position – initial condition]
SL
solid-liquid layer
vi
LIST OF TABLES
Table 1 Initial Data for Ikawa Program #1 ....................................................................... 9
Table 2 Initial Data for Ikawa Program #2 ..................................................................... 10
Table 3 Initial Data for Forward Finite Difference Program .......................................... 11
Table 4 Properties in Sub-Domain, Meshing, and Time Steps....................................... 13
Table 5 Function: cini(x) for various  ........................................................................ 14
Table 6 Function:  eff (ceff) for various  .................................................................... 15
Table 7 Liquid Region Closure Time Summary............................................................. 29
Table 8 Initial Data: Effective and Actual Concentration for Various  .................... 35
vii
LIST OF FIGURES
Figure 1 Ni-P binary phase diagram [10] ......................................................................... 6
Figure 2
(a) Schematic Diagram of Ni-P Binary Alloy System, and (b)/(c)
Concentration Profile of P at Bonding Region [9] ............................................................ 8
Figure 3 Schematic Diagram of COMSOL Model .......................................................... 13
Figure 4 Liquid Layer Half-Width vs. Bonding Time for Li’s Model [1] ..................... 17
Figure 5 Concentration vs. x for Ikawa Program 1 ........................................................ 19
Figure 6 Liquid Layer Half Width vs. Time for h=1.25x10-5 m in Ikawa ...................... 21
Figure 7 Liquid Layer Half Width vs. log Time for h=1.25x10-5 m in Ikawa................ 21
Figure 8 Liquid Layer Half-Width vs. Time for a*=340 in Forward Finite Difference . 24
Figure 9 Liquid Layer Half-Width vs. log Time for a*=340 in Forward Finite Difference
......................................................................................................................................... 24
Figure 10 Liquid Layer Half-Width vs. log Time for a* in Forward Finite Difference . 26
Figure 11 Liquid Layer Half-Width vs. log Time for   0.006%,0.004% ............... 27
Figure 12 Liquid Layer Half Width vs. log Time – All the Methods ............................. 29
Figure 13 Actual vs. Effective Concentration for   0.006% ................................... 36
Figure 14 Actual vs. Effective Concentration for   0.004% .................................... 36
Figure 15 Concentration vs. x for Varying Diffusivity .................................................. 46
Figure 16 Liquid Layer Half-Width vs. Time for Varying Diffusivity .......................... 47
Figure 17 Liquid Layer Half-Width vs. Time for h=2.5x10-5 m in Ikawa ...................... 47
Figure 18 Liquid Layer Half-Width vs. log Time for h=2.5x10-5 m in Ikawa ............... 48
Figure 19
Liquid Layer Half-Width vs. log Time for Varying Molar Volume at
-6
h=12.5x10 m in Ikawa ................................................................................................... 49
Figure 20 Liquid Layer Half-Width vs. log Time for Cs=0.16% in Forward Finite
Difference ........................................................................................................................ 50
viii
ACKNOWLEDGMENT
I would like to thank my family, friends, employer and professors for their support while
I worked on the Masters program.
This topic idea originated with Professor Ernesto Gutierrez-Miravete as one of his
research case studies. When I initially started the project, I had not heard of transient
liquid phase (TLP) bonding. I was entering the project with a slight knowledge of
diffusion from my undergraduate courses. Professor Gutierrez-Miravete was extremely
helpful throughout the project timeline.
Professor Sudhangshu Bose aided in creating a distinguishable final piece for the project.
For all their help and support, I am specially thanking Professor Gutierrez-Miravete and
Professor Bose as they advised me throughout the project.
ix
ABSTRACT
The main purpose of this project was to investigate the transient liquid phase (TLP)
bonding process through numerical modeling to develop a qualitative understanding of
the process and a quantitative characterization of the relationships that exist among
process parameters such as concentration distribution in the material system, time and
temperature. A material system of nickel phosphorous, Ni-P superalloy (19% atomic
percent P) with a pure nickel, Ni base metal was used in this study. A detailed review of
the relevant literature was first carried out. Then three different modeling methods were
investigated: an exact solution, a finite difference method and a finite element approach
using the COMSOL Multiphysics software.
The results obtained from the three
modeling methods were compared with previously published results and also among
themselves. All three models indicate that the liquid region ultimately disappears under
isothermal conditions when held for sufficiently long treatment times as P atoms diffuse
out of the region and into the base metal. However, the numerical approaches predict
initial growth of the liquid region before its disappearance.
While the qualitative
behavior predicted by all three methods was in agreement with prior work, there were
some discrepancies in the predicted liquid region closure times. The finite element
modeling approach needs to be further investigated.
.
keyword(s): TLP bonding process, diffusion-controlled process, Ikawa, Ni-P, repairing
techniques, joint healing, superalloy
x
1. INTRODUCTION
TLP bonding method is a joining process applied to metallic systems. The process heals
and repairs component imperfections through heat treatment. The most common
repairing material is superalloy, a set of high temperature capable alloys of Ni or Co.
The repairing technique produces a strong, interface-free joint with no residue of the
repairing superalloy. Major applications for TLP bonding is in aerospace, specifically
aircraft engines and rockets. Equipments such as jet engines components undergo
extreme thermal and stress cycling, and are subjected to significant deterioration.
Examples of such deterioration are erosion, cracking, corrosion, wearing and damage
due to external factors. If the engine components are inspected and declared to be
damaged, they must either be repaired or replaced. Usually these damages only manifest
as small cracks in the component, but will prove to be a liability if not repaired.
1.1 TLP Bonding Process
The TLP bonding process requires the superalloy to have the same chemistry as the
material component being repaired. The superalloy is blended with a pliable polymer
binder and is smeared onto the crack of the component where it is subjected to a heat
treatment at constant high temperature under a controlled inert atmosphere. The blend
melts under the heat creating a thin liquid layer (interlayer) within the crack, being held
by capillary action. The atoms diffuse into the component through the crack wall at a
diffusion rate of the given material. Not to be confused with diffusion bonding, the thin
liquid layer eliminates the need of high bonding or clamping force that is required in
diffusion bonding. As the atoms diffuses into the component, the phosphorus content in
the interlayer decreases causing the melting point to exceed the heat treatment
temperature, solidifying the liquid.
After a period of time, the content solidifies
following the microstructures of the component at the cracked wall. When the process is
complete, the interlayer becomes indistinguishable from the materials of the component
and the performance of the component is restored to be comparable to that of a brand
new component.
The TLP bonding method behaves as a diffusion process under a controlled atmosphere
subjected to extreme temperature conditions. Therefore, the liquid layer’s thickness,
1
time and temperature are important parameters in TLP bonding process, as they are for
natural diffusion. The liquid layer’s thickness determines liquid region closure time, as
diffusion is dependent on the total mass that needs to be transferred. Similarly, an
increase in temperature for the controlled atmosphere will increase the mass transfer
rate, resulting in a decrease of the process duration.
1.2 Types of Modeling
In literatures, modeling of TLP bonding conditions and superalloy compositions were
performed analytically and numerically. Many methods utilize the concept of Fick’s
diffusion equation to provide insight into the mechanisms of TLP bonding, and are used
to optimize the bonding conditions and compositions. However, analytical models are
used less frequently than numerical models due to the complexities of the analytical
approach. Also, analytical models are typically designed to be ideal and impractical for
actual applications, as they assume constants to keep the degrees of freedom to a
minimum. Numerical modeling is a time effective method, compared to analytical
modeling, for visualization of the relationships among initial design parameters. The
user only has to change the input data to get results for conditions they are trying to
pursue.
Although there are many articles in the literature on the modeling of the TLP bonding
process, each with their unique approach, the project of this report will only address
selected literatures.
Ikawa et al. [8,9] utilized theoretical and experimental approaches to clarify the bonding
mechanism and to obtain bonding conditions to ensure joint performance for Ni-base
superalloys. In his method, the solidification process of the liquid film at the bonding
region is theoretically determined to be controlled by the diffusion process [8].
Li et al. [1] approached the process through a volume controlled finite difference
method. The volume controlled finite difference method is a numerical method widely
used in fluid flow applying the concepts of conventional finite difference method with
the finite element method. The numerical model was designed to a centered singlemeshed fixed-grid source based method using the diffusion equation to characterize and
simulate the different phases and moving interfaces through implicit time integration.
2
The fixed-grid source-based method was originally developed to simulate the
temperature fields for melting-solidification phase change processes.
The main
objective of the work was to demonstrate the validity and effectiveness of the method in
comparison with previous models, such as those of Illingworth et al. [7] and Zhou et al.
[12].
COMSOL Multiphysics is a finite element analysis, solver, and simulation
software/FEA software package for various physics and engineering applications [11].
The current project is modeled through COMSOL’s diffusion application mode as a onedimensional (1-D) model to simulate the behavior of mass transport. To the best of the
author’s knowledge, the COMSOL approach for TLP bonding process has not been
previously performed and published in literatures. The COMSOL model is verified by
comparison with the models presented in selected literatures.
1.3 Problem Statement
The current project performed a finite element analysis of the TLP bonding process by
using FORTRAN and COMSOL to explore the qualitative characterization of the
process and the quantitative relationships that exist among the following process
parameters:

Temperature

Time

Width of crack joint in a base metal

Elemental concentration profile
The finite element analysis is compared with the approximation methods from selected
literatures listed below:

A theoretical approach using the method established by Ikawa et al. [8, 9]

A numerical approach using a forward finite difference method under an explicit
scheme, and,

A fixed-grid numerical modeling by Li et al. [1]
3
1.3.1
Methods
The project explored the methods of selected literatures, with the exception of
COMSOL. The COMSOL model for TLP bonding process is currently not presented in
the literature.
The different methods from the literatures and COMSOL are then
compared to obtain a characterization of the bonding behavior and the relationships
among the process parameters.
1.3.1.1 Li’s Method
The forward finite difference method and Ikawa’s theoretical approximation model were
performed using FORTRAN coding subjected to Li’s input data [1]. The results in Li’s
fixed-grid source-based mesh were the basis of comparison between FORTRAN and
COMSOL modeling.
1.3.1.2 Ikawa’s Method
Ikawa’s equations are the basis for developing two FORTRAN programs to characterize
the behavior of the TLP process. The first program models the behavior of the solid and
liquid interface. The results show the concentration distribution for diffusion in a semiinfinite medium at various time durations, C vs. x.
The second program models the concentration behavior of the liquid layer width at
various time durations of the TLP bonding process. The results show the evolution of
the liquid layer half width of the crack at various times for a specified temperature range,
xSL vs t.
The programs were provided by Professor Gutierrez Miravete and are modified to fit the
requirements of the project.
1.3.1.3 Forward Finite Difference Method
Professor Gutierrez-Miravete, one of my project advisors, performed academic research
on the TLP bonding process and has provided, for the completeness of the project, a
FORTRAN program implementing forward finite different method under an explicit
scheme for TLP bonding process. In the model, a mesh was first introduced and the
diffusion equation was discretized using the standard second order accurate finite
4
difference formula for the space derivative and a forward finite difference formula for
the time derivative. The helped results model the concentration behavior of the liquid
layer width at different durations of the TLP bonding process, xSL vs. t.
1.3.1.4 COMSOL Method
COMSOL modeling is a numerical approach that has not been performed in previous
TLP studies. In order to create a model, a characterization of the TLP bonding process
was obtained from modeling based on existing literatures. Therefore, the project
explored whether the COMSOL modeling was capable of simulating the TLP bonding
process. The overall COMSOL results for the project modeled the concentration
behavior of the liquid layer width at different durations of the TLP bonding process, x SL
vs. t. It also monitored the behavior of the solid-liquid layer interface.
5
2. THEORY AND METHODOLOGY
2.1 Diffusion
Many numerical models of the TLP bonding process follow the atomic level of natural
diffusion taking into consideration temperature and time.
At higher temperatures,
diffusion occurs at a faster rate, as demonstrated by the diffusion coefficient equation, D.
D   D0 e

Q
RT
.
(Equation 1)
2.2 Assumptions and Initial Conditions
Figure 1 Ni-P binary phase diagram [10]
The project’s focus is the system consisting of a Ni-P superalloy and a base metal of
pure Ni. Under isothermal conditions of high temperature in a control volume, the P
atoms diffuse into the base metal of pure Ni. The concentration is limited to the solute
characteristics of the Ni-P binary phase diagram, Figure 1. More specifically, the project
used Ni-P with a 19% atomic or 11.02% weight percentage of P at an initial temperature
of T=1473K. At this temperature, the liquid solution is super-saturated with a Cl and Cs
values of 10.223% and 0.166%, respectively. As the P atoms diffuse into the base metal,
6
the solid-liquid interface moves towards the inwards into the liquid layer shown in
Figure 2(b) and (c), and the width of the crack becomes thinner and thinner. As the
system becomes fully diffused, the concentration of P in the liquid is equal to that of the
solid, creating a repaired yet "identical" component as long as the base metal is
concerned.
To develop the model of the system, the following overall assumptions were introduced:
1. One dimensional (1-D) diffusion-controlled model
2. At the solid-liquid interface, the concentration gradient of P in the solid phase is
assumed to be constant even if the solid-liquid interface moves
3. Only phosphorus diffusion is calculated due to its rapid diffusivity
4. The diffusion coefficient is independent of composition
5. The crack is in the form of a wide deep slot so that only the diffusion of the
phosphorus perpendicular to plane of crack and into component is considered
6. Displacement of the solid-liquid interface (moving boundary layer) occurs
The model also has the following initial conditions:
1. The crack is symmetrical. The region for analysis will be the region from the
center of the crack (x=0).
 C

2. Concentration in the material inside the crack is uniform 
 0  and the
 x

concentration in the surround material (far away from the crack) is specification
value (C  C ) .
These initial conditions and assumptions are used in all methods.
2.3 Fixed-Grid Numerical Modeling
An in-depth description of the TLP bonding process using a fixed grid numerical model
was explored in Li et al [1]. His work provides a complete formulation of the method
using mathematical modeling. His model involves a 1-D diffusion-controlled, solidliquid phase and moving interface. The results were only used to validate and for
comparison with COMSOL, Ikawa, and forward finite difference models. The input
parameters for the models were provided by Li’s input parameters for Ni/Ni-P.
7
2.4 Ikawa’s Method
The following is a summary of the major concepts and equations used in the FORTRAN
programs for Ikawa’s method. His methodology is further explained in his publications
[8].
2.4.1
Equations
The equation formulations start by introducing the chemical flux of P atoms through the
solid-liquid interface, J,
 D  dC 
J   DC   

 Vs  dx 
(Equation 2)
As the solid-liquid interface moves dx in the time interval dt due to P atoms diffusing
into the solid phase, J is modified and expressed by the following:
 dx  C C 
J    l  s 
 dt  Vl Vs 
(Equation 3)
Where Ci and Cs are values obtained from the Ni-P phase diagram [10] as displayed in
the schematic diagram shown in Figure 2(a).
Figure 2 (a) Schematic Diagram of Ni-P Binary Alloy System, and (b)/(c) Concentration Profile of P
at Bonding Region [9]
Assuming that the condition of the solid phase at the solid-liquid interface remains
constant, the P concentration at an arbitrary position x in the solid phase is expressed by
the following:
C C

Vl Vl

 x  h   
* 1  erf 
 
 2 Dt  

(Equation 4)
8
and at x=h
C
 dc 
   s
Dt
 dx  x  h
(Equation 5)
Using Equation (2), (3), and (5) at an arbitrary time t, the thickness of the liquid film, 2x,
is expressed as the following:
 4Cs
2 x  2h  
 Vs 
  Cl Cs 
*  
 V V 
s 
  l
1
Dt
(Equation 6)
Under the 1-D model case, solute atom concentration in the bonding layer at arbitrary
position and time, x and t, is expressed as the following:
C
C ( x, t )  Co   o
 2
 x  h  
  h  x 
 erf 
  erf 
 
  2 Dt 
 2 Dt  
(Equation 7)
And solute atom concentration at the center of the bonding layer is expressed by the
following:

 h   
C (t )  Co 1  erf 
   Ci
2
Dt



(Equation 8)
Equations (6), (7), and (8) are the main equations necessary to develop the FORTRAN
programs utilizing the Ikawa method to model the solid-liquid layer interface, and
concentration behavior of the liquid layer width at different durations of the TLP
bonding process.
2.4.2
Initial Design Parameters
Table 1 and 2 represent the design parameters for the FORTRAN programs. The design
parameters are the input data for running the programs.
Table 1 Initial Data for Ikawa Program #1
Cl
Cs
H
D
V
Dx
10.223%
0.166%
0.0000125 m
1.8x10-11 m2/s
7.344x10-6 m3/mol
0.000001 m
9
Dt
3600 s
0.001

Although Li prescribed a Ci = 19%, Ikawa’s program utilizes Cl as its Ci. Also,  is an
arbitrary value used to track precisely the moving boundary (the solid-liquid layer
interface) for the first program.
Table 2 Initial Data for Ikawa Program #2
Cl
Cs
H
D0
Q
V
Tlow
Thigh
Dt
10.223%
0.166%
0.0000125 m
0.0228 m2/s
248000 J
7.344x10-6 m3/mol
1423 K
1523 K
150 s
The diffusivity, D0, was determined from using known values[1, 8, 12] of D, Q, R, and T
at 1473K, and the second program’s temperature range for analysis is from Tlow=1423K
to Thigh=1523K.
2.5 Forward Finite Difference
The following is a summary of the major concepts and equations used in the FORTRAN
program for the forward finite difference method.
2.5.1
Equations
The governing equation for the forward finite difference in TLP bonding follows the
transport equation:
C  2u
 2
t
x
(Equation 9)
where C(x,t) is the solute atom concentration in the bonding layer at arbitrary time and u
is defined as
10
'

Dl (a  a*) if a  al  a *
u '

Ds (a  a*) if a  a s  a *
(Equation 10)
a * is an arbitrary reference value of the activity. al and a s are chemical activities and
are expressed as follows:
a  al  
C
Cl
(liquid) and a  a s   s (solid)
kl  l
k s s
(Equations 11, 12)
where k l and k s are obtained from the phase diagrams for liquidus and solidus phase
slopes from Figure 1,
kl 
1
1
and k s 
,
dTl
dTs
dC
dC
(Equations 13, 14)
Dl' and Ds' are the modified diffusion coefficients,
Dl'  Dl k l  l and Ds'  Ds k s  s .
(Equations 15, 16)
The diffusion equation is discretized using the standard second order accurate finite
difference formula for the space derivative and a forward finite difference formula for
the time derivative.
Since the method is explicit, stability of the computation restricts the size of the time
step. The maximum time step is expressed by the Courant-Fridrich-Lewy condition
t min 
x 2
(Equation 17)
2Dmax
These equations are the basis for the development of the FORTRAN program.
2.5.2
Initial Design Parameters
Table 3 represents the initial input design parameters (input data files) for the
FORTRAN program.
Table 3 Initial Data for Forward Finite Difference Program
Dl
Ds
kl
500x10-12 m2/s
18x10-12 m2/s
-0.02684
11
ks
Cl
Cs
Ci
Cb
a*
Xil
Xtotal
ttotal
-0.0014
10.223%
0.5%
19%
0%
340
0.0000125 m
0.003 m
40000 s
2.5.2.1 Effect of a* and ks
Under the forward finite difference model, the slope of the solidus line on the Ni-P
binary phase diagram is an initial data parameter that was not referenced in Li,
Illingworth, Ikawa, or Zhou and North [1, 7, 8, 9, 12]. Since an exact value was not
available, an approximation was made. The value used affected the program’s input
value for ks and a*, as both depend on an exact value for a*. To complete this project,
an analysis of the varying a* was analyzed to determine the behavior on the bonding
process. Since a* is an arbitrary reference value of the activity, and based on Equation
10, a* was assumed to have a value close to as. Assuming the value of ks in Table 3,
small values of a* ranging from 160 to 353.5 were used in the calculation. The allowed
maximum theoretical liquid layer half width, hmax, was 2.32x10-5m from Li et al. As a
reference, the corresponding a* was determined to match the results in selected
literatures.
2.6 COMSOL Theory and Design
The simulation in COMSOL was an implicit 1-D diffusion application model under
transient analysis.
2.6.1
Governing Equations
The COMSOL model follows the generic diffusion equation which has the same
structure as the heat conduction equation and is governed by the following mass transfer
equation (without convection):
 ts
C
  DC   G
t
(Equation 18)
12
Assuming the G=0, the equation simplifies to
 ts
C
  DC   0
t
(Equation 19)
This equation is similar to Equation 9, the governing transport equation for the TLP
bonding process.
2.6.2
Initial Design Parameters
The initial design parameters used for COMSOL were dependent on the data provided in
Li et al [1] and the assumptions listed in Section 2.2. Table 4, 5, and 6 represent the
initial input design parameters for the COMSOL program.
A schematic diagram of the 1-D model is shown in Figure 3 as a cross-section along the
x-direction, referencing a symmetrical assumption (liquid layer half-width). The system
is separated into three parts, the liquid layer, the solid-liquid layer interface, and the solid
layer, with two boundary conditions located at x=0 and x=L, identifying that the system
is insulated.
Figure 3 Schematic Diagram of COMSOL Model
2.6.2.1 COMSOL Properties
Table 4 Properties in Sub-Domain, Meshing, and Time Steps
 eff (ceff)
Time Scaling Coefficient,  ts
Diffusion Coefficient, D
Reaction Rate, G
Initial Concentration, c(t0)
Initial Liquid Layer Half Width, h
Total Length of System [from (1) to (2)], L
13
1.80x10-11 m2/s
0 mol/(m3-s)
cini(x)
12.5x10-6 m
250x10-6 m
5.00x10-7 m
0.1 s
50000 s
Maximum Element Size, dx
Maximum Time Stepping Size, dt
Total Time
The definition of the functions,  eff (ceff) and cini(x), are addressed in Section 2.6.2.2.
The total length of the system provided in COMSOL was significantly less than the
value applied in the other methods. Through testing of the other models prior to the
testing of COMSOL, Figures 5 through 9 demonstrate that the only requirement was for
the total length was to be a value such that the system was significantly larger than the
initial liquid layer half width. It did not need to meet the requirement specified by Li et
al. By performing this cut in the properties, it allowed a significant reduction of the
computational processing time to run the COMSOL model. This advantage will allow a
finer mesh on either or both the time and element spacing.
2.6.2.2 Functions
The time scaling coefficient,  ts , and the initial concentration distribution, c(t0), were
created as linear piecewise functions defined by Table 5 and 6.
2.6.2.2.1 Function: cini(x)
Interpolation Method: Linear
Table 5 Function: cini(x) for various 
x, m
0
1.25x10-5
1.251x10-5
2.5x10-4
 = 0.006%
cini(x)
0.08949
0.08949
0
0
 = 0.004%
cini(x)
0.08947
0.08947
0
0
Under the initial conditions, the liquid layer has 100% of the P, corresponding with
Ci=19% and the solid layer is pure Ni (0% of the P or Cb=0%). The calculation for cini(x)
is discussed in Appendix A.
14
2.6.2.2.2 Function:  eff (ceff)
Interpolation Method: Linear
Table 6 Function:  eff (ceff) for various

 = 0.006%
 eff (ceff)
ceff
0
-3
1.59x10
1.60x10-3
1.72x10-3
1.73x10-3
8.95x10-2
 = 0.004%
 eff (ceff)
ceff
1
1
838.08
838.08
1
1
0
-3
1.61x10
1.62 x10-3
1.70 x10-3
1.71 x10-3
8.95x10-2
1
1
1257.13
1257.13
1
1
The function  eff (ceff) calculates the effects of the time scaling coefficient identified by
the COMSOL program.
Typically a time scale of 1 is used for general diffusion
problems, but under the case of the project, the solid-liquid interface is considered to be
a very thin film compared to the solid and liquid region. Therefore, in order to account
for the scenario of the solid-liquid layer interface,  eff (ceff) is developed. The values
obtained for  eff (ceff) are explained in Appendix A.
2.6.2.3 Sub-domain and Boundary Conditions
The sub-domain properties utilize the functions  eff (ceff) and cini(x) to describe the time
scaling coefficient and the initial concentration.
The functions are mainly a
transformation to represent the solid-liquid layer interface, the separation between the
solid (Ni) and liquid (superalloy Ni-P) layer, and the initial concentration distribution.
As mentioned before, the boundary conditions are from point 1 (x=0m) to point 2 (x=L),
which refer to the middle of the crack to the end of the base metal, representing the total
length. Since the system is under a controlled environment both boundary conditions are
insulated points.
15
2.6.2.4 Meshing/Time Steps
To mesh the model, the free mesh parameters option was used and meshed according to
the maximum element size. The maximum element size determines the minimum
number of nodes for the model. The element spacing must minimize the computer
memory’s usage while making certain the values obtained are reliable. With an initial
crack width of h=12.5x10-6m, a maximum element size of dx=0.05x10-6m, appears to be
a suitable value. This maximum element size corresponds to 5000 elements. A smaller
value of the maximum element size will result in more computer memory’s usage and
larger file size for the COMSOL model to store the data at the nodes.
The more important factor for the transient model is the maximum time stepping size for
the model. As a time dependent COMSOL model, the maximum time stepping size
determines the smoothness and accuracy of the data. Equation 17 can be used to
determine the maximum time stepping size, even though the equation is designed for the
stability in an explicit scheme. Under the implicit method, the time step can be more
lenient while still retaining the data. Using a maximum element size of 0.05x10-6 m and
Equation 17, a time stepping size of 6.944x10-5 s was calculated, resulting in a
recommended time step of less than 1s for the implicit scheme, 0.1s to be exact. If a
lower maximum time step value for the model is used, the model simulation will take
more time to solve. Moreover, the value chosen was set such that the program does not
take longer than expected to run. Based on the results of the other methods and Li et al,
an expected liquid region closure time would be between time t=0 to 1.0x105s.
16
3. RESULTS AND DISCUSSION
In order to evaluate and validate the results obtained in FORTRAN and COMSOL,
results from Li et al [1] must be explored. Following that, results through Ikawa’s
approximation, forward finite difference approximation and COMSOL modeling are
compared, respectively.
3.1 Li’s Method
Figure 4 Liquid Layer Half-Width vs. Bonding Time for Li’s Model [1]
Among the many data plots of results investigated by Li et al, the main result to evaluate
is Figure 4. The findings of Li compared with those of Zhou and North’s experiment
and mathematical model [12], and Illingworth et al’s mathematical model [7]. It shows
that the results published in his paper using fixed-grid finite numerical modeling were
comparable to the results obtained in the models described in [7] and [12]. A main
assumption in all the discussed models above is a constant molar volume. As a point of
reference, Illingworth et al’s method used a variable grid discretization model. Similarly,
Zhou and North’s mathematical model used a fixed-grid discretization model with an
equation describing the motion of the moving interface while solving for the diffusion
equations for the different phases simultaneously.
As shown in Figure 4, a graph of the liquid layer half-width against bonding time, the
liquid layer half width increased between time 0 to 10s prior to decreasing, which can
easily be explained by the concept of super-saturation.
17
In the model, the initial
concentration was Ci=19%, a value far away from its Cl at T=1473K (Cl=10.223%). The
theoretical maximum liquid layer half width, hmax, is 2.32x10-5m [1]. hmax is calculated
from the equilibrium concentration of the solute, 10.233 at.% P, by ignoring any
diffusion of the solute into the solid with the data provided in Li et al. hmax represents the
conservation of the solute.
Zhou and North’s experimental and theoretical exceeded hmax, while Illingworth’s and
Li’s models did not exceed hmax. Illingworth explains the behavior by two possibilities:

Incorrect constant molar volume in the liquid

The fluid flow during the experiments may have affected the thickness of the
liquid layer and diffusion.
with a higher probability of the issue being the incorrect constant molar volume. Li
further explores this issue with a 50% reduction in the molar volume of phosphorous
while keeping the other conditions unchanged, resulting in a solution similar to that of
the experimental noted in Zhou and North. Li concludes that the molar volume of
phosphorous is probably the real reason for the experimental values to be greater than
hmax based on constant molar volume. Constant molar volume will be investigated as a
side study for the project, shown in Appendix D.
As shown in Figure 4, the overall time needed to perform a complete TLP bonding
process, liquid region closure time, is  1.0x105s. This value will be referenced in all the
other model simulation.
3.2 Ikawa’s Method
3.2.1
Program 1
The results provided by the first program evaluate several behaviors within the TLP
bonding process. The first FORTRAN program models the behavior at the solid-liquid
interface. The program evaluates the concentration distribution for diffusion in a semiinfinite medium (base metal) for various time durations, C vs. x. The program’s code,
the program’s input data (Table 1), and output data are summarized in Appendix B,
respectively.
18
Figure 5 Concentration vs. x for Ikawa Program 1
3.2.1.1 Solid –Liquid Layer Interface
Utilizing the initial conditions at t=0s and T=1473K, the concentration of the liquid
superalloy, Ni-P, and the base metal, Ni, is separated at 1.25x10-5m as shown in Figure
5. The behavior of the solid-liquid layer interface matches Figure 2(c) with Cs and Cl
identified by Figure 1 and corresponding with Figure 2(a) using Li’s input parameters.
Cs and Cl play a dominate role of defining the interface as P atoms diffuses into the base
metal while the concentration is limited to Cs at T=1473K. The P atoms move from a
high concentration to a low concentration as defined by diffusion. The concentration of
P remains at 100% from its total initial concentration in the liquid phase. As the system
is exposed for a longer period of time (1 hr to 4 hrs), the solid-liquid layer interface
remains intact. However, the diffusion of the P atoms continues into the solid phase
from the liquid layer. In order to compensate for the diffusion of P atoms, the length of
the liquid layer half width decreases, and the boundary layer (solid liquid layer
interface), moves with the decreasing liquid layer half width. Also, the concentration of
P continues to diffuse further into the base metal with continued temperature exposure
and is displayed in the output data files in Appendix B. The process will continue until
19
the width of the crack is small, similar to the size of the solid-liquid layer interface. At
this point, the melting temperature rises resulting in a decrease of the Ni-P superalloy in
the interlayer. The decrease in the superalloy forces the melting temperature to become
higher than the heat treatment temperature.
Subsequently, the solution solidifies,
completing the TLP bonding process. The overall concentration of P in the system,
regardless of the time position of the TLP bonding process, remains the same. This
behavior explains that the system is a theoretical solution located in a diffusion
controlled environment.
3.2.1.2 Rate of Diffusion in the Liquid Layer Half Width
The final point to notice is the rate of the liquid laer half width change at each value of
time during the process. The liquid layer half width is directly proportional to the
diffusion. The rate of diffusion is implicitly shown to be quicker during the initial
period of time compared to the period towards the end of the TLP bonding process. The
behavior is explained by Ikawa’s formulation, specifically, Equation 6, and is further
analyzed in the second FORTRAN program by Ikawa.
3.2.2
Program 2
The second FORTRAN program models the concentration behavior of the liquid layer
half width at different durations of the TLP bonding process. This program will compile
and displayed results similar to Li et al [1] (i.e. liquid layer half width vs. time). The
program is an extension of the first program to provide a further study of the diffusion
controlled TLP bonding process using Ikawa’s method. The results in the second
program match the results in the first program, and expand on the behavior within the
TLP bonding process. The program’s code and the input (Table 2) and output data are
available in Appendix B, respectively. Other simulations were performed in regards to
different design parameters, i.e. different initial crack width and molar volume values.
They are provided as case studies in Appendix D and will be briefly discussed in the
report.
20
Figure 6 Liquid Layer Half Width vs. Time for h=1.25x10-5 m in Ikawa
Figure 7 Liquid Layer Half Width vs. log Time for h=1.25x10-5 m in Ikawa
21
3.2.2.1 Rate of Diffusion
As shown implicitly and explicitly in Figure 5 and 6 for T=1473K, the rate of the liquid
layer half-width change is at the highest in the beginning durations of the exposure.
Further into the exposure, the rate of the diffusion would slow down until the process
becomes completely diffused (xSL=0).
3.2.2.2 Effects of Varying Temperatures
Varying the temperature had a substantial effect on the duration for completing the TLP
bonding process. An increase in temperature while keeping the other input data
parameters constant, forced the TLP bonding process to complete within a shorter
duration. A decrease in temperature resulted in a longer bonding process time. The
temperature trend was confirmed by the characteristics of diffusivity, Equation 1. It was
expected that if the diffusivity of the material increased from 1.8x10-11m2/s to
5x10-10m2/s, the liquid region closure time shortened. The modeling of the diffusivity is
shown in Appendix D: Case Study 1 using both Ikawa programs.
3.2.2.3 Effects of Changing Initial Liquid Layer Half Width
The liquid region closure time for the TLP process at T=1473K is t=25024s. The time
obtained did not correlate with the completion times provided in Li et al [1]. However,
if the model was designed with an initial parameter h=2.5x10-5m versus h=1.25x10-5m, a
duration for completing the TLP process at T=1473K was t=100094s, keeping all other
initial data parameters constant. The modeling with 2.5x10-5m is shown in Appendix D:
Case Study 2 using the second program.
The only parameter assumption was the value of T=1473K as Li et al did not clearly
state T=1473K was used. A temperature value of T=1383K provides a result of
t=93465s, see Appendix B, which is comparable to the values specified in Li. The effect
of temperatures was described in Section 3.2.2.2, as the program was highly dependent
on the diffusion equation, Equation 1.
22
3.2.2.4 Behavior of the Liquid-Layer Half Width
Ikawa’s model does not take into account the super-saturation levels as in clear form
comparing the results in Figure 4 and 7. According to Figure 4, the liquid layer half
width increased prior to decreasing whereas, in Figure 7, the initial portion of the time
up to 10s was ignored. Ikawa’s theoretical formulation [8] does not factor in the effect of
the super saturation levels. Under Ikawa’s method, the initial data for Cs was not
approached from Ci=19% but rather from Ci=Cl=10.223%.
Also, under Ikawa’s experimental approach [9], there is a sharp peak of P atoms right
before concentration in the bonded interlayer becomes nearly equal to that in the base
metal.
This experimental data confirms isothermal solidification.
However, his
theoretical approach ignores this observation, as the data shows a smooth rate of
diffusion.
3.2.2.5 Effects of Molar Volume Change
In the current project, a constant molar volume was assumed. The modeling of the
effects of molar volume change is shown in Appendix D: Case Study 3 using the second
program. The results demonstrated that the molar volume does not affect the behavior
and completion time of the bonding process.
3.3 Forward Finite Difference
The FORTRAN program modeled the concentration behavior of the liquid layer half
width at different durations of the TLP bonding process using a forward finite difference
under an explicit scheme. The output results were compared to Figure 4 and 7 for
a*=340. This value is explained in Section 3.3.3.
The program’s code, the input (Table 3) and output data are available in Appendix C,
respectively.
23
Figure 8 Liquid Layer Half-Width vs. Time for a*=340 in Forward Finite Difference
Figure 9 Liquid Layer Half-Width vs. log Time for a*=340 in Forward Finite Difference
24
3.3.1
Case Study 1 – Divergence of Cs
It appears that the program did not converge with an input parameter of Cs=0.166%,
Appendix E. However, when an initial design parameter of Cs=0.5% was used, the
results converged. A debugging of the program was performed, but the divergence
continued to exist for very small values of Cs. A resolution of applying Cs=0.5% for the
convenience of the project was performed on all tests.
3.3.2
Behavior of the Liquid Layer Half Width
The forward finite difference model behavior was comparable to the behavior described
in Li et al [1] and in Ikawa’s paper. The results provided an initial super-saturation as
defined by Ci=19%. From that point, the liquid layer half-width jumped from
1.25x10-5m to 2.25x10-5m at the beginning of the process, t=1s. If the step sizes were
smoothed out, similarities appear as shown by comparison of Figure 6 and 7 to 8 and 9,
respectively. The maximum liquid layer half width appeared when t=10s and started
declining at t=100s, similar to the models described in Li et al [1]. However, the
forward finite difference model completed the bonding process (xSL=0m) before Li’s
model. One reason is the use of a forward finite difference versus the method applied in
Li. The second reason is the effect of a*.
3.3.3
Effect of Arbitrary Reference Value of the Activity [a*] and Reciprocals of
Slopes on Binary Phase Diagram [ks]
As explained previously, the input parameters for the arbitrary reference value of the
activity and the reciprocals of the slopes, a* and ks, respectively, were unknown. ks was
dependent on the solidus line of the phase diagram and was a difficult parameter to
accurately interpolate, resulting in an input parameter of ks=-0.0014 through
interpolation of Figure 1. A proper a* was determined using the theoretical liquid layer
half width under constant molar volume of hmax=2.32x10-5m, as mentioned in Li et al
[1]. The results of back-testing a* parameter showed, Figure 10, that a*  340 are
suitable arbitrary reference values assuming ks=-0.0014.
25
Figure 10 Liquid Layer Half-Width vs. log Time for a* in Forward Finite Difference
As a* decreased, the liquid region closure time increased. Similarly, as a* increased, the
liquid region closure time decreased. Ikawa’s liquid region closure time (t=25000s)
corresponds with arbitrary reference value, a*  200 (txa*[email protected]=0=23905s). Similarly,
the maximum liquid layer half width was 2.4x10-5m.
3.4 COMSOL
Through the use of COMSOL, the concentration behavior of the liquid layer half width
at different durations of the bonding process was evaluated. The results obtained in
Figure 4, 7, and 9 were used to compare with the COMSOL model.
The plot of the liquid layer half-width against time was plotted for various  values,
Figure 11.
26
Figure 11 Liquid Layer Half-Width vs. log Time for
3.4.1
  0.006%,0.004%
Effect of Infinitesimal Difference [  ]
 played a very critical role in the COMSOL model.
A large  caused the
concentration slope of the solid and liquid phase to be extremely large as shown from the
behavior expressed in Figure 5. A large difference is not an accurate representation of
the solid-liquid interface behavior. In terms of the user-defined COMSOL functions,
cini(x) and  eff (Ceff) were affected by  . The cini(x) change is not significant to warrant
any actual change in the COMSOL model. However, the effective reciprocal activity
coefficient  eff (Ceff), was largely affected, as it represented the transformation
concentration slope of the solid phase, liquid phase, and the solid-liquid layer interface.
A slight change of  from 0.006% to 0.004%, caused a change in the  eff value at the
solid-liquid layer interface from 838.08 to 1257.13. As shown in Figure 11, this change
has increased the liquid region closure time from 2000s to 6600s. Similarly, with a
smaller  , a finer mesh is required for precision. A coarse mesh will affect the results,
increasing the chance of inaccurate data value during its iterations.
27
3.4.2
Behavior of the Liquid Layer Half Width
The behavior monitored in the COMSOL model is similar to the models described in
Li’s and forward finite difference. The initial concentration, Ceff,i=8.949% or Ci=19%,
demonstrated that the system starts out in super-saturated levels, causing the liquid layer
half-width to increases initially before the effect of diffusion kicks in to decrease the
liquid layer half-width. For the COMSOL model with   0.006% , the total time
needed for the bonding process is t=2000s. For   0.004% , the total time needed was
t= 6600s.
3.4.3
Behavior of the Effective Concentration Values
The C vs x was compared with the results provided in Ikawa. In Figure 5, the behavior
of the solid-liquid layer interface followed the binary phase diagram, Figure 1; more
specifically, Cs and Cl held their solid and liquid concentration separation as the liquid
layer half width decreased. The data provided in Appendix A for the COMSOL models
clearly showed the value of Ceff,l retaining its phase diagram functionality up until a time
period right before the bonding process completes. At that moment, the concentration
quickly dropped to zero, and the bonding process is completed. The value of Ceff,s does
not have the same characteristics as described in Ikawa. At points significantly away
from the interface, the concentration profile spikes, caused by irregularities in the
interpolation data. The error occurs when the chart records an erroneous value in its
interpolation data, i.e. negative value for concentration at any relative x position along
the system. However, the value of Ceff,s never exceeds Ceff,l at values close to the
interface until the process was complete.
3.4.4
COMSOL Limitations
The results obtained for the behavior of the effective concentrations while comparing
with existing literatures questioned the validity of the model. The accuracy of the results
obtained in COMSOL is limited by the computer specifications as accuracy requires a
long processing time and large memory to store large file due to finer mesh sizes
(smaller maximum element size) and smaller time stepping sizes.
28
3.5 Overall Discussion
The details of the results obtained in all the methods were not in correlation with each
other for the liquid region closure times as displayed in Table 7 and Figure 12. Many
discrepancies were discussed above, and all affected the liquid region closure time. The
results of each of the method were plotted together along with the initial data
extrapolated from Li’s graphs. The end results of all the methods starting at Ci=19%
demonstrated an increase in the liquid layer half width before decreasing.
Table 7 Liquid Region Closure Time Summary
Method
Li
Ikawa
Forward Finite Difference with a* = 200
Forward Finite Difference with a* = 340
COMSOL with   0.006%
COMSOL with   0.006%
Time to Complete Bonding, s
10000
25024
23905
9080
2000
6600
Figure 12 Liquid Layer Half Width vs. log Time – All the Methods
29
4. CONCLUSIONS
The project was successful in providing qualitative understanding of the TLP process
and quantitative characterizations of the relationships that exist among the process
parameters, such as concentration distribution in the material system, time and
temperature. All three models indicate that the liquid region ultimately disappears under
isothermal conditions when held for sufficiently long treatment times as P atoms diffuse
out of the region and into the base metal. However, the numerical approaches predicted
initial growth of the liquid region before its disappearance.
While the qualitative
behavior predicted by all three methods was in agreement with prior work, there were
some discrepancies in the predicted liquid region closure times.
The Ikawa’s model did not match Li et al’s model [1] for the initial behavior of the
liquid layer half width and the liquid region closure time. However, the behavior
provided in Ikawa was similar to that of Li’s model. The initial behavior of the liquid
layer half width in Li’s model was not displayed in the Ikawa’s model mainly because
the initial concentration value assumed in Ikawa was Ci=Cs=10.223%, which was the
concentration at the liquid layer interface. Ikawa’s program was designed to match its
theoretical limitations [8]. The initial atomic concentration value of P in the liquid used
in Li et al was Ci=19%, a value indicating super-saturated levels at the T=1473K. Since
the forward finite difference model and the COMSOL model utilized Ci=19% as the
initial concentration value, the results demonstrated a similar initial behavior of the
liquid layer half-width in Li’s model.
In the forward finite difference model, the liquid layer half-width quickly jumped up to
its super-saturation level at its corresponding temperature, T=1473K.
As the
concentration of P atom diffused from the liquid layer to the solid layer, the movement
was similar to that provided in Ikawa and Li. Like Ikawa’s method, the forward finite
difference model did not match Li’s model in liquid region closure times.
The
inconsistent quantitative results were due to the initial input parameters for Cs, a* and ks.
The forward finite difference program diverged if a value of Cs=0.16% was inputted
resulting in the use of Cs=0.5%. Since a* and ks was unknown, an approximation of the
solidus line was obtained using the Figure 1, resulting in a ks=-0.0014. This input
parameter was used to back-test the corresponding arbitrary reference value, a*. The
30
three unknown parameters were factors in changing the liquid region closure times for
the forward finite difference program. The value of a* was determined to be, 340 and
200, for hmax=2.32x10-5m and a liquid region closure time of t=25000s, respectively. As
we decrease a*, we can obtain a value similar to the results provided in previously
presented results.
The Ikawa’s model displayed a clear characterization of the bonding process by showing
the behavior of the solid-liquid layer interface as it closely resembles the properties of
the binary phase diagram for Ni-P, Figure 1. The interface remained in-tact until the
liquid layer half width becomes that of the interface width, infinitesimal small value.
From there, the melting point rises causing the remaining liquid layer to solidify,
representing the completion of the TLP bonding process.
Several side studies were performed for the Ikawa model. If the initial half width of
crack was increased to 2.5x10-5m (double), versus the 1.25x10-5m, the results obtained
were comparable to the results provided in Li. Likewise, the Ikawa model also
demonstrated that the assumption of a constant molar volume has very little bearing on
the liquid region closure time resulting hardly in any change in the overall bonding
process.
The COMSOL model was the most difficult model to use in analyzing the TLP process,
as there were no existing references to demonstrate all of COMSOL’s capabilities. The
model in COMSOL completed the liquid region closure much quicker than any of the
models described before. It exhibited a higher liquid layer half width, far more than its
recommended theoretical value [1].
Additionally, the results provided by the model
were limited to the capabilities of the computer. A decrease in  will result in a
longer liquid region closure times, a result that would be desired. However, for the
decrease of  , the calculation was limited by the capabilities of software; i.e. element
and time spacing. Even though runs were performed with very thin step sizes, the data
did not populate as expected from the other models. Secondly, the solid-liquid layer
interface behavior in the COMSOL model did not match the behavior represented in
Ikawa’s first program. Although the Ceff,l held its phase diagram properties until the
completion of the bonding process, Ceff,s and the concentration values far away from the
crack in the solid layer did not yield similar behavior as displayed in Ikawa’s program.
31
Although the post-processing from the COMSOL model did not match in exactly the
results of the previous methods, the overall qualitative and quantitative characteristics of
the TLP bonding process were successful. Many diffusion-controlled behaviors were
apparent in the Ikawa model and forward finite difference. Referencing the models of
Ikawa [8,9], forward finite difference and COMSOL to Li et al was difficult, as the
initial input parameters were limited to those that were provided in their articles [1] and
any other references, [7, 10, 12]. These issues aided in the inaccuracies of the results.
4.1 Future Work
While the COMSOL model qualitatively predicted the trend, it did not match the exact
results of the previous models, with an unknown resolution being the issue. The model
needs additional investigation of its diffusion application model under transient analysis.
It is unclear whether TLP bonding process can be modeled in COMSOL accurately as
there are no existing literatures to verify its capabilities.
32
5. REFERENCES
1. Li, J. F.; Agyakwa, P. A.; Johnson, C. M., A fixed-grid numerical modeling of
transient liquid phase bonding and other diffusion-controlled phase changes, Springer
Science & Business Media LLC., J Mater Sci (2010) 45:2340-2350
2. Arafin, M. A.; Medraj, M.; Turner, D. P.; Bocher P., Transient liquid phase bonding
of Inconel 718 and Inconel 625 with BNi-2: Modeling and experimental investigations,
Elsevier B.V., Materials Science and Engineering A 447 (2007) 125-133
3. Pouranvari, M.; Ekrami, A.; Kokabi, A. H., Effect of bonding temperature on
microstructure development during TLP bonding of a nickel base superalloy, Elsevier
B.V., Journal of Alloys and Compounds 469 (2009) 270-275
4. Lo, P.-L.; Chang, L.-S.; Lu Y.-F., High strength alumina joints via transient liquid
phase bonding, Elsevier ltd and Techna Group S.r.l., Ceramics International 35 (2009)
3091-3095
5. Wikstrom, N. P.; Egbewande, A. T.; Ojo, O. A., High temperature diffusion induced
liquid phase joining of a heat resistant alloy, Elsevier B.V, Journal of Alloys and
Compounds 460 (2008) 379-385
6. AlHazaa, A.; Khan, T. I.; Haq, I., Transient liquid phase (TLP) bonding of AL7075 to
Ti-6Al-4V alloy, Elsevier Inc., Materials Characterization 61 (2010) 312-317
7. Illingworth, T. C.; Golosnoy, I. O.; Clyne, T. W., Modeling of transient liquid phase
bonding in binary systems-A new parametric study, Elsevier B. V, Materials Science and
Engineering A 445-446 (2007) 493-500
8. Ikawa, H.; Nakao, Y.; Isai, T., Theoretical Considerations on the Metallurgical
Process in T.L.P. Bonding of Nickel-Base Superalloys, Trans. of the Japan Welding
Society, Vol. 10, No. 1 (1979) 24-29
9. Ikawa, H.; Nakao, Y., Transient Liquid Phase Diffusion Bonding of Nickel-Base Heat
Resisting Alloys, Trans. of the Japan Welding Society, Vol. 8, No. 1 (1977) 3-8
10. Journal of Phase Equilibrium Volume 21, Number 2, 210, DOI
10.1361/105497100770340345
11. Crowley, A. B. and Ockendon, J. R., On the Numerical Solution of an Alloy
Solidification Program, HMT Vol. 22, No 6-K (1978)
12. Zhou, Y. and North, T. H., Kinetic modeling of diffusion-controlled, two-phase
moving interface problems, Model Simulation Material Science (1993) 1(4):505
33
13. Wikipedia contributors. "COMSOL Multiphysics." Wikipedia, The Free
Encyclopedia. Wikipedia, The Free Encyclopedia, 16 Nov. 2010. Web. 5 Dec. 2010.
14. COMSOL Multiphysics Help Desk Documentation, Modeling Guide [modeling.pdf]
34
6. APPENDIX A – COMSOL Reference
This appendix will cover all the information needed for COMSOL and will
supply/reference the necessary data.
6.1 Developing the Model: Initial Concentration [cini(x)], Effective
Reciprocal Activity Coefficient [  eff (Ceff)]
In developing the COMSOL model, the initial concentration of the model and the solidliquid layer interface are two important characteristics that need to be addressed. The
actual concentrations are known factors, as mentioned in Li et al [1] and used in the
other models. The separation of the solid and liquid layer is initially at 12.5x10-5m. This
value corresponds with the values Cs=0.166% and Cl=10.223%. In order to create a
model in COMSOL to incorporate both values, an infinitesimal difference, 2 in the
concentration is created to mimic an approximation. The center of the difference would
be the Cs=0.166% value and the value for  is determined.
Under the ideal case, the Ceff correlates with C at the solid and liquid phases in a direct
one-to-one linear relationship. Applying the above conditions Ceff,b, Ceff,s, Ceff,l, and Ceff,i
are established for Cb, Cs, Cl, and Ci and are provided in Table 8 for the case of
  0.006% and   0.004% with Figure 13 and 14 to demonstrate the linear
distribution of   0.006% and   0.004% respectively.
Table 8 Initial Data: Effective and Actual Concentration for Various
Positions, x, m
2.5x10-4
1.25x10-5
1.25x10-5
0

Ceff
0%
0.16%
0.172%
8.949%
 eff
6E-05
C
0%
0.166%
10.223%
19%
838.0833
35

Ceff
0%
0.162%
0.17%
8.947%
 eff

4E-05
C
0%
0.166%
10.223%
19%
1257.125
Figure 13 Actual vs. Effective Concentration for   0.006%
Figure 14 Actual vs. Effective Concentration for   0.004%
6.1.1
Sample Calculation for Effective Concentration [Ceff,b] and Effective
Reciprocal Activity Coefficient [  eff ]
Ceff,b=8.949% was determined using an effective concentration transformation assuming
an ideal case s=1 (one-to-one relationship)
s
0.19  0.10223
 1  x  0.08949
x  0.0172
Similarly,  eff the effective reciprocal activity coefficient, is formulated by
 eff 
Cl  C s
0.10223  0.00166

 838.08 which represent the slope of the
Ceff ,l  Ceff , s 0.00172  0.00166
concentration.
36
6.2 COMSOL Data
6.2.1
6.2.1.1
Models
 = 0.004%:
d0.004%-overall.mph
6.2.1.2
 = 0.006%:
d0.006%-overall.mph
d0.006%-records1s.mph
6.2.2
6.2.2.1
Output Data
 = 0.004%:
Folder: d0.004%
6.2.2.2
 = 0.006%:
Folder: d0.006%
6.2.2.3 Liquid Layer Half Width vs. xSL
Folder: Summary Data
37
7. APPENDIX B – FORTRAN: Ikawa Method
7.1 First Program
7.1.1
C
C
C
C
C
C
C
C
C
C
Code
ikawa.f
PROGRAM FOR THE CALCULATION OF CONCENTRATION
DISTRIBUTION FOR DIFFUSION IN A SEMI-INFINITE
MEDIUM FROM AN EXTENDED SOURCE OF LIMITED EXTENT.
THE MOVING BOUNDARY MOTION IS APPROXIMATELY
TAKEN INTO ACCOUNT ACCORDING TO IKAWA'S METHOD
PARAMETER (NX=101)
PARAMETER (NT=5)
C
C
NEXT PARAMETER REQUIRED FOR PLOTTING !!!
PARAMETER (MX=103)
C
C
DIMENSION C(NX,NT),X(MX),X2(NX),T(NT),HWC(NT)
DIMENSION CO(MX),C1(MX),C2(MX),C3(MX),C4(MX)
DIMENSION XXP(2),YYP(2)
C
C
OPEN (55,FILE='ikawa_revised.dat')
OPEN (66,FILE='ikawa.res')
C
C
C
READ INPUT DATA
READ(55,*) CL,CS,HWCO,D,VOLM,DXDIST,DTIME,
1 EPSI
C
C
C
C
C
C
C
C
C
C
C
C
INTRODUCE AND COMPUTE PRELIMINARY PARAMETERS
PI=3.141592654
A1=0.0705230784
A2=0.0422820123
A3=0.0092705272
A4=0.0001520143
A5=0.0002765672
A6=0.0000430638
RGAS=8.3144
OPE=1.0 + EPSI
OME=1.0 - EPSI
D=DO*EXP(- Q/(RGAS*TABS))
write(6,*) D
XTOTAL=DXDIST*FLOAT(NX-1)
TTOTAL=DTIME*FLOAT(NT-1)
XXP(1)=HWCO
XXP(2)=HWCO
YYP(1)=0.0
YYP(2)=1.0
NEXT PARAMETERS REQUIRED FOR PLOTTING !!!
N1=NX+1
N2=NX+2
COMPUTE THE CONCENTRATION PROFILES FOR VARIOUS TIMES
START TIME MARCHING
38
C
TIME=0.0
T(1)=TIME
HWC(1)=HWCO
DO 500 IT=2,NT
TIME=TIME + DTIME
T(IT)=TIME
XCAR=2.0*SQRT(D*T(IT))
FACT1=2.0*CS/(SQRT(PI)*VOLM)
FACT2=(CL/VOLM - CS/VOLM)
HWC(IT)=HWCO - (FACT1/FACT2)*XCAR/2.0
IF(HWC(IT).LT.0.0) HWC(IT)=0.0
XXP(1)=HWC(IT)
XXP(2)=HWC(IT)
XDIST=0.0
DO 100 IX=1,NX
X(IX)=XDIST
IF(X(IX).LT.OME*HWC(IT)) C(IX,IT)=CL
IF(X(IX).LT.OME*HWC(IT)) GO TO 75
25 CONTINUE
IF(X(IX).GE.OME*HWC(IT).AND.X(IX).LE.OPE*HWC(IT))
1 C(IX,IT)=CS
IF(X(IX).GE.OME*HWC(IT).AND.X(IX).LE.OPE*HWC(IT))
1 GO TO 75
50 CONTINUE
IF(X(IX).GT.OPE*HWC(IT)) CONTINUE
X2(IX)=X(IX) - HWC(IT)
XX2=X2(IX)/XCAR
ATX2=A1*XX2
ATX2=ATX2+A2*XX2*XX2
ATX2=ATX2+A3*XX2*XX2*XX2
ATX2=ATX2+A4*XX2*XX2*XX2*XX2
ATX2=ATX2+A5*XX2*XX2*XX2*XX2*XX2
ATX2=ATX2+A6*XX2*XX2*XX2*XX2*XX2*XX2
C
ERF2=1.0-1.0/(1.0+ATX2)**16
C
C(IX,IT)=CS*(1.0 - ERF2)
C
75 CONTINUE
XDIST=XDIST + DXDIST
100 CONTINUE
500 CONTINUE
C
C
DO 600 IX=1,NX
C
C0(IX)=C(IX,1)/CL
C1(IX)=C(IX,2)/CL
C2(IX)=C(IX,3)/CL
C3(IX)=C(IX,4)/CL
C4(IX)=C(IX,5)/CL
C
write(55,*) float(IX-1)*DXDIST,C0(IX)
write(16,*) float(IX-1)*DXDIST,C1(IX)
write(17,*) float(IX-1)*DXDIST,C2(IX)
write(18,*) float(IX-1)*DXDIST,C3(IX)
write(19,*) float(IX-1)*DXDIST,C4(IX)
600 CONTINUE
C
c
write(66,990) (X(IX),IX=1,NX)
c
write(66,990) (C1(IX),IX=1,NX)
c 990 format(1x,5(e12.6,2x))
C
C
C
700 CONTINUE
C
C
STOP
END
39
7.1.2
Input Data
10.223 0.166 12.5e-6 1.8e-11 7.344e-6 1e-6 3600.0 0.001
7.1.3
Output Data
Ikawa – First Program/DATA-12.5micro-meters/initial
Ikawa – First Program/DATA-12.5micro-meters/1hr
Ikawa – First Program/DATA-12.5micro-meters/2hrs
Ikawa – First Program/DATA-12.5micro-meters/3hrs
Ikawa – First Program/DATA-12.5micro-meters/4hrs
7.2 Second Program
7.2.1
Code
C
C
C
C
C
C
C
ikawa1.f
PROGRAM FOR THE CALCULATION OF THE HALF WIDTH
OF CRACK BEING REPAIRED (IKAWA'S METHOD)
PARAMETER (NT=400001)
PARAMETER (NTEMP=11)
C
C
C
NEXT PARAMETER REQUIRED FOR PLOTTING !!!
PARAMETER (MT=NT+2)
C
C
DIMENSION T(MT),HWC(NT,NTEMP),TABS(NTEMP)
DIMENSION C1(MT),C2(MT),C3(MT),C4(MT),C5(MT),
1 C6(MT),C7(MT),C8(MT),C9(MT),C10(MT),C11(MT)
C
C
OPEN (55,FILE='ikawa1_revised.dat')
C
C
C
READ INPUT DATA
READ(55,*) CL,CS,HWCO,DO,Q,VOLM,
1 TLOW,THIGH,DTIME
C
C
C
INTRODUCE AND COMPUTE PRELIMINARY PARAMETERS
PI=3.141592654
RGAS=8.3144
DTEMP=(THIGH-TLOW)/FLOAT(NTEMP-1)
TTOTAL=DTIME*FLOAT(NT-1)
C
C
C
C
CC
C
C
NEXT PARAMETERS REQUIRED FOR PLOTTING !!!
N1=NT+1
N2=NT+2
FOR EACH TEMPERATURE COMPUTE THE HALF WIDTH OF THE CRACK
40
C
C
TEMP=TLOW
DO 200 ITEMP=1,NTEMP
TABS(ITEMP)=TEMP
D=DO*EXP(- Q/(RGAS*TABS(ITEMP)))
write(6,*) TABS(ITEMP), D
TIME=0.0
T(1)=TIME
HWC(1,ITEMP)=HWCO
DO 100 IT=2,NT
TIME=TIME + DTIME
T(IT)=TIME
XCAR=2.0*SQRT(D*T(IT))
FACT1=2.0*CS/(SQRT(PI)*VOLM)
FACT2=(CL/VOLM - CS/VOLM)
HWC(IT,ITEMP)=HWCO - (FACT1/FACT2)*XCAR/2.0
IF(HWC(IT,ITEMP).LT.0.0) HWC(IT,ITEMP)=0.0
100
CONTINUE
TEMP=TEMP + DTEMP
200 CONTINUE
C
C
DO 600 IT=1,NT
C1(IT)=HWC(IT,1)
C2(IT)=HWC(IT,2)
C3(IT)=HWC(IT,3)
C4(IT)=HWC(IT,4)
C5(IT)=HWC(IT,5)
C6(IT)=HWC(IT,6)
C7(IT)=HWC(IT,7)
C8(IT)=HWC(IT,8)
C9(IT)=HWC(IT,9)
C10(IT)=HWC(IT,10)
C11(IT)=HWC(IT,11)
write(26,*) float(IT-1)*DTIME,C1(IT)
write(27,*) float(IT-1)*DTIME,C6(IT)
write(28,*) float(IT-1)*DTIME,C11(IT)
600 CONTINUE
C
C
700 CONTINUE
C
C
STOP
END
7.2.2
Input Data
10.223 0.166 12.5e-6 0.01121 248000 17.344e-6
1423.0 1523.0 1
7.2.3
Output Data
Ikawa – First Program/DATA-12.5micro-meters/1383K
Ikawa – First Program/DATA-12.5micro-meters/1423K
Ikawa – First Program/DATA-12.5micro-meters/1473K
Ikawa – First Program/DATA-12.5micro-meters/1523K
41
8. APPENDIX C – FORTRAN: Forward Finite Difference Method
8.1 Code
C
C
C
C
C
C
C
C
C
C
C
tlp.f
PROGRAM FOR CALCULATION OF MOVING BOUNDARY MOTION AND
CONCENTRATION PROFILES DURING TRANSIENT LIQUID PHASE (TLP)
BONDING
INPUT DATA : DL,DSS,EKAL,EKAS,CLIQ,CSOL,
CI,CB,ASTAR,XIL,XTOTAL
INPUT GRID : NX,NT,KOUNTP
PARAMETER(NX=2001)
C
C
C
C
C
DIMENSION VARIABLES
DIMENSION AO(NX),CO(NX),UO(NX),A(NX),C(NX),U(NX)
DIMENSION FS(NX)
DIMENSION XX(NX),C1(NX)
C
C
C
C
READ INPUT DATA
OPEN(55,FILE='tlp.dat')
OPEN(66,FILE='tlp.res')
C
READ(55,*)DL,DS,EKAL,EKAS,CLIQ,CSOL,CI,CB,ASTAR,XIL,
1 XTOTAL,TTOTAL,KOUNTP
C
C
C
C
COMPUTE AUXILIARY QUANTITIES
DPL=- DL*EKAL
DPS=- DS*EKAS
DPALS=(DPL + DPS)/2.0
C
NXM1=NX-1
C
C
C
C
C
C
C
C
NEXT PARAMETERS REQUIRED FOR PLOTTING !!!
N1=NX+1
N2=NX+2
COMPUTE MINIMUM DT REQUIRED FOR STABILITY
NXM1=NX-1
DX=XTOTAL/FLOAT(NXM1)
NXIL=INT(XIL/DX) + 1
DX2=DX*DX
DTIMEL=DX2/(2.1*DL)
DTIMES=DX2/(2.1*DS)
DTCFL=AMIN1(DTIMEL,DTIMES)
dt1=1.0
DT=amin1(dtcfl,dt1)
NT=INT(TTOTAL/DT)+1
NNXX=(NT-1)/KOUNTP + 1
C
C
C
INITIALIZE VARIABLES
TIME=0.0
DXSL=0.0
42
X=0.0
DO 100 IX=1,NX
XX(IX)=X
IF(X.LT.XIL) AO(IX)=-CI/EKAL
IF(X.GE.XIL) AO(IX)=-CB/EKAS
IF(X.LT.XIL) UO(IX)=DPL*(AO(IX) - ASTAR)
IF(X.GE.XIL) UO(IX)=DPS*(AO(IX) - ASTAR)
IF(X.LT.XIL) CO(IX)=CI
IF(X.GE.XIL) CO(IX)=CB
X= X + DX
100 CONTINUE
C
do ix=1,nx
write(67,*) DX*float(ix-1),CO(ix)
enddo
C
write(6,*) time,1.0,nxil,xil
c
write(6,*) time,dxsl
c
write(61,990) nnxx
c
write(61,991) time,dxsl
C
C
WRITE(6,1000) TIME
C
WRITE(6,2000) (CO(IX),IX=1,200)
C
C
ADVANCE THE SOLUTION ONE TIME STEP
C
TIME=0.0
KOUNT=0
DO 500 IT=1,NT
TIME= TIME + DT
KOUNT= KOUNT + 1
DO 250 IX=2,NXM1
IXP1=IX+1
IXM1=IX-1
D2UDX2= (UO(IXP1)-2.0*UO(IX)+UO(IXM1))/(DX*DX)
C(IX)= CO(IX) + DT*D2UDX2
IF(C(IX).GT.CLIQ) U(IX)=(-C(IX)/EKAL - ASTAR)*DPL
IF(C(IX).LE.CLIQ.AND.C(IX).GE.CSOL) U(IX)=0.0
IF(C(IX).LT.CSOL) U(IX)=(-C(IX)/EKAS - ASTAR)*DPS
IF(U(IX).GE.0.0) A(IX)=U(IX)/DPL + ASTAR
IF(U(IX).LE.0.0) A(IX)=U(IX)/DPS + ASTAR
250
CONTINUE
C
C
BOUNDARY CONDITIONS
C
(AT X=0.0)
C
IF(CO(1).GT.CLIQ.AND.C(2).GT.CLIQ) DP12=DPL
IF(CO(1).LE.CLIQ.AND.C(2).GE.CLIQ) DP12=DPALS
IF(CO(1).LE.CSOL.AND.C(2).GE.CLIQ) DP12=DPALS
IF(CO(1).LE.CSOL.AND.C(2).GT.CSOL) DP12=DPALS
IF(CO(1).LT.CSOL.AND.C(2).LT.CSOL) DP12=DPS
c
C(1)=CO(1) + 2.0*DT*(UO(2) - UO(1))/DX2
C(1)=C(2)
IF(C(1).GT.CLIQ) U(1)=(-C(1)/EKAL - ASTAR)*DPL
IF(C(1).LE.CLIQ.AND.C(1).GE.CSOL) U(1)=0.0
IF(C(1).LT.CSOL) U(1)=(-C(1)/EKAS - ASTAR)*DPS
IF(U(1).GE.0.0) A(1)=U(1)/DPL + ASTAR
IF(U(1).LE.0.0) A(1)=U(1)/DPS + ASTAR
C
C
(AT X=XTOTAL)
C
A(NX)=A(NXM1)
U(NX)=U(NXM1)
C(NX)=C(NXM1)
C
C
PRINT RESULTS
C
C
IF(KOUNT.GE.KOUNTP) WRITE(6,1000) TIME
C
IF(KOUNT.GE.KOUNTP) WRITE(6,2000) (C(IX),IX=1,200)
C
IF(KOUNT.GE.KOUNTP) WRITE(6,2000) (A(IX),IX=1,NX)
43
C
C
C
IF(KOUNT.GE.KOUNTP) WRITE(6,2000) (FS(IX),IX=1,NX)
IF(KOUNT.GE.KOUNTP) WRITE(6,2000) (U(IX),IX=1,NX)
IF(KOUNT.NE.1) GO TO 350
300 CONTINUE
C
C
TRACK DOWN RATE OF DISPLACEMENT OF CRACK WALL
C
DO 325 IX=1,NXM1
IXP1=IX+1
C1(IX)=CO(IX)/CI
C1(IXP1)=CO(IXP1)/CI
IF(C1(IX).GT.CSOL/CI.AND.C1(IXP1).LE.CSOL/CI) IXSL=IX
IF(C1(1).LE.CSOL/CI) IXSL=1
IF(IT.EQ.1) IXSL=NXIL
DXSL=(XX(IXSL) - XIL)*1e6
XSLL=XX(IXSL)*1e6
IF(IT.EQ.1) DXSL=0.0
325 CONTINUE
C1(NX)=CO(NX)/CI
C
write(6,*) time,c1(1),ixsl,xx(ixsl)
C
write(66,991) time,dxsl
write(86,*) time,xsll*1e-6
990 format(1x,i10)
991 format(1x,2(f12.6,2x))
C
C
C
350 CONTINUE
IF(KOUNT.GE.KOUNTP) KOUNT=0
C
C
C
UPDATE VARIABLES
C
DO 400 IX=1,NX
IF(C(IX).GT.CLIQ) FS(IX)=0.0
IF(C(IX).LT.CSOL) FS(IX)=1.0
IF(C(IX).LE.CLIQ.AND.C(IX).GE.CSOL)
1
FS(IX)=(CLIQ-C(IX))/(CLIQ-CSOL)
AO(IX)=A(IX)
UO(IX)=U(IX)
CO(IX)=C(IX)
400
CONTINUE
500
C
CONTINUE
do ix=1,nx
write(68,*) DX*float(ix-1),C(ix)
enddo
C
1000
2000
2500
FORMAT(1X,F20.9)
FORMAT(1X,5F14.6)
FORMAT(1X,5E14.6)
STOP
END
8.2 Input Data
5e-10 1.8e-11 -0.02684 -0.0014 10.223 0.5 19.0 0.0 340
0.0000125 0.003 20000.0 60
8.3 Output Data
Forward Finite Difference/Case Studies/ a_star-160
44
Forward Finite Difference/Case Studies/ a_star-200
Forward Finite Difference/Case Studies/ a_star-250
Forward Finite Difference/Case Studies/ a_star-340
Forward Finite Difference/Case Studies/ a_star-353.6
45
9. APPENDIX D – Ikawa Case Studies
9.1 Case Study 1: Effect of Diffusivity
The results of concentration vs x are compared for temperature T=1473K at t=60s, for
D=1.8x10-11m2/s and D=5.0x10-10m2/s using both Ikawa programs. The results are
plotted in Figure 15 and 16.
Figure 15 Concentration vs. x for Varying Diffusivity
An increase in the diffusivity causes the processing to complete quicker. A diffusivity of
D=5.0x10-10m2/s had a completion time of t=900s. The data files are provided below for
Figure 15 and 16.
9.1.1
Data Files
Ikawa - First Program/Case Study 1/Program 1/ [email protected]
Ikawa - First Program/Case Study 1/Program 1/ [email protected]
Ikawa - First Program/Case Study 1/Program 2/ [email protected]
Ikawa - First Program/Case Study 1/Program 2/ [email protected]
46
Figure 16 Liquid Layer Half-Width vs. Time for Varying Diffusivity
9.2 Case Study 2: Liquid Layer Half Width Change [h=2.5x10-5 m]
Figure 17 Liquid Layer Half-Width vs. Time for h=2.5x10-5 m in Ikawa
47
Figure 18 Liquid Layer Half-Width vs. log Time for h=2.5x10-5 m in Ikawa
The results obtained in this Case Study shows that at T=1473K for h=2.5x10-5m, the
total time for completing the TLP bonding process is comparable to that of Figure 4. The
data files are provided below for Figure 17 and 18.
9.2.1
Data Files
Ikawa – Second Program/DATA-25micro-meters/T1423K
Ikawa - Second Program/DATA-25micro-meters/T1473K
Ikawa - Second Program/DATA-25micro-meters/T1523K
Ikawa - Second Program/DATA-25micro-meters/6hrs
Ikawa - Second Program/DATA-25micro-meters/12hrs
Ikawa - Second Program/DATA-25micro-meters/18hrs
Ikawa - Second Program/DATA-25micro-meters/24hrs
9.3 Case Study 3: Molar Volume Change
The results obtained in this Case Study shows that the molar volume does not have a
significant effect on the duration for complete TLP bonding process. The values used for
the molar volumes were random and varied in orders of magnitude, but as shown by the
48
data provided below in Figure 19, it clearly shows that molar volume is not a huge factor
in altering the results.
Figure 19 Liquid Layer Half-Width vs. log Time for Varying Molar Volume at h=12.5x10-6 m in
Ikawa
9.3.1
Data Files
Ikawa – Second Program/Case Study - Molar Volume/6.334e-6
Ikawa – Second Program/Case Study - Molar Volume /7.334e-6
Ikawa – Second Program/Case Study - Molar Volume /12.334e-6
Ikawa – Second Program/Case Study - Molar Volume /17.334e-6
49
10.APPENDIX E – Forward Finite Difference Case Studies
10.1 Case Study 3: Concentration of Solid Input Parameter [Cs=0.16%]
The results obtained for Cs=0.16% in Figure 20 demonstrated that there are some errors
to the program, as the liquid layer half-width did not converge. The cause of the error
can be the result of having a very small Cs. The program was debugged but a resolution
using Cs=0.16% was improbable.
Figure 20 Liquid Layer Half-Width vs. log Time for Cs=0.16% in Forward Finite Difference
10.1.1 Data Files
Forward Finite Difference/Case Studies/Cs0.0016
50
Fly UP