How to iterate until convergence?
30 views (last 30 days)
Show older comments
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...
0 Comments
Accepted Answer
Torsten
on 12 Feb 2019
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.
0 Comments
More Answers (1)
Torsten
on 18 Dec 2018
Edited: Torsten
on 18 Dec 2018
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
See Also
Categories
Find more on Bodies in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!