How to iterate until convergence?

I need some help with the following iteration problem. I am trying to calculate the dendrite tip radius (R), which depends on several parameters (PT and PC for example) that also contain the variable R. I want to iterate until the value for R converges. Right now, I have to manually update the initial value for R with the new value that is calculated using the equation.

To add to the confusion, some of the variables needed for the dendrite tip radius calculation (k, mv, ksi_L, etc.) depend on the sample solute concentration (C0) which are in column vector form. For example, C0(1,1) is the solute concentration for Chromium, so k(1,1) should be the non-equil. partition correction for Chromium. This is the same for each of the solutes. In the equation for R, the column vectors should be multiplied together and summed.

My code so far:

% Undercooling calculation 
clear
clc

% Constants: Material Properties

alpha = 5.46E-6;    % thermal diffusivity
D0 = 2.0E-9;    % solute diffusivity
a0 = 2.358E-10;    % atomic spacing in liquid
sigma = 0.319;  % solid liquid interface energy
delta_Sf = 1020485;    % entropy of fusion
omega = 7.8384E-6;    % molar volume 
CLp = 5749190;   % heat capacity
C0 = [0.2118; 0.1868; 0.0312; 0.0181; 0.0115; 0.0107; 0.0015]; % sample solute concentration 
n = length(C0);
mL = [-2.74; -0.29; -20.48; -9.42; -16.06; -4.68; -11.0]; % equil. liquidus slope
ke = [0.961; 1.084; 0.183; 0.638; 0.406; 0.843; 0.123]; % equil. partition coefficient
delta_Hf = [20500; 13800; 26800; 36000; 18700; 10700; 105000]; % heat of fusion

% Inset Initial Guess:

V = 1.05; % growth front velocity
R = 1.54E-7; % dendrite tip radius (initial guess)

% Dendrite tip radius iteration calculation:

sigma_star = 1/(4*pi^2); % tip radius stability constant
Vd = D0/a0; % atomic diffusive speed
gamma = sigma/delta_Sf; % Gibbs-Thomas coefficient
PT = (V*R)/(2*alpha); % thermal peclet number
PC = (V*R)/(2*D0); % solutal peclet number
ksi_L = 1 - 1/(sqrt(1+(1/(sigma_star*PT^2)))); % tip radius stability parameter

fun = @(u) exp(-u)./u;
upper = Inf;
Iv_PC = PC*exp(PC)*integral(fun,PC,upper); % Ivantsov function of the termal Peclet number

k = ones(7,1); % non-equil. partition correction 
mv = ones(7,1); % non-equil. liquidus slope correction
ksi_C = ones(7,1); % tip radius stability parameter

for i = 1:n
    k(i,1) = (ke(i,1)+(V/Vd)) ./ (1-((1-ke(i,1)).*C0(i,1))+(V/Vd));
    mv(i,1) = mL(i,1) * (1+(ke(i,1)-k(i,1)*(1-log(k(i,1)/ke(i,1)))))/(1-ke(i,1));
    ksi_C(i,1) = 1 + (2*k(i,1))/(1-(2*k(i,1))-sqrt(1+(1/(sigma_star*PC^2)))); % tip radius stability parameter
    
end 

 R = (gamma/sigma_star) / ...
     ( sum( (PT*ksi_L.*delta_Hf) / (omega*CLp) ) - ...
     sum( (2*PC.*mv.*C0.*(1-k).*ksi_C) ./ (1-(1-k)*Iv_PC) ) ) % dentrite tip radius 
 

If the above code is hard to follow, or is not very clear, I have attached a pdf outlining what I am trying to do with the equations written out. Note, if you look at the pdf, I just need help with the first part, the dendrite tip radius calculation. I did the initial calculation, the calculation shown in the pdf, using a program called Mathcad, and I am certain that the equations are correct

Thanks for the help!

P.S. I am a little paranoid when it comes to calculations, so I like to use lot of parties...

Answer 1

[Torsten]

The structure of your code should look like:

