Numerical Method
In our model, we are implementing two concomitant processes, diffusion and migration on the surface.
Let P be the free energy reduction associated with per unit volume of matter added to per unit surface area of the solid, and j be the volume of matter added to per unit area of the solid surface per unit time. The specific evaporation-condensation rate, L, is a migration kinetic parameter.
Let F be the driving force and J be the volumetric diffusion flux. They are proportional to each other. The diffusion rate, M, is diffusion kinetic parameter.
Let 
be the normal velocity of the surface. If diffusion is the only kinetics, we have
due to mass conservation. For concurrent interface migration and diffusion, mass conservation gives
Let 
be the volume of matter added to per unit area of the solid surface, and be the mass displacement. The two kinds of virtual changes together cause the surface to move in the normal direction by distance
Associated with the virtual motion, the free energy changes by 
. According to the definition of the two kinds of driving forces, F and P, matter relocation and exchange on the surface reduce the free energy, namely
In formulating a finite element method, it is convenient to use 
and 
as the basic variables. Eliminating F and P in the above equation, we obtain
By using finite element method, we can describe forces action on nodes by an element.
Using interpolation on an element, each element has four degrees of freedom 
to describe the motion of the element, and three degrees of freedom 
to describe diffusion on the element.
One can first construct a global H matrix with all zero elements. Then count all the elements, calculate each 7×7 local H matrix, and add them to the correct locations in the global H matrix.
Global H matrix
Property calculation
-Variables as a function of temperature
In our model, we want to observe the surface morphology change due to annealing temperature. The diffusion mobility M and the migration mobility L is also function of temperature, and it was calculated as follows. It is assumed that the mobility is isotropic and irepresented by a number related to the self-diffusivity by the Einstein relation. Au (gold) deposited thin film is used for our simulation
[1]Arcidiacono, S., Bieri, N. R., Poulikakos, D., and Grigoropoulos, C. P., "On the Coalescence of Gold Nanoparticles," International Journal of Multiphase Flow, 30 7-8 (2004) 979-994.
[2] Egry, I. , Lohoefer, G., Jacobs, G., “Surface tension of liquid metals: results from measurements on ground and in space” Physical Review Letters, 75 (1995) 4043-4046
Temperature profile
Annealing time was 30 minutes including 5 minute temperature increase and 5 minute temperature decrease. Aneealing temperature was 300K, 350K, 400K and 459K.
Initial surface geometry
Surface roughness quantification
Groove depth and surface roughness of the gold thin film was calculated using Root Mean Square.
Nodal point update
Step1: new nodal points are calculated using Euler method
Step2: nodal points are updated using linear interpolation to resolve the problem induced by irregular element size
