plot2d of romberg-integrating functions?
- Subject: plot2d of romberg-integrating functions?
- From: Robert Gloeckner
- Date: Fri, 16 Dec 2005 00:35:50 +0100
hello,
here is the full code. i can plot most functions.
but i have problems plotting
n_RXI(t, TRATE)
nRTXI(T, TRATE)
in 2d and in 3d ...
so, i would be very happy, if you could give me hints on that, too ;)
thanks,
Robert
/* $Header:
/home/cvs/SigmApi/examples/Api/Crystallization/Moneke_maxima.txt,v 1.6
2005/12/15 23:18:19 robert Exp $ */
/******************************************************************************/
/* 2005-12-16
*/
/* Robert Gloeckner, r134gloeck4514ner@dki.4tu-darm45stadt.de (no
numbers) */
/*
*/
/* these equations calculate crystallisation degree of a cooling melt of
*/
/* polypropylene. the melt initially is 170°C (T_Max) and is cooled down
to */
/* 20°C (T_Min) with a cooling speed of T_Rate. there are also
calculations */
/* for variable cooling-rates.
*/
/* equations and constants were published by m. moneke in his phd 2000
*/
/* original derivation by kolmogorov et al.
*/
/*
*/
/******************************************************************************/
killall();
/**************** material constants **************************/
NTEMP : 273.15;
R : 1.0;
p : 30.0;
T0m : 187.0 + NTEMP; /* . */
T0g : -3.0 + NTEMP; /* . */
G0 : 2.8e3; /* . */
C : 265; /* needed for eq4_11, van Krevelen-factor */
C1 : 1443.3; /* . */
C2 : 51.6; /* . */
Kg : 3.0e5; /* . */
ag : 0.03; /* . */
am : 0.0205; /* . */
AN : 2e18; /* . */
BN : 15.0; /* . */
TC : 57.0 + NTEMP; /* . */
Cd : 5.0; /* eq4_10 (Arrhenius) */
MGF : 0.0; /* 0.310; #PP: enthaelt auch Russanteil */
MC_inf : 0.47; /* . :: c_max */
cp_a1 : 2.13; /* . */
cp_a2 : 3.27e-3; /* . */
cp_c1 : 0.966; /* . */
cp_c2 : 1.354e-2; /* . */
C_c : 0.0894; /* . is taitc */
v_cryst_1 : 1.0343E-3; /* . */
v_cryst_2 : 0; /* . */
v_cryst_3 : 0; /* . */
v_cryst_4 : 8912; /* . */
v_cryst_5 : 0; /* . */
v_amorph_1 : 1.1605E-3;/* . */
v_amorph_2 : 8.26e-7; /* . */
v_amorph_3 : 1355; /* . */
v_amorph_4 : 4.01e-3; /* . */
lamda_cryst_1 : 0.54; /* . */
lamda_amorph_1 : 0.121; /* . */
lamda_amorph_2 : 1.25e-5; /* . */
lamda_amorph_3 : 1.88e-4; /* . */
lamda_amorph_4 : 7e-8; /* . */
cp_GF : 795; /* . */
vGF : 1/2540.0; /* . -> glass rho */
lamda_GF : 1.035; /* . */
/**************** heat capacity *******************************/
MC(T, _p) := 0.468; /* ist eine Schaetzung, eigentlich Tabelle
B.4 -> Funktion!!! */
cp_amorph(T) := cp_a1 + cp_a2 * T;
cp_cryst(T) := cp_c1 + cp_c2 * T;
cp_semi(T, _p) := MC(T, _p) * cp_cryst(T) + (1 - MC(T, _p)) *
cp_amorph(T);
eq5_4(T, _p) := MGF * cp_GF + (1 - MGF) * cp_semi(T, MC(T, _p));
/***************** specific volume ****************************/
B_cryst(T) := v_cryst_4 * exp( -v_cryst_5 * T);
v_cryst_0(T) := v_cryst_1 + v_cryst_2 * T + v_cryst_3 * T * T;
v_cryst(T, _p) := v_cryst_0(T) * (1 - C_c * log(1 + (_p / B_cryst(T)
) ) );
B_amorph(T) := v_amorph_3 * exp( -v_amorph_4 * T);
v_amorph_0(T) := v_amorph_1 + v_amorph_2 * T;
v_amorph(T, _p) := v_amorph_0(T) * (1 - C_c * log(1 + (_p /
B_amorph(T) ) ) );
v_semi(T, _p) := MC(T, _p) * v_cryst(T, _p) + (1 - MC(T, _p)) *
v_amorph(T, _p);
v(T, _p) := MGF * vGF + (1 - MGF) * v_semi(T, _p);
/***************** heat conduction ****************************/
lamda_cryst(T, _p) := lamda_cryst_1;
lamda_amorph(T, _p) := lamda_amorph_1 + lamda_amorph_2 * _p +
(lamda_amorph_3 + lamda_amorph_4 * _p) * T;
V_C(T, _p) := MC(T, _p) * v_cryst(T) / v_semi(T, _p);
V_GF(T, _p) := MGF * vGF / v_semi(T, _p);
lamda_semi(T, _p) := 1 / ( ( V_C(T, _p) / lamda_cryst(T, _p) )
+ ( ( 1 - V_C(T, _p) ) / lamda_amorph(T, _p) ) );
lamda(T, _p) := 1 / ( (V_GF(T, _p) / lamda_GF) + ( ( 1
- V_GF(T, _p) ) / lamda_semi(T, _p) ) );
/************* nucleation rate dN/dt **************************/
eq4_7(T, _p) := ((AN) / (sqrt(2 * %pi * (BN^2)))) * (exp( (-(T -
TC)^2) / (2 * (BN^2)) ) );
/************ nucleation grow-rate ****************************/
Tmp(_p) := T0m + _p * am;
Tgp(_p) := T0g + _p * ag;
f(T) := (2 * T) / (Tmp(p) + T);
eq4_12a(T, _p) := (2 * T) / (Tmp(_p) + T);
eq4_12(T, _p) := (Kg) / (eq4_12a(T, _p) * T * (Tmp(_p) - T));
eq4_11(T, _p) := ( C * (Tmp(_p)^2) ) / ( (T^2) * (Tmp(_p)
- T) );
eq4_10(T, _p) := (Cd * (Tmp(_p)^2) ) / ( (T^2) * (Tmp(_p)
- T) );
eq4_9(T, _p) := (C1) / (R * (C2 + T - Tgp(_p)));
eq4_8__4_9_4_11(T, _p) := G0 * exp(-eq4_9(T, _p) - eq4_11(T, _p) );
/* Vogel-Fulcher */
eq4_8__4_10_4_11(T, _p) := G0 * exp(-eq4_10(T, _p) - eq4_11(T, _p) );
/* Arrhenius */
eq4_8__4_9_4_12(T, _p) := G0 * exp(-eq4_9(T, _p) - eq4_12(T, _p) );
/* Hoffmann-Lauritzen */
/**************** graphing borders ****************************/
T_Min : 20 ;
T_Max : 170 ;
T_Rate: -40;
p_Min : 1;
p_Max : 220;
/*************** contant cooling rate definition **************/
T_(t) := T_Max + NTEMP + t * T_Rate; /* linear cooling t[s] */
time(T) := ( T - T_Max ) / T_Rate; /* for time-Temp conversion */
t_Min : float( min( time(T_Min), time(T_Max)) );
t_Max : float( max( time(T_Min), time(T_Max)) );
/************** constant cooling rate: 1D-functions ***********/
G(T) := float(eq4_8__4_9_4_12(T, p_Min)); /* grow rate of nucleus at a
constant pressure */
A(T) := float(eq4_7(T, p_Min)); /* creation of nucleus at a
constant pressure */
Q1(TT) := romberg( G(T1), T1, T_Min, TT);
Q2(TT) := romberg( A(T1) * ( Q1(T1)^3 ), T1, T_Min, TT);
/* this is the symbolic kolmogorov equation */
SG(t, z) := integrate( G(T_(tau)), tau, z, t); /* inner integral */
SPHI(t) := integrate( A(T_(z)) * SG(t, z)^3 , z, t_Min, t); /* outer
integral */
XI(t) := 1 - exp( (-4 * %pi / 3) * SPHI(t)); /* degree of
crystallinity at time t [s] */
/* numerical integration method: romberg */
rombergtol : 1E-1; /* very low accuracy */
n_SG(t, z) := romberg( G(T_(tau)), tau, z, t);
n_SPHI(t) := romberg( A(T_(z)) * n_SG(t, z)^3 , z, t_Min, t);
n_XI(t) := 1 - exp( (-4 * %pi / 3) * n_SPHI(t)); /* degree of
crystallinity at time t [s] */
n_T_XI(T) := n_XI(float(time(T))); /* degree of crystallinity over T */
compile(all);
/**************************************************************/
/****** temperature-(cooling)-rate dependent functions ********/
/**************************************************************/
RT(t, TRATE) := T_Max + NTEMP + t * TRATE; compile(RT); /*
Temperature(time, dT/dt */
Rtime(T, TRATE) := ( T - T_Max ) / TRATE; compile(Rtime); /*
time(Temperature, dT/dt */
timemax(TRATE) := float( max( Rtime(T_Min, TRATE), Rtime(T_Max,
TRATE))); compile(timemax);
timemin(TRATE) := float( min( Rtime(T_Min, TRATE), Rtime(T_Max,
TRATE))); compile(timemin);
/* minimum and maximum of time( dT/dt) */
n_RSG(t, z, TRATE) := romberg( G(RT(tau, TRATE)), tau, z, t);
compile(n_RSG); /* sum of nucleation growing-function */
n_RSPHI(t, TRATE) := romberg( A(RT(z, TRATE)) * n_RSG(t, z, TRATE)^3 ,
z, timemin(TRATE), t); compile(n_RSPHI); /* double-sum (nucleation
creation * sum of nucleation grow^3 */
n_RXI(t, TRATE) := 1 - exp( (-4 * %pi / 3) * n_RSPHI(t, TRATE));
compile(n_RXI); /* crystallisation degree dependent on time and dT/dt */
nRTXI(T, TRATE) := n_RXI(float(Rtime(T, TRATE)), TRATE); compile(nRTXI);
/* crystallisation degree as a function of T and dT/dt */
maxXI(TRATE) := n_RXI(float(timemax(TRATE)), TRATE); compile(maxXI); /*
maximum crystallisation degree as a function of dT/dt (implicit
depending on T_Min and T_Max also) */
/******************* testing & displaying **********************/
G(T_(0)); /* constant cooling rate: testing functions */
G(T_(1));
A(T_(0));
A(T_(1));
n_SG(0, 1), numer; /* constant cooling rate: testing numerical
integral-calculation */
n_SG(0, 2), numer;
n_SPHI(0), numer;
n_SPHI(1), numer;
n_XI(0), numer;
n_XI(1), numer;
plot2d(T_, [t, t_Min, t_Max]); /* display Temperature-rante */
plot2d(G(T_(t)), [t, t_Min, t_Max]); /* constant cooling rate:
nucleation grow */
plot2d(A(T_(t)), [t, t_Min, t_Max]); /* constant cooling rate:
nucleation creation rate */
plot2d(n_XI, [t, t_Min, t_Max], [nticks, 10]); /* constant cooling rate:
crystall-degree(t) */
plot2d(n_T_XI, [T, T_Min, T_Max], [nticks, 5]); /* constant cooling
rate: crystall-degree(T) */
plot2d(maxXI, [TRATE, -100, -1], [nticks, 10]); /* maximum
crystall-degree as a function of cooling rate */
--
https://gna.org/projects/mipisti - (microscope) picture stitching
T_a_k_e__c_a_r_e__o_f__y_o_u_r__R_I_G_H_T_S.
P_r_e_v_e_n_t__L_O_G_I_C--P_A_T_E_N_T_S
http://www.ffii.org, http://www.nosoftwarepatents.org