Thermodynamics: Chemistry Optimization

The following is an excerpt from Dewar - Characterization and Evaluation of Aged 20Cr32Ni1Nb Stainless Steels, and references a stainless steel used in this thesis. The content of this article can also be applied to a general case, and should not be limited to to the presented example.

1. Batch ThermoCalc Executables

In the present study we are interested in seeing how nitrogen, and titanium affect the resulting microstructure, as well as their interactions with the other components in a 2032Nb alloy. From the ASTM specification for 20Cr-32Nb-Nb the compositional matrix in Table 1 can be used see the ranges for each element that will be explored. In ThermoCalc the stepping command will be used over the temperature range Tm - Trt for finite compositions of the alloy. Therefore, certain levels in the compositional ranges for each element will need to be decided where then equilibrium will be calculated for all possible combinations of these finite chemistries. If a maximum, median, and minimum values are chosen for each range than ThermoCalc will need to run stepping procedures for 37 = 2187 chemical compositions. Writing these ThermoCalc files out manually would be quite time consuming, so two batch programs were written in Visual C# to automatically write and run numerous ThermoCalc log files.

Table 1: Range of 20-32Nb Stainless Steel compositions investigated in the ThermoCalc test matrix. Compositional ranges correspond to ASTM A351, grade CT15C.
Ni Cr Nb Si Mn C N/Ti W,Mo,Ti,Zr
Compositional Range
31-34 19-21 0.5-1.5 0.5-1.5 0.15-1.5 0.05-0.15 0-0.3 < 0.05

The first program titled BatchTCC_Comp shown in Figure 1 takes the compositional range specified, increments each component by the # of iterations set by the user, and then compiles .tcm files for each chemistry combination and runs the scripts. The y-axis is specified by the user where depending on their choice either a stepping or mapping procedure will take place. For this program the x-axis is suppose to be a component, but other properties like temperature and pressure can be used as well. The program outputs both the postscript, and the excel files for the phase or property diagram, and also graphs the composition for each stable phase in the calculation. For this script the relevant phases must be specified, which can be determined by manually doing a test run of the system through ThermoCalc beforehand. The second program titled BatchTCC_Temp is a much simpler version of the first program where only the temperature can be specified as the x-axis for a stepping procedure. If all of the alloy element boxes are used, and the # of iterations is limited to three, the process should be expected to take at least 30 hrs if the ‘# of simultaneous processes’ is set to 3-5. Both the executable installations, and the source code Visual Studio solutions can be downloaded from

Fig. 1: GUI interface for batch running a compositional matrix through ThermoCalc where the x-axis is a component of the system, and the y-axis is either the phase fraction, another component of the system, or temperature.

Fig. 2: GUI interface for batch running a compositional matrix through ThermoCalc where the x-axis is a temperature range.

1.1. Preparing for Compilation

After the process has been completed it is important to copy all of the excel file outputs into a separate folder, as some of the subsequent scripts can not handle searching for files through subfolders. Secondly, Thermocalc outputs the excel files in .xls format which conflicts with some of the Matlab code that can only open .xlsx formatted excel files. The first step is to go into the top folder of the directory where the batch ThermoCalc process saved the files. The batch process program saves all of the steping, or mapping plots in the format “< element >< composition >...”, where for a 20Cr32Ni alloy the file name would be “Cr20Ni32...”. For Windows 7 operating systems, find the search bar in the file explorer, and type “Cr*.xls”. This will find all the excel files that start with the characters ‘Cr’. Copy and save all of these files to another folder. Next, download the files named xls2xlsx.rar from Extract, and open the excel file. In the first section type in the path of the copied .xls files in the “Original File Path” cell, and make a new folder for the .xlsx files, and type its path in the “Destination File Path” cell. Click the “Convert to xlsx” button, and all the .xls file should be saved as .xlsx files in the new folder. The same can be done for any compositional data for a specific phase by specifying the cells under the “Compostion xls File” title.

2. Matlab Code - Compiling Mass ThermoCalc Data To Excel

Once all of the ‘.xlx’ files have been converted to ‘.xlsx’ format, and are copied to a single file (eg. ‘_xlsx’), the Matlab distribution files for the batch ThermoCalc application can be employed to extract, and compile all of the relevant information needed to produce matrix arrays, and linear regression fitting. From the volume fraction plots, there are specific points on the curves for each phase that can be useful when trying to quantify the optimization of an alloy based on equilibrium microstructure.

2.1. Finding the solubility temperature, and the terminal phase fractions for intermetallics and Chrome Carbides

