|
DataMuseum.dkPresents historical artifacts from the history of: RC4000/8000/9000 |
This is an automatic "excavation" of a thematic subset of
See our Wiki for more about RC4000/8000/9000 Excavated with: AutoArchaeologist - Free & Open Source Software. |
top - metrics - download
Length: 18432 (0x4800)
Types: TextFile
Names: »mprofile15«
└─⟦00964e8f7⟧ Bits:30007478 RC8000 Dump tape fra HCØ.
└─⟦b2ec5d50f⟧
└─⟦this⟧ »mprofile15«
begin
comment MPROFILE, VER 15, 1-OCT-81, Kaare Lund Rasmussen, RC4000.;
comment VER 15 : no plotting, undercooling assumed in simulation
one-sided growth only. Local Ni and P for input of
calculated loacl P as phosfor;
comment VER 15 : phase diagram Romig and Goldstein (1980).
Corrected solution of differential equation.
Phosphorous included in phasediagram, both alfa and gamma
and in diffusion coefficients, according to Moren
and Goldstein (1979).;
comment Mainprogram for cooling simulation above 320 deg C;
integer i,j,k,maxno,plotseq,ordn,x1,x2,x2a,start,ch,nomeas,noalfa,
alfainit,oldq,agrellgrow,oldx2,oldx2a,chp,chpsat,aga;
real t,tempstart,tempslut,bulknipc,q,qa,dx1,dx2,xbound1,xbound2,pl,
dt,r1,r2,rod,dxg1,dxg2,kk,aa,bb,cc,skt,xzero1,xzero2,xbound,stable,
erral,err1,err2,gamw,alfapre,alfanex,theta,kappa1,kappa2,phosfor,
xmpc,mpc,tzero,undercool,understart,nialfapre,dxal,dxgal,
nipcalint,phosalfa,phosgam,nipcalfa,nipcgam,phosalsat,
nialsat,phosgasat,nigasat,nial0,niga0,phosal0,phosga0,
phoschem,localp,tee1,tee2,tee,secinyear,wfe,wni;
real array nipc,nipci,nipc2,xx1,xg,s,nipcal,xa(1:200),
text(1:4),xmeas,nimeas(1:60),p,talfa,phgam,phalfa,
nialph0,nialphinf,nigaph0,nigaphinf(0:20),
const(1:10);
real procedure dgam(nipc);
value nipc;
real nipc;
begin
comment this procedure computes the diffusion
koefficient of ni in gamma iron at
temperature t;
real ni;
ni := nipc*100.0/(wni*(nipc/wni+(100.0-nipc)/wfe));
dgam := exp(0.0519*ni + 1.15)*(1.0 + 9.3*phosgam)*
exp(-(76400 - 11.6*ni)/(1.986*t))*secinyear;
end of real procedure dgam;
real procedure dgamP(temp,nipc,pres);
value temp,nipc;
real temp,nipc,pres;
begin
comment Pressure correction to diffusion constant
in the gamma phase. Deltavol = activation volume.;
real d,deltavol;
deltavol := 0.0;
d := dgam(nipc);
d := d*exp(-pres*deltavol/(1.986*temp));
dgamP := d;
end of real procedure dgamP;
real procedure dalfa;
begin
comment this procedure computes the diffusion
koefficient of nickel in alpha iron
at temperature t, according to Romig
and Goldstein (1979).;
dalfa := 10.5*exp(-64300/(1.986*t))*secinyear*
(1.0 + 1.2665*phosalfa + 0.6231*phosalfa*phosalfa);
end of real procedure dalfa;
real procedure dalfaP(temp,pres);
value temp;
real temp,pres;
begin
comment Pressure correction to diffusion constant
in the alpha phase. Deltavol = activation volume.;
real d,deltavol;
deltavol := 0.0;
d := dalfa;
d := d*exp(-pres*deltavol/(1.986*temp));
dalfaP := d;
end of real procedure dalfaP;
procedure phasestart(phosfor);
value phosfor;
real phosfor;
begin
comment This procedure initialize the phasediagram.;
noalfa := 7;
talfa(0) := 1187; talfa(1) := 973; talfa(2) := 873;
talfa(3) := 773; talfa(4) := 748; talfa(5) := 723;
talfa(6) := 673; talfa(7) := 573;
comment alpha-phase;
nialph0(0) := 0.0; nialph0(1) := 4.1; nialph0(2) := 5.4;
nialph0(3) := 6.05; nialph0(4) := 6.1; nialph0(5) := 6.0;
nialph0(6) := 5.5; nialph0(7) := 4.1;
nialphinf(0) := 0.00; nialphinf(1) := 4.5; nialphinf(2) := 6.1;
nialphinf(3) := 7.5; nialphinf(4) := 7.64; nialphinf(5) := 7.8;
nialphinf(6) := 7.2; nialphinf(7) := 4.7;
comment gamma-phase;
nigaph0(0) := 0.0; nigaph0(1) := 10.4; nigaph0(2) := 20.2;
nigaph0(3) := 33.3; nigaph0(4) := 35.4; nigaph0(5) := 39.6;
nigaph0(6) := 46.0; nigaph0(7) := 54.6;
nigaphinf(0) := 0.0; nigaphinf(1) := 9.3; nigaphinf(2) := 15.7;
nigaphinf(3) := 27.3; nigaphinf(4) := 29.6; nigaphinf(5) := 32.5;
nigaphinf(6) := 37.5; nigaphinf(7) := 46.0;
end of procedure phasestart;
procedure phasedia(temp,ordn,start,undercool);
value temp,undercool;
real temp,undercool;
integer ordn,start;
begin
comment This procedure evalueates the phasediagram
at temperature temp. The equilibrium tieline of
Moren and Goldstein (1979) are used.;
real ni,ni0,niinf,pho,phinf,ta1,ta2,a,b,c,ecxp;
integer i,n;
comment Saturation and zero Ni and P values are found.;
n := 0;
for i := 1 step 1 until noalfa do
if temp < talfa(i) then n := i;
if n = noalfa then n := n-1;
ta1 := talfa(n); ta2 := talfa(n+1);
nial0 := (nialph0(n+1) - nialph0(n))/(ta1 - ta2);
nial0 := nialph0(n) + nial0*(ta1 - temp);
if nial0 < 0.0 then nial0 := 0.0;
nialsat := (nialphinf(n+1) - nialphinf(n))/(ta1-ta2);
nialsat := nialphinf(n) + nialsat*(ta1-temp);
if nialsat < 0.0 then nialsat := 0.0;
niga0 := (nigaph0(n+1) - nigaph0(n))/(ta1-ta2);
niga0 := nigaph0(n) + niga0*(ta1-temp);
if niga0 > 50.0 then niga0 := 50.0;
nigasat := (nigaphinf(n+1) - nigaphinf(n))/(ta1-ta2);
nigasat := nigaphinf(n) + nigasat*(ta1-temp);
phosalsat := 10**(-2440.0/temp + 2.346);
phosgasat := 10**(-2790.0/temp + 2.29);
comment The real values are found now;
if chpsat = 0 then
begin
c := localp*(niga0 - nial0);
b := (phosgasat*(nial0-bulknipc)
+localp*(nigasat-niga0))
+(phosalsat*(bulknipc-niga0)
-localp*(nialsat-nial0));
a := (phosgasat*(nialsat-nial0)
-phosalsat*(nigasat-niga0));
tee2 := b*b - 4.0*a*c;
if tee2 < 0.0 then tee1 := -b
else tee1 := - b - sqrt(tee2);
tee := tee1*0.5/a;
if tee > 1.0 then tee := 1.0;
nipcgam := niga0 + tee*(nigasat-niga0);
nipcalfa := nial0 + tee*(nialsat-nial0);
phosgam := tee*phosgasat;
if phosgam < 0.0 then phosgam := 0.0;
phosalfa := tee*phosalsat;
if nipcgam < bulknipc then nipcgam := bulknipc;
if tee = 1.0 and x2 > 5 then chpsat := 1.0;
comment Maximum P obeyed;
ecxp := xbound1*phosalfa + (xzero1-xbound1)*phosgam;
if ecxp/xzero1 > localp then
begin
ecxp := localp*xzero1/ecxp;
phosgam := phosgam*ecxp;
phosalfa := phosalfa*ecxp;
end
else
if phosgam > localp then phosgam := localp;
comment Local bulk P is set at Kamacite nucleation;
if localp > phosfor then localp := phosfor;
if start = 1 and chp = 0 then
begin
localp := phosfor;
chp := 1;
write(out,<:<10>Kamacite nucleated at temp = :>,
<<dddd.ddd>,t,<:<10>Local bulk Ni,P = :>,
nipcgam,localp,<:<10>:>);
end;
comment undercooling included, start of growth check;
if nipcgam >= bulknipc then
begin
if understart = 0 then understart := t;
if understart - t < undercool then nipcgam := bulknipc;
end;
if nipcgam <= bulknipc and start = 0 then
begin
nipcgam := bulknipc;
start := 0;
end
else
start := 1;
if start = 1 and chp = 0 then aga := 1 else aga := 0;
end
else
begin
nipcgam := nigasat;
nipcalfa := nialsat;
phosgam := phosgasat;
phosalfa := phosalsat;
end;
end of procedure phasedia;
real procedure prescorrect(temp,start,pres);
value temp,start,pres;
real temp,pres;
integer start;
begin
comment Pressure correction to the phase diagram;
prescorrect := pres*1'-11;
end of real procedure prescorrect;
real procedure time(temp,const);
real temp;
real array const;
begin
comment time as a function of temperature
ie the cooling history;
time := -const(1) - temp*const(2);
end of real procedure time;
real procedure temp(time,const);
real time;
real array const;
begin
comment temperature as a function of time
ie cooling history;
temp := (-time-const(1))/const(2);
end of real procedure temp;
real procedure coolrate(temp,const);
real temp;
real array const;
begin
comment this procedure gives the instant cooling rate,
ie the first derivative of temp(time), as a function
of temperature;
coolrate := 1/const(2);
end of real procedure coolrate;
procedure writesection;
begin
comment This procedure writes out all information.;
integer i;
real r1,r2;
mpc := nipc(x2);
write(out,<<-d.dddddddddd'-dd>,
<:<10><10>MPROFILE VER 15<10>:>,
<:<10>Coolingrate = :>,coolrate(t,const),
<:<10>Bulk nipc = :>,bulknipc,
<:<10>P initial = :>,phosfor,
<:<10>Local P in taenite at kamacite nucleation = :>,localp,
<:<10>Undercooling = :>,undercool,
<:<10>Total width of bands :>,xzero1,
<:<10>Resulting gamma bandwidth = :>,
xzero1-xbound1,
<:<10>MPC = :>,mpc,
<:<10>nipcgam = :>,nipcgam,
<:<10>nipcalfa = :>,nipcalfa,
<:<10>phosgam = :>,phosgam,
<:<10>phosalfa = :>,phosalfa,
<:<10>nial0 = :>,nial0,
<:<10>nialsat = :>,nialsat,
<:<10>niga0 = :>,niga0,
<:<10>nigasat = :>,nigasat,
<:<10>phosalsat = :>,phosalsat,
<:<10>phosgasat = :>,phosgasat,
<:<10>temp = :>,t);
end of procedure writesection;
procedure growth(t,xbound,nipc,dx,xzero,err);
real array nipc;
real t,xbound,dx,xzero,err;
begin
comment This procedure solves the differential equation
and calculates the distance which the phase-
boundary will be moved in this timestep.;
dt := abs(time(t,const) - time(t-skt,const));
nipc(1) := nipci(1) := nipcgam;
stable := 1.0;
comment now the diffusion equation is solved;
comment Gamma phase;
for i := 2 step 1 until x2-1 do
begin
r1 := nipc(i-1) + nipc(i+1) - nipc(i)*2;
r2 := stable*dgam(nipc(i))*dt/(dx**2);
if r2 >= 0.5 then
begin
if x2 > 5 then
begin
write(out,<:<10>Too heavy growth in gamma phase<10>:>);
goto bounderr;
end
else
r2 := 0.5;
end;
nipci(i) := nipc(i) + r1*r2;
if nipci(i) > nipci(i-1) and x2 <= 5 then
begin
writesection;
write(out,<:<10>Growth error in gamma at high temp<10>:>);
goto bounderr;
end;
end;
r1 := 2*(nipc(x2-1)-nipc(x2))*stable*
dgam(nipc(x2))*dt/(dx**2);
nipci(x2) := nipc(x2) + r1;
if nipci(x2) > nipci(x2-1) then nipci(x2) := nipci(x2-1);
for i := 2 step 1 until x2 do
nipc(i) := nipci(i);
comment Alpha phase;
if xbound > 0.1*xzero then
begin
nipcal(x2a) := nipci(x2a) := nipcalfa;
if start > 0 and alfainit = 0 then
begin
for i := 1 step 1 until x2a-1 do nipcal(i) := nipcal(x2a);
for i := 1 step 1 until x2a do xa(i) := (i-1)*0.1*xzero/x2a;
alfainit := 1;
end;
for i := x2a-1 step -1 until 2 do
begin
r1 := nipcal(i-1) + nipcal(i+1) - 2*nipcal(i);
r2 := stable*dalfa*dt/(dxal**2);
if r2 > 0.5 then
begin
if x2a > 5 then
begin
write(out,<:<10>Too heavy growth in alpha<10>:>);
goto bounderr;
end
else
r2 := 0.5;
end;
nipci(i) := nipcal(i) + r1*r2;
end;
r1 := 2*(nipcal(2)-nipcal(1))*stable*
dalfa*dt/(dxal**2);
nipci(1) := nipcal(1) + r1;
if nipci(1) > nipci(2) then nipci(1) := nipci(2);
for i := 1 step 1 until x2a-1 do nipcal(i) := nipci(i);
end
else
begin
for i := 1 step 1 until x2a do nipcal(i) := nipcalfa;
end;
comment growth of alpha phase;
grow:
if start > 0 and xbound >= 0.1*xzero then
begin
nipcalint := 0.0;
for i := 1 step 1 until x2a-1 do
nipcalint := nipcalint + dxal*(nipcal(i) + nipcal(i+1))/2;
nipcalint := nipcalint + erral*nipcal(1);
end
else
nipcalint := nipcalfa*xbound;
if start = 0 then goto endgrow;
kk := 0.0;
for j := 2 step 1 until x2 do
kk := kk + nipc(j) + nipc(j-1);
kk := kk*dx/2 +err*nipc(x2);
aa := (nipc(1) - nipc(2))/(2*dx);
bb := nipcal(x2a) - nipc(1);
cc := nipcalint + kk - xzero*bulknipc;
kk := bb**2 - 4*aa*cc;
if nipcal(x2a) - nipc(1) >= 0.0 then
rod := 0.0
else
if aa < 1'-10 then
rod := - cc/bb
else
if kk < 0.0 then rod := dx
else
rod := (-bb - sqrt(kk))/(2*aa);
if rod >= dx then
begin
kk := 0.0;
for j := 3 step 1 until x2 do
kk := kk + nipc(j) + nipc(j-1);
kk := kk*dx/2 + err*nipc(2);
aa := (nipc(2) - nipc(3))/(2*dx);
bb := nipcal(x2a) - nipc(2);
cc := nipcalint + dx*nipcal(x2a) + kk - xzero*bulknipc;
kk := bb**2 - 4*aa*cc;
if aa < 1'-10 then rod := -cc/bb +dx
else
if kk < 0.0 then rod := 2*dx
else
rod := dx + (-bb-sqrt(kk))/(2*aa);
if rod < 2*dx then write(out,<:<10>Enhanced growth at temp = :>,
<<dddd.ddd>,t,<: x2 = :>,<<ddd>,x2);
if rod >= 2*dx then
begin
kk := 0.0;
for j := 4 step 1 until x2 do
kk := kk + nipc(j) + nipc(j-1);
kk := kk*dx/2 + err*nipc(x2);
aa := (nipc(3) - nipc(4))/(2*dx);
bb := nipcal(x2a) - nipc(3);
cc := nipcalint + 2*dx*nipcal(x2a) + kk - xzero*bulknipc;
kk := bb**2 - 4*aa*cc;
if aa < 1'-10 then rod := -cc/bb + 2*dx
else
if kk < 0.0 then rod := 3*dx
else
rod := 2*dx + (-bb-sqrt(kk))/(2*aa);
write(out,<:<10>Heavy enhanced growth at temp = :>,
<<dddd.ddd>,t,<: x2 = :>,<<ddd>,x2);
end;
end;
if rod >= 3*dx then
begin
write(out,<:<10>rod gtr 3 dx at temp = :>,<<-dddd>,t,
<<-d.dd'-d>,<: rod = :>,rod,<: dx = :>,
dx,<: x2 = :>,<<ddd>,x2);
goto bounderr;
end;
xbound := xbound + rod;
if xbound < 0.0 then xbound := 0.0;
endgrow:
end of procedure growth;
procedure move(x,xbound,xzero,nipc,dxg,dx,dxgal,dxal);
real array x,nipc;
real xbound,dxg,dx,dxgal,dxal,xzero;
begin
comment This procedure makes the interpolation
between the old points and generates the new
points.;
comment Gamma phase;
xg(1) := x(1);
x(1) := xbound;
for i := 2 step 1 until x2 do
begin
xg(i) := x(i);
x(i) := x(i-1) + dx;
end;
j := 2;
for i := 1 step 1 until x2 do
begin
for k := i while x(i) > xg(j) and j < oldx2 do j := j + 1;
if j <= oldx2 then
begin
aa := (nipc(j) - nipc(j-1))/dxg;
nipci(i) := aa*(x(i) - xg(j-1)) + nipc(j-1);
end
else
nipci(i) := nipc(oldx2);
end;
for i := 1 step 1 until x2 do nipc(i) := nipci(i);
comment Alpha phase;
if xbound < 0.1*xzero then goto endmove;
for i := 1 step 1 until oldx2a do
xg(i) := xa(i);
for i := 1 step 1 until x2a do
xa(i) := (i-1)*dxal + erral;
j := 2;
for i := 1 step 1 until x2a do
begin
if xa(i) > xg(j) and j < oldx2a then j := j + 1;
if abs(dxgal) < 1'-5 then dxgal := dxal;
if j <= oldx2a then
begin
aa := (nipcal(j) - nipcal(j-1))/dxgal;
nipci(i) := aa*(xa(i) - xg(j-1)) + nipcal(j-1);
end
else
nipci(i) := nipcal(oldx2a);
end;
for i := 1 step 1 until x2a do nipcal(i) := nipci(i);
endmove:
end of procedure move;
comment Headsection starts.;
comment initialization and read in block starts;
err1 := err2 := 0.0;
start := 0;
alfainit := 0;
xbound1 := xbound2 := 0.0;
for i := 1 step 1 until 200 do
begin
xx1(i) := xa(i) := 0.0;
xg(i) := 1.0;
end;
x1 := 1;
understart := 0.0;
agrellgrow := 1;
nialfapre := -1.0;
readstring(in,text,1);
read(in,tempstart,tempslut,maxno,q,const(2),bulknipc,
gamw,alfapre,alfanex,kappa1,kappa2,phosfor,undercool);
tzero := tempstart;
phasestart(phosfor);
nomeas := 0;
const(2) := 1000000/const(2);
comment the length unit is cm and the time unit is years.;
gamw := gamw*'-4;
alfapre := alfapre*'-4;
alfanex := alfanex*'-4;
xzero1 := alfapre*kappa1;
xzero2 := alfanex*kappa2;
xzero1 := 0.5*(xzero1 + xzero2);
for i := 1 step 1 until 200 do
nipc2(i) := nipc(i) := nipcal(i) := bulknipc;
dx1 := xzero1/4;
dx2 := xzero2/4;
dxal := dxgal := xzero1/((x2a-1)*10);
x2 := oldq := x1 := x2a := oldx2 := oldx2a := 5;
const(1) := -tempstart*const(2);
localp := phoschem := phosfor;
chpsat := 0;
chp := 0;
secinyear := 31556926.08;
wni := 58.71;
wfe := 55.847;
skt := -0.5;
comment initialization block end;
for t := tempstart step skt until tempslut do
begin
if time(t,const) > 6'9 then goto now;
phasedia(t,ordn,start,undercool);
if aga = 1 then phasedia(t,ordn,start,undercool);
if nipcalfa < nialfapre then agrellgrow := 0;
nialfapre := nipcalfa;
growth(t,xbound1,nipc,dx1,xzero1,err1);
if xbound1 >= xzero1 then goto bounderr;
x1 := x2;
oldx2a := x2a;
oldx2 := x2;
dxg1 := dx1;
dxg2 := dx2;
dxgal := dxal;
dx1 := sqrt(-skt*dgam(nipc(2))/(coolrate(t,const)*q));
dxal := sqrt(-skt*dalfa/(coolrate(t,const)*q));
if dx1 > (xzero1 - xbound1)/4 then
begin
dx1 := (xzero1 - xbound1)/4;
x2 := 5;
end
else
if dx1*(maxno-1) < xzero1 - xbound1 then
begin
dx1 := (xzero1 - xbound1)/(maxno-1);
x2 := maxno;
end else
x2 := entier((xzero1 - xbound1)/dx1) + 1;
err1 := xzero1 - xbound1 - dx1*(x2-1);
if xbound1 >= 0.1*xzero1 then
begin
if dxal > xbound1/4 then
begin
dxal := xbound1/4.0;
x2a := 5;
end
else
if dxal*(maxno-1) < xbound1 then
begin
dxal := xbound1/(maxno-1);
x2a := maxno;
end
else
x2a := entier(xbound1/dxal) + 1;
end;
erral := xbound1 - dxal*(x2a-1);
if oldq+5 < x2 then
begin
write(out,<:<10>x2 = :>,
<<ddd>,x2,<: and x2a = :>,x2a,<: at temp = :>,<<dddd.ddddd>,t);
oldq := x2;
end;
move(xx1,xbound1,xzero1,nipc,dxg1,dx1,dxgal,dxal);
if tzero - t > 1900.0 then
begin
tzero := t;
writesection;
end;
end of temperature loop;
now:
writesection;
goto slut;
bounderr:
write(out,<:<10>ERROR IN SIMULATION.<10><10>First aid turn q down:>,
<:<10><10>The results beneath are in error.<10>:>);
writesection;
goto slut;
slut:
end of program;
▶EOF◀