function main

  V = linspace(0.01,0.25,1000);
  R = zeros(numel(V),1);
  R0 = 1.0e-2;

  for m = 1:numel(V)
    Veloc = V(m);
    R(m) = fzero(@(R)fun_iter(R,Veloc),R0);
    R0 = R(m)
  end

end


function res = fun_iter(R,Veloc)

% Constants: Material Properties 

    alpha = 9.7E-5;    % thermal diffusivity
    D0 = 2.0E-9;    % solute diffusivity
    sigma = 0.167;  % solid liquid interface energy
    delta_Sf = 1.13E06;    % entropy of fusion

    G = 2.5E07;   %thermal gradient
    C0 = [0.37369; 9.540786]; % sample solute concentration 
    n = length(C0);
    mL = [-5.1599; -6.5792]; % equil. liquidus slope
    ke = [0.30974; 0.1227]; % equil. partition coefficient
    delta_C = [0.83; 10.97];   %-1.37; -2.92; 2.59

    % Dendrite tip radius iteration calculation:

    gamma = sigma/delta_Sf; % Gibbs-Thomas coefficient
    PT = (Veloc*R)/(2*alpha); % thermal peclet number
    PC = (Veloc*R)/(2*D0); % solutal peclet number
    Iv_PC = PC*exp(PC)*expint(PC); % Ivantsov function of the termal Peclet number
    Iv_PT = PT*exp(PT)*expint(PT); % Ivantsov function of the termal Peclet number

    Gi = ones(2,1); % non-equil. liquidus slope correction
    ksi_i = ones(2,1); % tip radius stability parameter

    constant = sqrt((1+((4*pi*D0)/(R*Veloc))^2)^(1/2));   % constant for the tip radius stability parameter equation, used to simplify

    for i = 1:n
        Gi(i,1) = (-Veloc.*delta_C(i,1)) / ( D0.*(1-(1-ke(i,1)).*Iv_PC) );  % composition gradient of the solute i ahead of the interface
        ksi_i(i,1) = 1 - ((constant) ./ (1 - constant - 2*ke(i,1)) ); % tip radius stability parameter  
    end

    R1 = 2*pi*(sqrt( gamma / ( sum(mL.*Gi.*ksi_i)-G ) )); % dentrite tip radius 
    res = R - R1;

end

But check your equations ; I get negative values as argument to the sqrt in the calculation of R1.

Best wishes

Torsten.

Answer 2

[Torsten]

function main
R0 = 1.0e-6;
R = fzero(@fun_iter,R0)
end

function res = fun_iter(R)
% Undercooling calculation 

% Constants: Material Properties

alpha = 5.46E-6;    % thermal diffusivity
D0 = 2.0E-9;    % solute diffusivity
a0 = 2.358E-10;    % atomic spacing in liquid
sigma = 0.319;  % solid liquid interface energy
delta_Sf = 1020485;    % entropy of fusion
omega = 7.8384E-6;    % molar volume 
CLp = 5749190;   % heat capacity
C0 = [0.2118; 0.1868; 0.0312; 0.0181; 0.0115; 0.0107; 0.0015]; % sample solute concentration 
n = length(C0);
mL = [-2.74; -0.29; -20.48; -9.42; -16.06; -4.68; -11.0]; % equil. liquidus slope
ke = [0.961; 1.084; 0.183; 0.638; 0.406; 0.843; 0.123]; % equil. partition coefficient
delta_Hf = [20500; 13800; 26800; 36000; 18700; 10700; 105000]; % heat of fusion

% Inset Initial Guess:

V = 1.05; % growth front velocity
%R = 1.54E-7; % dendrite tip radius (initial guess)

% Dendrite tip radius iteration calculation:

sigma_star = 1/(4*pi^2); % tip radius stability constant
Vd = D0/a0; % atomic diffusive speed
gamma = sigma/delta_Sf; % Gibbs-Thomas coefficient
PT = (V*R)/(2*alpha); % thermal peclet number
PC = (V*R)/(2*D0); % solutal peclet number
ksi_L = 1 - 1/(sqrt(1+(1/(sigma_star*PT^2)))); % tip radius stability parameter