The first Matlab script will be used for compiling all the relevant information about any intermetallic phases, or chrome carbides that may have precipitated during long term aging. For example, the G-phase curve in  Figure 3a has three definable points that can be analyzed for all the chemistries when determining the minimization of G-phase. The first definable point is the temperature at which the phase starts to precipitate, referred to as the stability temperature. The chemical formula of G-phase is Ni16Nb6Si7, where niobium will be needed to facilitate the precipitation of this intermetallic. The precipitation of G-phase is known to occur through a transformation mechanism with NbC [156]; However, if there is excess niobium at the dendrite boundaries the G-phase will be able to precipitate without consuming NbC precipitates. Excess niobium precipitating out as G-phase is indicative of the plateau region shown in  Figure 3a, where NbC remains unaffected. Once a certain temperature is reached M23C6 becomes stables, promoting the dissolution of NbC, freeing up the niobium to precipitate even more G-phase. Once all of the NbC has been dissolutioned and all the niobium has been exhausted in solution and precipitated as G-phase; the maximum, and terminal phase fraction for G-phase has been reached, denoted as point 3 on  Figure 3a. Point 3 must also be defined as the terminal phase fraction, meaning that for any subsequent decreases in temperature there will be only minor fluctuations in the equilibrium phase fraction. This will avoid any situations where the driving force of the examined phase would start to decrease with decreasing temperature, where it could eventually become metastable, and replaced with another phase.

For Matlab to be able to define these points for each chemistry in the design matrix they must be formulated as mathematical expressions. The stability temperature can be defined as mα > 0, where m is the phase fraction (mol%) for phase α. Identifying the Plateau region can be achieved by finding where ∂mα∕∂T 0, and 2mα∕∂T2 0. Finding the maximum terminal phase fraction is done in the same way as the plateau region with the added conditions: max(mα) , and at max(mα), mα∕∂T2 0T ∈{TxTrt}, where Tx is the temperature where max(mα).



Fig. 3: The volume fraction of each phase over a temperature range contains three notable regions; the stability temperature, the plateau region, and the maximum volume fraction. These three points are outlined for a) G-Phase, b) and M23C6 of an alloy with the composition Cr19-Ni31-Si0.5-Nb0.5-C0.05-Mn0.825-Ti0.05.

The Matlab distribution files for this section can be found at, under ‘precipitate_script’. The Matlab files are given separately, but you can download the full distribution by downloading the ThermoCalc_Script1.rar file. The two main scripts that run the program are allxls.m, and findeq.m. ‘allxls.mis used as a batch processing of all the .xlsx files produced by ThermoCalc, where each file in the directory listed on Line 1 of the script will go through the ‘findeq.mscript and then appended to a excel file for further use. ‘findeq.mis the script that extracts the information described above, shown in  Figure 3a, and  Figure 3b.