Iv_PC = PC*exp(PC)*expint(PC); % Ivantsov function of the termal Peclet number

k = ones(7,1); % non-equil. partition correction 
mv = ones(7,1); % non-equil. liquidus slope correction
ksi_C = ones(7,1); % tip radius stability parameter

for i = 1:n
    k(i,1) = (ke(i,1)+(V/Vd)) ./ (1-((1-ke(i,1)).*C0(i,1))+(V/Vd));
    mv(i,1) = mL(i,1) * (1+(ke(i,1)-k(i,1)*(1-log(k(i,1)/ke(i,1)))))/(1-ke(i,1));
    ksi_C(i,1) = 1 + (2*k(i,1))/(1-(2*k(i,1))-sqrt(1+(1/(sigma_star*PC^2)))); % tip radius stability parameter  
end 

 R1 = (gamma/sigma_star) /...
     ( sum( (PT*ksi_L.*delta_Hf) / (omega*CLp) ) -...
     sum( (2*PC.*mv.*C0.*(1-k).*ksi_C) ./ (1-(1-k)*Iv_PC) ) ); % dentrite tip radius 

res = R - R1;

end

[Meelap Coday]

Thanks for the help so far, but I need some more help with the code. The above code you have written works great when considering only a single value for V. I have been trying to modify the code so that I can input a range of values for V, the solidification rate, and output a range of values for R corresponding to the values of V. I want to create a graph that shows the Solidification Rate vs. Dendrite Tip Radius (Solidification rate on the Y-axis and Dendrite Tip Radius on the X-Axis). Here is what I have so far:

function main_V2(R0)

R = ones(1,1000);
lenR = length(R);

R0 = 1.0e-8;

for m = 1:lenR
    R(1,m) = fzero(@fun_iter,R0)
end

end

function res = fun_iter(R)
% Constants: Material Properties

% ------------------------ growth front velocity --------------------------
% Enter in a value for the velocity:
% -------------------------------------------------------------------------
V = linspace(0.01,0.25,1000);
lenV = length(V);

for k = 1:lenV
    Veloc = V(1,k); 

    alpha = 9.7E-5;    % thermal diffusivity
    D0 = 2.0E-9;    % solute diffusivity
    sigma = 0.167;  % solid liquid interface energy
    delta_Sf = 1.13E06;    % entropy of fusion
    
    G = 2.5E07;   %thermal gradient
    
    C0 = [0.37369; 9.540786]; % sample solute concentration 
    n = length(C0);
    mL = [-5.1599; -6.5792]; % equil. liquidus slope
    ke = [0.30974; 0.1227]; % equil. partition coefficient
    delta_C = [0.83; 10.97];   %-1.37; -2.92; 2.59
    
    % Dendrite tip radius iteration calculation:
    
    gamma = sigma/delta_Sf; % Gibbs-Thomas coefficient
    PT = (Veloc*R)/(2*alpha); % thermal peclet number
    PC = (Veloc*R)/(2*D0); % solutal peclet number
    Iv_PC = PC*exp(PC)*expint(PC); % Ivantsov function of the termal Peclet number
    Iv_PT = PT*exp(PT)*expint(PT); % Ivantsov function of the termal Peclet number
    
    Gi = ones(2,1); % non-equil. liquidus slope correction
    ksi_i = ones(2,1); % tip radius stability parameter
    
    constant = sqrt((1+((4*pi*D0)/(R*Veloc))^2)^(1/2));   % constant for the tip radius stability parameter equation, used to simplify
    
    for i = 1:n
        Gi(i,1) = (-Veloc.*delta_C(i,1)) / ( D0.*(1-(1-ke(i,1)).*Iv_PC) );  % composition gradient of the solute i ahead of the interface
        ksi_i(i,1) = 1 - ((constant) ./ (1 - constant - 2*ke(i,1)) ); % tip radius stability parameter  
    end
    
    R1 = 2*pi*(sqrt( gamma / ( sum(mL.*Gi.*ksi_i)-G ) )); % dentrite tip radius 
    
     res = R - R1;
 
end
end

Although the core of the code has not changed since the version shown above, I have changed some of it. The only difference is that I changed the equation I used to calculate R. Additionally, I have also changed some of the constants used as a results of a change in the R equation.

Ultimately, I need to use the values of R and V to find a temperature value, and create another graph that shows the relationship between the solidification rate and the temperature (Solidification rate on the Y-axis and Temperature on the X-Axis). When I use the function, it does not save any of the value in my workspace, only displaying a value of R in the command window. Furthermore, when I try to call the function in a different file of code, it give me the following error:

R = main_V2(1.08E-8); %This is what I input in a different file of code.

Error using main_V2
Too many output arguments.

Error in AlSi10Mg_Undercooling (line 10)
R = main_V2(1.08E-8);

Here is the code I have that I used to calculate and display the values for the temperature. Since the constants and equations used in the radius equation are not saved to the workspace, I have to re-input those values and equations. Is there a way that I can combine the radius calculation and temperature calculation into one code, or at least call the radius function in the temperature code?

R = main_V2(1.08E-8);
V = 0.17149;

alpha = 9.7E-5;    % thermal diffusivity
D0 = 2.0E-9;    % solute diffusivity
sigma = 0.167;  % solid liquid interface energy
delta_Sf = 1146146;    % entropy of fusion
omega = 9.99E-6;    % molar volume 
CLp = 2443050;   % heat capacity
V0 = 2347.53;   % kinetic rate parameter
Rbar = 8.314;    % universal gas constant
Tm = 853;    % melting temperature
G = 2.5E07;   % thermal gradient

C0 = [0.37369; 9.540786]; % sample solute concentration 
n = length(C0);
mL = [-5.1599; -6.5792]; % equil. liquidus slope
ke = [0.30974; 0.1227]; % equil. partition coefficient
delta_Hf = [13200; 50200]; % heat of fusion
delta_C = [0.83; 10.97];   %-1.37; -2.92; 2.59


% -------- Miscellaneous calculations for to find undercooling: -----------

gamma = sigma/delta_Sf; % Gibbs-Thomas coefficient
PT = (V*R)/(2*alpha); % thermal peclet number
PC = (V*R)/(2*D0); % solutal peclet number
Iv_PC = PC*exp(PC)*expint(PC); % Ivantsov function of the termal Peclet number
Iv_PT = PT*exp(PT)*expint(PT); % Ivantsov function of the termal Peclet number


% ---------------------- Undercooling calcualtion -------------------------

% thermal undercooling:
delta_TT = (delta_Hf*Iv_PT)/(omega*CLp); 
delta_TT = sum(delta_TT);

% curvature undercooling:
delta_TR = 2*(gamma/R); 

% solutal undercooling:
Cl_i = ones(2,1);
for i = 1:n
    Cl_i(i,1) = C0(i,1) + delta_C(i,1)*(Iv_PC / (1-Iv_PC*(1-ke(i,1)))); 
    delta_TC = -mL(i,1).*(Cl_i(i,1)-C0);
end
delta_TC = sum(delta_TC);

% kinetic undercooling:
mu = (delta_Hf*V0)/(Rbar*Tm^2); 
delta_TK = V./mu;
delta_TK = sum(delta_TK);

% total undercooling:
delta_T = delta_TT + delta_TR + delta_TC + delta_TK;


% -------------------------- Display results ------------------------------

% format short

text_R = sprintf('The dendrite  tip radius is: %0.4d m.\n',R);
disp(text_R)

text_T = sprintf('The thermal undercooling is: %1.4d K.\nThe curvature undercooling is: %.4d K. \nThe solutal undercooling is: %.4d K. \nThe kinetic undercooling is: %.4d K. \n\nThe total undercooling is : %.4d K' , delta_TT, delta_TR, delta_TK, delta_TC, delta_T);
disp(text_T)