The first part of the script extracts the relevant data, and filters out all the unnecessary phases. Line 23 & 24 gets the row index for the row closest to T=400C, and T=100C. Line 29 extracts the data between this temperature range for the current phase in the for loop. The temperature range 100-400C was chosen as the intermetallic phases and the the M23C6 phase have been observed to be stable within this temperature range with ThermoCalc. This is mainly due to the underestimation in the driving force of G-phase predicted by ThermoCalc in the TTNI8 database. Once G-phase is added to an iron database (TTFE*), this range is likely to change to above 800C, as G-phase is seen to be stable in niobium stainless steel alloys at this temperature [52324]. The If statement in line 31 passes each phase to the next step if the average phase fraction between 100-400K is less than 7%, the maximum phase fraction is less than 10%, and the mean lower threshold is above zero. This is assuming an upper bound on the phases under question, which is a reasonable assumption for the 2032Nb precipitates as literature has never cited a phase fraction for any intermetallics (e.g. G-phase) or chromium carbides above 10 vol%. Setting these bounds also helps filter out strange FCC iterations (eg. FCC_A1#3, FCC_A1#4) that ThermoCalc predicts for very low temperatures, which have a very high driving force. Experimentally no other FCC phases other than MC carbides, and the solution phase were found, so it is appropriate to disregard these other FCC iterations. It should be noted that if this script is to be applied to other systems these conditions may need to be changed. The last conditions filters out the (Nb,Ti)(C,N) phase, or FCC_A1#2 phase, which is addressed in the second, separate script. finphase on line 34 collects the phase names of the phases that pass the criteria in Line 32. The next loop goes through each phase in finphase.


for phase = phases
 [~,y]=find(ismember(headertext, phase(1,1)));
 %get submatrix of phase from 400C to rt.

 if mean(subphase) < 0.07 && mean(subphase) > 6E-4 ...
	 && max(subphase)< 0.10 ...
	 && strcmp(phase(1,1), 'NP(FCC_A1#2)')==0 
		 finphase = [finphase phase(1,1)];

For each phase the derivative, and second derivative of the curve are calculated using the forwards difference method (i i + 1). The slope between points i, and i + a was also calculated where a is an arbitrary value used to determine if the curve has become constant for further decreases in temperature. For the experiments discussed in this thesis the value of a was set at 8.

% Gets index of current phase for data
[~,y]=find(ismember(headertext, phase(1,1)));
datPhase = data(:,y);
dn = length(datPhase)-1;
dphase = zeros(dn,1);
dphaseI = zeros(dn,1);
ddphase = zeros(dn-1,1);
dChange = 8;

for i=1 : dn-dChange
for i=1 : dn-dChange
for i=1 : dn-dChange
npm = data(:,y);
dtemp = temp(2:length(temp));

Line 63 creates an index n, for counting the number of points in dphase(i), which calculated the derivative of the original curve between i, and i + 8. mxval is the maxium phase fraction attained by the original curve, while mxIdx is the index of this value, and stabTemp is the index of first instance of mα > 0.

 n = length(dphase)-1;
 mxVal = max(npm);
 mxIdx = find(npm==mxVal, k,'first');
 stabTemp = find(npm,k,'first');

Next, the script needs to go through the phase fraction curve, and its derivative, and find the maximum terminal phase fraction, and if there is a plateau region. This can be achieved by going along each curve, first determining if the driving force for precipitation remains constant, or the phase fraction remains constant, between i and i + a, and then if this point is in a plateau region, or if it is the terminal phase fraction.

if (dphase(k,1) < 2E-5) && (abs(npm(k+1,1)) > 6E-4) ... 
	&& mean(npm(k:k+stbRng,1)) <= mxVal ...
	&& mean(npm(k:k+stbRng,1)) > mxVal-0.003 ...
	&& mean(ddphase(k:k+stbRng,1))< 2E-6...
	&& (dtemp(k) < 800)
elseif(dphase(k,1) < 2E-5) && (abs(npm(k+1,1)) > 6E-4) ... 
	&& mean(ddphase(k:k+stbRng,1))< 2E-6...
	&& (dtemp(k) < 800) && (dtemp(k) > 250)
	platArrPh = [platArrPh; {data(k-1,y)}];
	platArrT = [platArrT; {dtemp(k)}] ;
	isPlateau = 1; 

The If statement in line 78 determines if the driving force remains constant between k (which is just another iterative index similar to i), and k + stbRng, where stbRange is the a parameter, and if the phase fraction is terminal. This will pass if the first derivative, or the slope of the curve is close to zero (2E - 5), and if the phase fraction is within a threshold between mxV al (Maximum phase fraction), and an arbitrary lower threshold, mxV al - 0.003. The other conditions of this If statement determines if the second derivative, or the curvature between k, and k + stbRng (i i + a) is close to zero, and if the temperature is less than a certain temperature (eg. 800). The T < 800C was chosen from examination of ThermoCalc outputs, where G-phase, and M23C6 were never stable at temperatures greater than 800ºC. Once G-Phase is included in the TCFE* database this condition may not be true, as the actual transformation temperature for G-phase in creep resistant stainless steels is 852 ºC[25]. If the phase fraction passes these conditions it will be labeled as Stable, meaning a terminal phase is found, and it will break out of the For loop.

If these conditions are not met the phase fraction will go through another If statement on Line 103, which is almost the same as the If statement on Line 78 without the conditions that the phase fraction needs to be close to the maximum value mxV al. If these conditions pass they will be added to an array platArrPh which will contain the phase fractions for all points on the curve that meet these conditions, and another array platArrT which will contain the temperature. A boolean value called ‘isPlateauwill be given to this phase. Once the terminal phase fraction is found on Line 78, it will go through a gate labeling the phase as PlateauStableto indicate that a plateau was found, and the average of platArrPh, and platArrT will be taken and reported as the point where the plateau of the curve is found. If no terminal phase fraction is found, the maximum of the phase will simply be reported, and the phase will be labeled NotStableif no plateau is found, and PlateauNotif a plateau is found.

This script will be executed for each composition in the ’allxls.m’ file, and will output the results to a .xlsx file. An example of the resulting table is shown in Table 2. This table can now be reordered by grouping the data based on each phase for comparative purposes.

Table 2: Example of the output from the Matlab script ’findeq.m’, which looks at the intermetallic, and chrome carbide phases that might precipitate during long term aging. NPM represents the mole fraction of a phase.
NamePhaseT(Max Ph.Frac.)NPM(Max Ph.Frac.)Stability Temp Stability Ph.Frac.LabelT(Plateau)NPM(Plateau)
Cr19Ni31Si0.5Nb0.5C0.05Mn0.15Ti0.05'NP(G_PHASE)'206.850.017219565.850Plateau Stable450.850.002926
Cr19Ni31Si0.5Nb0.5C0.05Mn0.15Ti0.175'NP(G_PHASE)'212.94120.024289645.850Plateau Not450.850.010134
Cr19Ni31Si0.5Nb0.5C0.05Mn0.15Ti0.175'NP(MC)'287.58160287.58160Not Stable
Cr19Ni31Si0.5Nb0.5C0.05Mn0.15Ti0.3'NP(G_PHASE)'211.850.028125687.05550Plateau Not440.850.017198
Cr19Ni31Si0.5Nb0.5C0.05Mn0.15Ti0.3'NP(M23C6)'211.850.008444338.56160Not Stable
Cr19Ni31Si0.5Nb0.5C0.05Mn0.825Ti0.05'NP(G_PHASE)'232.38480.01602565.850Plateau Stable455.850.002844

2.2. Finding Maximum (Nb,Ti)(C,N) Phase fraction, composition, and terminal temperature

The second Matlab script will be used to compile useful information on (Nb,Ti)(C,N) for each chemistry in the compositional matrix. Since (Nb,Ti)(C,N) is of a cubic NaCl structure, which is the same structure as the solution phase (austenite), ThermoCalc creates the (Nb,Ti)(C,N) phase as a new composition set of the FCC phase. Therefore, in ThermoCalc the solution phase is denoted as FCC_A1#1, and (Nb,Ti)(C,N) is denoted as FCC_A1#2. For each chemistry iteration, FCC_A1#2 contains the end members (Nb,Ti)1(C,N)1, but the site fractions, or constituents fractions are uncertain. Furthermore, differentiation between NbC, and TiC is not made, and will be grouped into one phase called FCC_A1#2. This information cannot give us distinct phase fractions of NbC, and TiC, as the solubility of titanium in NbC, and niobium in TiC is unknown. For the purpose of this study it will be assumed that these two constituents (niobium, and titanium) are mutually exclusive to their respective phases, and a volume fraction of each will be presented. To determine the composition of the MC carbides, a plot consisting of the mole fraction of constituents of the FCC_A1#2 phase along the temperature range, Trt Tm, must be specified.  Figure 4b illustrates an example of a compositional plot of FCC_A1#2, which shows that the composition of FCC_A1#2 can fluctuate based on the solubility of the elements, and the stable phases in the system. To see how the solubilities of titanium and niobium in (Nb,Ti)(C,N) change with variations in alloy composition, the element fractions in FCC_A1#2 can be extracted at specific temperatures for comparison. As illustrated in  Figure 4b Matlab should extract the temperatures and the constituent fraction of FCC_A1#2 at the solubility limit of titanium, the maximum phase fraction of FCC_A1#2, and the dissolution temperature of FCC_A1#2. Due to an error in ThermoCalc reporting constituent fractions of a phase, these fractions do not do not drop to zero once the phase becomes unstable as seen in  Figure 4b.



Fig. 4: For the (Nb,Ti)(C) phase, or FCC_A1#2, a) the maximum phase fraction, and the dissolution temperature are extracted, along with b) specific composition of FCC_A1#2 at these regions, and the solubility limit of titanium.

The Matlab distribution files for extracting and compiling the NbC precipitates can be found here at, under ‘NbC_script. Odd anomalies from the composition files makes this script much more particular then the previous script, and should be reviewed before trying to modify it to work with other systems. Some tricks are used to avoid any errors in the data that might be unaccounted for, but some instances will still cause the script to error. If this script is providing trouble, and not completing, it is advised to cut certain sections that might not be needed. For example, there are four parts to this script that extract data from the given .xlsx files. The first part extracts the dissolution temperature and the maximum FCC_A1#2 fraction from the phase fraction data, the second part determines the solubility limit of titanium, the third part takes the maximum phase fraction and determines the respective amounts of NbC and TiC, and the fourth part finds the maximum amount of titanium in FCC_A1#2 and calculates the fractions of NbC and TiC at this point. The fourth part, finding the maximum titanium concentration, is not used in the presented study, but was included for analysis. This fourth part flags most of the errors, so it can be removed if it does not provide critical information to the study.

Tx = find(datPhase,k,'first'); % Finds first non-zero
Tend = find(datPhase==0, k,'last'); % Finds last zero for stability range
datPhaseSub = datPhase;
if(datPhase(Tend-1) > 0.01)
 datPhaseSub = datPhaseSub(1:Tend-1);
 Tend = find(datPhaseSub==0, k,'last');
if (Tend < Tx+30)
 Tend = length(datPhase);

Lines 14-24 determines the suitable range to extract the FCC_A1#2 data. This in an attempt to filter out any noise, or aberrations in the data. Tx finds the first non-zero row, assuming this is the stability temperature for FCC_A1#2, and Tend finds the temperature of the last zero row. Sometimes after dissolution of a phase, ThermoCalc predicts this phase to reprecipitate at low temperatures with a drastically different, often nonsensical compositions. It is important to filter out any regions in the FCC_A1#2 that do not have a composition of (Nb,Ti)(C). Setting Tend to the last zero row will filter out this data, unless FCC_A1#2 drops back down to zero after it has precipitated the second compositionally different phase. The condition on Line 18 will filter out any secondarily precipitated FCC_A1#2 phases by testing if the dissolutioned phase is above 0.01 mole fraction before dissolution, as these new phases usually have a very high driving force. If the FCC_A1#2 is stable until the end of the data, Tend is set to the last temperature in the data as per line 23. Variables dtemp, and npm are extracted between this temperature range (Tx Tend), which are the temperature array, and the phase fraction array respectively.

mxVal = max(npm);
mxIdx = find(npm==mxVal, k,'first');
mxTemp = dtemp(mxIdx-1);
if(Tend == length(datPhase))
 stabIdx = length(npm);
 stabIdx = find(npm==0,k,'first');  
stabTemp = dtemp(stabIdx-1);

Lines 31-39 will get the maximum FCC_A1#2 phase fraction, its index, and its corresponding temperature, as well as the dissolution temperature of FCC_A1#2 labelled as stabTemp. Niobium, and titanium mole fraction arrays are extracted and labeled as nbTot, and TiTot. The solubility limit of titanium is found on line 52, where the titanium concentration is below a lower threshold. A While loop introduced on Line 55 is used to find out if the found value is the actual solubility limit of titanium in FCC_A1#2 or if it was picking up some random fluctuation in the data. If the value picked up is a false positive it will be subtracted from the data, and a new solubility limit will be found until the criteria of the While loop is met. The catch statement between lines 59-74 will catch any errors that might occur in this while loop (e.g. indicies out of range), save the error to a log file, and break out of the script, and will continue running with the next file iteration. Lines 31-39 will get the maximum FCC_A1#2 phase fraction, its index, and its corresponding temperature, as well as the dissolution temperature of FCC_A1#2 labelled as stabTemp. Niobium, and titanium mole fraction arrays are extracted and labeled as nbTot, and TiTot. The solubility limit of titanium is found on line 52, where the titanium concentration is below a lower threshold. A While loop introduced on Line 55 is used to find out if the found value is the actual solubility limit of titanium in FCC_A1#2 or if it was picking up some random fluctuation in the data. If the value picked up is a false positive it will be subtracted from the data, and a new solubility limit will be found until the criteria of the While loop is met. The catch statement between lines 59-74 will catch any errors that might occur in this while loop (e.g. indicies out of range), save the error to a log file, and break out of the script, and will continue running with the next file iteration.

% Finds last non-zero for stability range
solIdx = find(tiTot >= 0.5*10^-3, k,'last');
tiTotSub = tiTot;
 while(tiTotSub(solIdx-1) < 0.5*10^-3)
  tiTotSub = tiTotSub(1:solIdx-1);
  solIdx = find(tiTotSub >= 0.5*10^-3, k,'last');
catch err
 %open file
 fid = fopen('logFile','a+');
 % write the error to file
 % first line: message

 % following lines: stack
 for e=1:length(err.stack)
  fprintf(fid,'%sin %s at %i\n',txt,err.stack(e).name,err.stack(e).line);

 % close file
solIdx = find(tiTot == tiTotSub(solIdx), k,'last');
tempSol = temp(solIdx);   % temp at the sol limit of Ti
phSol = datPhase(solIdx); % phase fraction at the sol limit of Ti.

Lines 79-84 splits the maximum phase fraction found before into NbC, and TiC by finding the concentrations of titanium and niobium in FCC_A1#2, and determining their respective phase fractions, assuming there is no solubility for titanium in NbC, and niobium in TiC.

%Max Phase fraction
tiMx = tiTot(mxIdx);
nbMx = nbTot(mxIdx);
totCon = tiMx + nbMx;
tiPh = tiMx/totCon*mxVal; %TiC Phase fraction
nbPh = nbMx/totCon*mxVal; % NbC Phase fraction

The last part of the script finds the maximum concentration of titanium, checking the validity of the number through a While statement, first checking if the found value is just a large fluctuation in the data, and that the value is below an upper threshold, knowing that the maximum concentration of titanium is never above 0.4. Lines 118-125 determine the respective NbC, and TiC concentrations at the found maximum titanium concentration. An example of the resulting output are shown in Table 3.

%Max Ti
tiTot = tiTot(1:length(tiTot)-1);
mxTi = max(tiTot);
mxTiIdx = find(tiTot==mxTi, k,'first');
tiTotSub = tiTot;
 while(abs(mxTi - tiTotSub(mxTiIdx+1)) > 1*10^-3 || mxTi > 0.5)
	if(mxTiIdx < length(tiTotSub)/2)
	 tiTotSub = tiTotSub(mxTiIdx+1:length(tiTotSub));
	 tiTotSub = tiTotSub(1:mxTiIdx-1);
	mxTi = max(tiTotSub);
	mxTiIdx = find(tiTotSub==mxTi, k,'first');
catch err

Table 3: Example of the output from the Matlab script ’findeqNbC.m’, which looks at the FCC_A1#2 phase or, (Nb,Ti)(C) phase, which extracts the maximum FCC_A1#2 phase fraction, and dissolution temperature, and determines the relative NbC and TiC phase fractions, as well as the solubility limit for titanium.
NameT(Mx.Ph. Frac.)NbC(Mx.Ph. Frac.)TiC(Mx.Ph.Frac.)T(Dissol)T T(Ti.Sol. Lmt.)NbC(Ti.Sol. Lmt.)T(Mx.Ti.)NbC(Mx.Ti.)TiC(Mx.Ti.)

2.3. Compiling and organizing data in Excel.

Once the Matlab scripts in section 2.1, and section 2.2 have successfully run and output their respective .xlsx files, download the excel file batchTCC.rarfrom this link, Once this file is open, go to the developer tab, and click the Visual Basic Button. In the sortMatlab script, lines 11-12 lines 16-17 refer to the excel files compiled by Matlab. This script will go through the specified files, reorganize, and categorize the data based on the phase and the type of data. The data compiled from section 8.1 is organized in the findPhase function, where the string array strPhase lists the phases that should be searched for. The data array stores the data found for each phase, where this data will be copied to its respective column in the excel table. The phase name, the maximum temperature, and maximum phase fraction, the stability temperature, and any plateau regions, are all stored in the data array. The next function findFCC does the same thing as the findPhase function, but applies it to the second excel file discussed in section 2.2. The organize, and organizeFCC functions take the data arrays returned from findPhase and findFCC, and sort them into their appropriate columns. To run the sortMatlab script press the the sort button on the excel table. Once that has completed the next step is to insert columns after the Name column for each element in the compositional matrix. In the Visual Basic Window click on Model 2 in the Modules folder located on the left-hand side project tree, and run the SplitName script. This takes each string in the Name column and splits the string as to fill the element columns just created. An example of the resulting excel file can be downloaded here, Finally select the header row click on the Data tab in the excel ribbon, and select the Filter option. Each column can now be arranged in ascending or descending order, or with various conditional statements. A supplemental graphing program is included in the .xlsm file, and can be executed by running the openGraphMaker function in Module 2. The X-,Y-,and Z- variables should be selected by clicking the header cells for the columns to be graphed. Up to three constant variables can be selected, by again clicking the header cell of the wanted column, and setting their values to the values you want to keep constant. A text file with the appropriate columns is output for further analysis with other plotting programs (e.g. GLE, or GNUPlot). Figure 4 shows how this plotting tool can be used with the resulting output.

Fig. 4: Excel output after running scripts in the batchTCC.rardistribution. Supplemental plotting tool is also displayed where data can be plotted quickly and easily for fast analysis.

3. Matrix Plots with GNUPlot

After the data is compile, and organized successfully (Columns represent the different phases, row represent the different chemistries), the data can now be used to output phase matrix plots, to visually represent the relationships between changes in composition and their equilibrium microstructure. Figure 5 is an example of a final matrix plot where influences from individual elements can be observed, made with GNUPlot. GNUPlot is an open source plotting program useful in converting large amounts of data, and displaying the information in various ways. GNUPlot for windows is distributed as an app package, meaning that for installation all that is needed is to extract the folder to somewhere in your computers directory. Then to open up GNUplot go to the Binary folder and double click the wgnuplot.exe executable. You can download all of the GNUplot matrix plot files discussed in this document by following the instructions in this link,

Fig. 5: A matrix plot looking at the maximum phase fraction of a 2032Nb microstructure with titanium additions, vs. compositional changes of individual elements. Relative effects from each component in the composition can be visually compared for analysis and discussion.

To run the GNUPlot (.gp) script, in the wgnuplot.exe environment type “cd path”, where path is the full directory where the .gp files are found (eg. C:\user\downloads). Then type “call filename”, where filename is the name of the .gp script you want to run (eg. The scripts should run and output a postscript file (.ps) of the resulting matrix plot. To view the code of the .gp file open it up in a text editor, or a coding environment such as Visual Studios, or Dreamweaver.

The first five lines of the GNUplot script set the output to a postscript file, and the height and width of the document. The name of the postscript file is “′′. To produce multiple plots in the same document the multiplot command is issued on line 7, creating a multiplot with 5 rows, and 8 columns. tmargin, bmargin, leftmargin, rmargin set the top, bottom, left, and right margins respectively for each plot. For GNUPlot to read the compiled excel data, the data should be saved as a .txt file (tab delimited) or a .csv file. Any columns that will not be used in the matrix plot should be deleted. It is also necessary to fill all empty empty cells with zeros, as empty cells cause indexing errors in GNUPlot. The zero cells can be filtered out in the script with conditional statements, as shown below.

set terminal postscript size 1650,750
set output ""
set term postscript enhanced font "Times-Roman" 12
unset key

set multiplot layout 5,8
 set tmargin 0.8
 set bmargin 0.8
 set lmargin 1.8
 set rmargin 1.8

The next code block will plot the first row of Figure 4. “MT_Phase.txt′′ is the excel data containing the maximum phase fractions, u is shorthand for ‘using’, where 1 : (100 * $$10) refers to the column1 : column2 indices of the .txt file. Post operations can be performed on the rows and columns where the data in row 10 is multiplied by 100 to get a percentage. Whenever doing a row or column operation the row and column must be referred to with $$. The command smoothuniquewithlines will plot a dotted line averaging each X-axis values which will be used to observe any trends in the data. Line 25 plots the Nb/C fraction as shown in Figure 4, where a line at Nb∕C = 7.7 is drawn by referencing a separate text file containing x, and y line coordinates. In this instance it is referring to line 1 {7.7,7.7}, and line 3 {0,2}.

 set format y "%g"
 set ytics 1 nomirror
 set ylabel "NbC(at%)"
 plot "MT_Phase.txt" u 1:(100*$$10), "" u 1:(100*$$10) smooth unique with lines
 set format y ""
 unset ylabel
 plot "MT_Phase.txt" u 2:(100*$$10), "" u 2:(100*$$10) smooth unique with lines
 plot "MT_Phase.txt" u 3:(100*$$10), "" u 3:(100*$$10) smooth unique with lines
 plot "MT_Phase.txt" u 4:(100*$$10), "" u 4:(100*$$10) smooth unique with lines
 plot "MT_Phase.txt" u 5:(100*$$10), "" u 5:(100*$$10) smooth unique with lines
 plot "MT_Phase.txt" u 6:(100*$$10), "" u 6:(100*$$10) smooth unique with lines
 plot "MT_Phase.txt" u 7:(100*$$10), "" u 7:(100*$$10) smooth unique with lines
 plot "MT_Phase.txt" u 8:(100*$$10), "Nb_C.txt" u 1:3 with l lt 1 lw 1

Line 13-25 are repeated until the final row where the x-axis needs to be generated. If there is a clear division in the data like that shown for the M6C row in Figure 5, the division can be analyzed by filtering the data in excel (explained above), and coloring the rows by selecting the row, and clicking the Conditional Formatting button in the Home tab, and choosing a Color Scales option. In this case there is a clear division when Nb∕C 25. To filter the data and draw lines for when Nb∕C > 25, and Nb∕C < 25 a conditional statement given on Line 57 is given, which says that if column 8 (Nb/C row) is greater than 25 and column 14 (M6C row) is greater than zero then print the value in column 14, else print 10, where 10 is an undefined value and will not be included in the interpolation line. In coding vernacular, ?, is another way of writing an If statement, and the proceeding : is an Else statement.

plot "MT_Phase.txt" u 1:($$14), "" u 1:((($$8 > 25) && ($$14 > 0))? $$14 : 1/0) smooth unique with lines, "" u 1:((($$8 < 25) && ($$14 > 0))? $$14 : 1/0) smooth unique with lines

For the last row of the plot the x-axis tics need to be set, along with the x-axis label, and the x-axis range.To end the script unset the x-axis, y-axis, and multiplot, and reset the default conditions. If the script should error midway in order to run the script again you may need to manually input the commands on lines 100-102 in order to get out of the multiplot environment. For students/members of the CCWJ more information about GNUPlot can be found here,, double clicking textbooks and then Gnuplot.

set format y "%g"
set ylabel "G-Phase (at%)"
set ytics 1 nomirror
set xtics 1 nomirror
set xlabel "Cr (wt%)"
plot "MT_Phase.txt" u 1:(100*($$15)), "" u 1:(100*($$15)) smooth unique with lines
set format y ""
unset ylabel
set xlabel "Ni (wt%)"
set xtics 31,1.5,34
plot "MT_Phase.txt" u 2:(100*($$15)), "" u 2:(100*($$15)) smooth unique with lines
	unset xlabel
	unset ylabel
unset multiplot
set key

4. References

[1]    Shi, S., Lippold, J.. Microstructure evolution during service exposure of two cast, heat-resisting stainless steels – hp-nb modified and 20-32nb. Mater Charact 2008;59(8):1029–1040.

[2]    Berghof-Hasselbcher, E., Gawenda, P., Schorr, M., Schtze, M., Hoffman, J.. Atlas of Microstructures. Materials Technology Institute; 2008.

[3]    Nishimoto, K., Saida, K., Inui, M., Takahashi, M.. Changes in microstructure of hp-modified heat-resisting cast alloys with long term aging. repair weld cracking of long term exposed hp-modified heat-resisting cast alloys. (report 2). Quarterly Journal of the Japan Welding Society 2000;18(3):449–458.

[4]    Nishimoto, K., Saida, K., Inui, M., Takahashi, M.. Mechanism of hot cracking in haz of repair weldments. repair weld cracking of long term exposed hp-modified heat-resisting cast alloys. (report 3). Quarterly Journal of the Japan Welding Society 2000;18(4):590–599.

[5]    Powell, D.J., Pilkington, R., Miller, D.A.. The precipitation characteristics of 20% cr/25% ni—nb stabilised stainless steel. Acta Metall 1988;36(3):713–724.

[6]    Chen, Q.Z., Thomas, C.W., Knowles, D.M.. Characterisation of 20cr32ni1nb alloys in as-cast and ex-service conditions by sem, tem and edx. Mater Sci Eng, A 2004;374(1-2):398–408.

[7]    Danielsen, H.K., Hald, J.. On the nucleation and dissolution process of z-phase cr(v,nb)n in martensitic 12%cr steels. Mater Sci Eng, A 2009;505(1-2):169–177.

[8]    Danielsen, H., Hald, J.. Influence of z-phase on long-term creep stability of martensitic 9-12%cr steels. In: 9th Liege Conference on Materials for Advanced Power Engineering. 2010, p. 310.

[9]    Sourmail, T., Bhadeshia, H.. Microstructural evolution in two variants of nf709 at 1023 and 1073 k. Metall Mater Trans A 2005;36(1):23–34.

[10]    Erneman, J., Schwind, M., Liu, P., Nilsson, J.O., Andrén, H.O., Ågren, J.. Precipitation reactions caused by nitrogen uptake during service at high temperatures of a niobium stabilised austenitic stainless steel. Acta Mater 2004;52(14):4337–4350.

[11]    Pickering, F.B., Keown, S.. Niobium in stainless steels. In: Stuart, H., editor. Niobium, Proceedings of the International Symposium. Warrendale, PA: Met. Soc. AIME; 1981, p. 1113–1142.

[12]    Sourmail, T.. Literature review precipitation in creep resistant austenitic stainless steels. Mater Sci Technol 2001;17(January):1–14. URL

[13]    Xiao, B., Xing, J.D., Feng, J., Zhou, C.T., Li, Y.F., Su, W., et al. A comparative study of cr 7 c 3 , fe 3 c and fe 2 b in cast iron both from ab initio calculations and experiments. Journal of Physics D: Applied Physics 2009;42(11):115415.

[14]    Jo, T.S., Lim, J.H., Kim, Y.D.. Dissociation of cr-rich m23c6 carbide in alloy 617 by severe plastic deformation. J Nucl Mater 2010;406(3):360–364.

[15]    Holman, K.L., Morosan, E., Casey, P.A., Li, L., Ong, N.P., Klimczuk, T., et al. Crystal structure and physical properties of mg6cu16si7-type m6ni16si7, for m = mg, sc, ti, nb, and ta. Mater Res Bull 2008;43(1):9–15.

[16]    Hans Lukas, S.G.F., Sundman, B.. Computational Thermodynamics. Cambridge University Press; 2007.

[17]    Hillert, M.. Phase Equilibria, Phase Diagrams and Phase Transformations: Their Thermodynamic Basis. Cambridge University Press; 2007.

[18]    Liu, Z.K.. Computational thermodynamics using thermo-calc; 2010. ThemoCalc & Dictra Training Course.

[19]    Song, Y.Y.. Thermodynamic study on b and fe substituted cr23c6 using first-principles calculations. Ph.D. thesis; Pohang University of Science and Technology; 2010.

[20]    Hoffman, J., Magnan, J.. Cast 20cr32ni1nb alloy aged mechanical property improvements via chemistry modifications. In: Corrosion 2003. no. 3469; NACE International; 2003,.

[21]    Vitek, J.. G-phase formation in aged type 308 stainless steel. Metallurgical and Materials Transactions A 1987;18(1):154–156.

[22]    Danielsen, H.K., Hald, J.. A thermodynamic model of the z-phase cr(v, nb)n. Calphad 2007;31(4):505–514.

[23]    Shi, S., Lippold, J., Ramirez, J.. Hot ductility behavior and repair weldability of service-aged, heat-resistant stainless steel castings. Weld J 2010;89(10):210–217.

[24]    Shibasaki, T., Mohri, T., Takemura, K.. Experience with cast material for steam reformer furnaces. Ammonia Plant Saf Relat Facil 1994;34:166–176.

[25]    Hoffman, J.. High temperature aging characteristics of 20cr32ni1nb castings. In: Corrosion 2000. no. 512; NACE International; 2000,.

[26]    Montgomery, D., Runger, G.. Applied Statistics and Probability for Engineers, 4th Edition. John Wiley & Sons; 2006.

[27]    Spliid, H.. Design and analysis af experiments with k factors having p levels; 2002.

[28]    Minitab, . Technical support document rank deficiency; 2010. /RankDeficiency.pdf.


5.Glossary *

Element, or species, that occupies a specific sublattice of a specific phase. A phase can also be considered as a constituent of the total system.. 5, 14


The final chemical formula of a stable or metastable phase whose sublattice(s) are occupied by single constituents. For example M23C6 is an end member of (Cr,Ni,Fe,Nb)23C6.. 8


The independent variable of a factorial design. 23


main effect
How much the change in an individual factor effects the change in the response variable of a factorial design.. 24


Independent repetition of a treatment in a factorial experiment. 26
response variable
The dependent variable of a factorial experiment, or a regression model.. 24


A computational thermodynamics program that can calculate equilibrium phase diagrams for multicomponent systems, as well as Scheil simulations, and various thermodynamic properties (Cp, ΔHm, ΔGm etc...). 3
A specific level of a factor in a factorial design.. 24

6. Acronyms *

Analysis of variance. 24, 25


Compound-energy formalism. 8


Long Range Ordering. 8


7. Nomenclature *

the effect of the ith level of factor ‘B’
random error component
for all instances of ...
in a set ...
Overall mean effect
Chemical potential of component or end-member i
the effect of the ith level of factor ‘C’
the effect of the ith level of factor ‘A’
total Gibbs energy; G = αmα ·Gmα
partial Gibbs energy of component i in phase α; Giα = (  α)
integral molar Gibbs energy of a phase
constituent array of order i
interaction parameter of compound I
fraction of a phase
moles of component i
gas constant, 8.314Jmol-1K-1
coefficient of multiple determination
molar entropy of a phase
Temperature (K)
total mol fraction of component i; xi = αmα · xiα
mole fraction of component i in phase α

© 2012 CCWJ. All Rights Reserved.