Overview
Comment:Replaced ^ with pow (required for -cc and -oct).
Downloads: Tarball | ZIP archive
Timelines: family | ancestors | descendants | both | origin/master | trunk
Files: files | file ages | folders
SHA3-256: 8983165fd09e197b9172bc1d1818f0f259b26ef9a96bb3973026dd4a0cc27fc2
User & Date: geraint@users.sourceforge.net on 2003-09-14 22:31:45.000
Other Links: branch diff | manifest | tags
Context
2003-09-14
23:05:47
Made pi a constant recognised by MTT for -cc and -oct.
Required by NonlinearMSD example.
Will cause problems for models which declare pi as a parameter.
check-in: b36bb9b4a5 user: geraint@users.sourceforge.net tags: origin/master, trunk
22:31:45
Replaced ^ with pow (required for -cc and -oct). check-in: 8983165fd0 user: geraint@users.sourceforge.net tags: origin/master, trunk
2003-09-13
22:26:39
Use std::pow instead of pow to avoid pow double/Complex ambiguity error. check-in: f222810d22 user: geraint@users.sourceforge.net tags: origin/master, trunk
Changes
1
2
3
4
5
6
7
8



9
10
11
12
13
14
15
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18








+
+
+







# Numerical parameter file (SimpleGasTurbine_numpar.txt)
# Generated by MTT at Tue Mar 31 12:15:00 BST 1998

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# %% Version control history
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# %% $Id$
# %% $Log$
# %% Revision 1.1  2000/12/28 16:55:29  peterg
# %% To RCS
# %%
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

#Dummies
 alpha = 1;
 c_v = 1;
 density = 1;
 ideal_gas = 1;
︙
48
49
50
51
52
53
54
55
56


57
58
59
60
61
62
63
64
65
66
67
68

69
70
71
72
73
74
51
52
53
54
55
56
57


58
59
60
61
62
63
64
65
66
67
68
69
70

71
72
73
74
75
76
77







-
-
+
+











-
+






m_c = (p_3*v_c)/(t_3*r);

%Equate pressures
p_4 = p_1;
p_2 = p_3;

%Compute ss temperatures (isentropic)
t_2 = t_1*(p_2/p_1)^alpha;
t_4 = t_3*(p_4/p_3)^alpha;
t_2 = t_1*pow((p_2/p_1),alpha);
t_4 = t_3*pow((p_4/p_3),alpha);

%Find the steady-state work output
w_0 = c_p*(t_3-t_4) - c_p*(t_2-t_1);

%Unit mass flow
mdot = 1;

%Corresponding shaft speed
omega_0 = mdot/k;

%Compute the corresponding load resistance (to absorb that work)
r_l = w_0/(omega_0)^2;
r_l = w_0/pow((omega_0),2);

%Compute shaft inertia to give unit time constant (j_s*r_l)
j_s = r_l;

%Find angular momentum to give shaft speed omega_0
mom_0 =  omega_0*j_s;
1
2
3
4
5
6
7
8



9
10
11
12
13
14
15
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18








+
+
+







# -*-octave-*- Put Emacs into octave-mode
# Input specification (Reactor_input.txt)
# Generated by MTT at Fri Mar  3 11:52:23 GMT 2000
###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## Revision 1.3  2003/06/06 06:38:02  gawthrop
## Made compatible with current MTT.
##
## Revision 1.2  2000/12/28 18:52:24  peterg
## Updated for new formats
##
## Revision 1.1  2000/12/28 17:09:55  peterg
## To RCS
##
###############################################################
︙
29
30
31
32
33
34
35
36

37
38
39
40
41
42
43
44
45
46
47
48
49
50
51

52
53
54
55
56
57
58
32
33
34
35
36
37
38

39
40
41
42
43
44
45
46
47
48
49
50
51
52
53

54
55
56
57
58
59
60
61







-
+














-
+








## Solve for the steady-state concentrations
## Solve for ca - a quadratic.
a 	= k_3*e_3;	#ca^2 
b 	= k_1*e_1 + f_s;	#ca^1 
c 	= -c_0*f_s;

c_a	= (-b + sqrt(b^2 - 4*a*c))/(2*a);
c_a	= (-b + sqrt(pow(b,2) - 4*a*c))/(2*a);

## solve for c_b
c_b 	= c_a*k_1*e_1/(f_s+k_2*e_2);


#States (masses)
x1 = c_a*v_r;
x2 = c_b*v_r;

#Thermal state
x3 = c_p*t_s*v_r;


#Steady-state input q needed to achieve steady-state t_s
q_s = -( (t_0-t_s)*c_p*f_s + e_1*h_1*k_1*x1 + e_2*h_2*k_2*x2 + e_3*h_3*k_3*x1^2);
q_s = -( (t_0-t_s)*c_p*f_s + e_1*h_1*k_1*x1 + e_2*h_2*k_2*x2 + e_3*h_3*k_3*pow(x1,2));

## The two inputs at steady-state
u1 = f_s;
u2 = q_s;


# Set the inputs
︙
1
2
3
4
5
6
7
8



9
10
11
12
13
14
15
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18








+
+
+







# -*-octave-*- Put Emacs into octave-mode
# State specification (Reactor_state.txt)
# Generated by MTT at Fri Mar  3 11:52:23 GMT 2000
###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## Revision 1.3  2003/06/06 06:38:31  gawthrop
## Made compatible with current MTT.
##
## Revision 1.2  2000/12/28 18:52:25  peterg
## Updated for new formats
##
## Revision 1.1  2000/12/28 17:09:55  peterg
## To RCS
##
###############################################################
︙
29
30
31
32
33
34
35
36

37
38
39
40
41
42
43
44
45
46
47
48
49
50
51

52
53
54
55
56
57
58
32
33
34
35
36
37
38

39
40
41
42
43
44
45
46
47
48
49
50
51
52
53

54
55
56
57
58
59
60
61







-
+














-
+








## Solve for the steady-state concentrations
## Solve for ca - a quadratic.
a 	= k_3*e_3;	#ca^2 
b 	= k_1*e_1 + f_s;	#ca^1 
c 	= -c_0*f_s;

c_a	= (-b + sqrt(b^2 - 4*a*c))/(2*a);
c_a	= (-b + sqrt(pow(b,2) - 4*a*c))/(2*a);

## solve for c_b
c_b 	= c_a*k_1*e_1/(f_s+k_2*e_2);


#States (masses)
x1 = c_a*v_r;
x2 = c_b*v_r;

#Thermal state
x3 = c_p*t_s*v_r;


#Steady-state state q needed to achieve steady-state t_s
q_s = -((t_0-t_s)*c_p*f_s + e_1*h_1*k_1*x1 + e_2*h_2*k_2*x2 + e_3*h_3*k_3*x1^2);
q_s = -((t_0-t_s)*c_p*f_s + e_1*h_1*k_1*x1 + e_2*h_2*k_2*x2 + e_3*h_3*k_3*pow(x1,2));

## The two inputs at steady-state
u1 = f_s;
u2 = q_s;

## Load up the states
## Removed by MTT on Thu Dec 28 18:46:20 GMT 2000: mttx(1) = x1;
︙
1
2
3
4
5
6
7
8



9
10
11
12
13
14
15
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18








+
+
+







# -*-octave-*- Put Emacs into octave-mode
# Input specification (ReactorTF_input.txt)
# Generated by MTT at Fri Mar  3 11:52:23 GMT 2000
###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## Revision 1.2  2003/06/06 06:38:44  gawthrop
## Made compatible with current MTT.
##
## Revision 1.1  2000/12/28 17:12:57  peterg
## To RCS
##
###############################################################

## Reduce steady-state parameter file (ReactorTF_sspar.r)
## as siso_sspar ecxept that inputs/states have different meaning
︙
26
27
28
29
30
31
32
33

34
35
36
37
38
39
40
41
42
43
44
45
46
47
48

49
50
51
52
29
30
31
32
33
34
35

36
37
38
39
40
41
42
43
44
45
46
47
48
49
50

51
52
53
54
55







-
+














-
+





## Solve for the steady-state concentrations
## Solve for ca - a quadratic.
a 	= k_3*e_3;	#ca^2 
b 	= k_1*e_1 + f_s;	#ca^1 
c 	= -c_0*f_s;

c_a	= (-b + sqrt(b^2 - 4*a*c))/(2*a);
c_a	= (-b + sqrt(pow(b,2) - 4*a*c))/(2*a);

## solve for c_b
c_b 	= c_a*k_1*e_1/(f_s+k_2*e_2);


#States (masses)
x1 = c_a*v_r;
x2 = c_b*v_r;

#Thermal state
#x3 = c_p*t_s*v_r;


#Steady-state input q needed to achieve steady-state t_s
q_s = -( (t_0-t_s)*c_p*f_s + e_1*h_1*k_1*x1 + e_2*h_2*k_2*x2 + e_3*h_3*k_3*x1^2);
q_s = -( (t_0-t_s)*c_p*f_s + e_1*h_1*k_1*x1 + e_2*h_2*k_2*x2 + e_3*h_3*k_3*pow(x1,2));

# Set the inputs
## Removed by MTT on Thu Jun  5 14:13:24 BST 2003: mttu(1) = q_s + 0.1*q_s*(t>0.01); # q (ReactorTF)
reactortf__t	=  q_s + 0.1*q_s*(t>0.01); # q (ReactorTF)
1
2
3
4
5
6
7
8



9
10
11
12
13
14
15
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18








+
+
+







# -*-octave-*- Put Emacs into octave-mode
# State specification (ReactorTF_state.txt)
# Generated by MTT at Fri Mar  3 11:52:23 GMT 2000
###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## Revision 1.2  2003/06/06 06:39:05  gawthrop
## Made compatible with current MTT.
##
## Revision 1.1  2000/12/28 17:12:57  peterg
## To RCS
##
###############################################################

## Reduce steady-state parameter file (ReactorTF_sspar.r)
## as siso_sspar ecxept that states/states have different meaning
︙
26
27
28
29
30
31
32
33

34
35
36
37
38
39
40
29
30
31
32
33
34
35

36
37
38
39
40
41
42
43







-
+








## Solve for the steady-state concentrations
## Solve for ca - a quadratic.
a 	= k_3*e_3;	#ca^2 
b 	= k_1*e_1 + f_s;	#ca^1 
c 	= -c_0*f_s;

c_a	= (-b + sqrt(b^2 - 4*a*c))/(2*a);
c_a	= (-b + sqrt(pow(b,2) - 4*a*c))/(2*a);

## solve for c_b
c_b 	= c_a*k_1*e_1/(f_s+k_2*e_2);


#States (masses)
x1 = c_a*v_r;
︙
1
2
3
4
5
6
7
8



9
10
11
12
13
14
15
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18








+
+
+







# -*-octave-*- Put Emacs into octave-mode
# Input specification (ReactorTQ_input.txt)
# Generated by MTT at Fri Mar  3 11:52:23 GMT 2000
###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## Revision 1.2  2003/06/06 06:39:20  gawthrop
## Made compatible with current MTT.
##
## Revision 1.1  2000/12/28 17:19:08  peterg
## To RCS
##
###############################################################

## Reduce steady-state parameter file (ReactorTQ_sspar.r)
## as siso_sspar ecxept that inputs/states have different meaning
︙
26
27
28
29
30
31
32
33

34
35
36
37
38
39
40
29
30
31
32
33
34
35

36
37
38
39
40
41
42
43







-
+








## Solve for the steady-state concentrations
## Solve for ca - a quadratic.
a 	= k_3*e_3;	#ca^2 
b 	= k_1*e_1 + f_s;	#ca^1 
c 	= -c_0*f_s;

c_a	= (-b + sqrt(b^2 - 4*a*c))/(2*a);
c_a	= (-b + sqrt(pow(b,2) - 4*a*c))/(2*a);

## solve for c_b
c_b 	= c_a*k_1*e_1/(f_s+k_2*e_2);


#States (masses)
x1 = c_a*v_r;
︙
1
2
3
4
5
6
7
8



9
10
11
12
13
14
15
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18








+
+
+







# -*-octave-*- Put Emacs into octave-mode
# State specification (ReactorTQ_state.txt)
# Generated by MTT at Fri Mar  3 11:52:23 GMT 2000
###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## Revision 1.2  2003/06/06 06:39:39  gawthrop
## Made compatible with current MTT.
##
## Revision 1.1  2000/12/28 17:19:08  peterg
## To RCS
##
###############################################################

## Reduce steady-state parameter file (ReactorTQ_sspar.r)
## as siso_sspar ecxept that states/states have different meaning
︙
26
27
28
29
30
31
32
33

34
35
36
37
38
39
40
29
30
31
32
33
34
35

36
37
38
39
40
41
42
43







-
+








## Solve for the steady-state concentrations
## Solve for ca - a quadratic.
a 	= k_3*e_3;	#ca^2 
b 	= k_1*e_1 + f_s;	#ca^1 
c 	= -c_0*f_s;

c_a	= (-b + sqrt(b^2 - 4*a*c))/(2*a);
c_a	= (-b + sqrt(pow(b,2) - 4*a*c))/(2*a);

## solve for c_b
c_b 	= c_a*k_1*e_1/(f_s+k_2*e_2);


#States (masses)
x1 = c_a*v_r;
︙
1
2
3
4
5
6
7
8
9



10
11
12
13
14
15
16
17
18
19
20

21
22
23
24
25
26
27
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

23
24
25
26
27
28
29
30









+
+
+










-
+







# -*-octave-*- Put Emacs into octave-mode
# Numerical parameter file (PPPCantileverBeam_numpar.txt)
# Generated by MTT at Mon Apr 19 06:24:08 BST 1999

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# %% Version control history
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# %% $Id$
# %% $Log$
# %% Revision 1.1  2000/12/28 17:27:26  peterg
# %% To RCS
# %%
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

# Parameters
N = 16;
BeamLength = 0.58;
BeamWidth = 0.05;
BeamThickness = 0.005;
Youngs = 1e6;
Density = 1e5;
Area = BeamWidth*BeamThickness;
AreaMoment = (BeamThickness*BeamWidth^2)/12;
AreaMoment = (BeamThickness*pow(BeamWidth,2))/12;


EI= 58.6957			# from Reza
rhoA= 0.7989			# from Reza
 
dz = BeamLength/N;		# Incremental length
dm = rhoA*dz;			# Incremental mass
︙
1
2
3
4
5
6
7
8
9



10
11
12
13
14
15
16
17
18
19
20

21
22
23
24
25
26
27
28
29
30

31
32
33
34
35
36
37
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

23
24
25
26
27
28
29
30
31
32

33
34
35
36
37
38
39
40









+
+
+










-
+









-
+







# -*-octave-*- Put Emacs into octave-mode
# Numerical parameter file (CantileverBeam_numpar.txt)
# Generated by MTT at Mon Apr 19 06:24:08 BST 1999

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# %% Version control history
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# %% $Id$
# %% $Log$
# %% Revision 1.1  2000/12/28 17:58:27  peterg
# %% To RCS
# %%
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

# Parameters
N = 21;
BeamLength = 0.58;
BeamWidth = 0.05;
BeamThickness = 0.005;
Youngs = 68.94e9;
Density =  2712.8;
Area = BeamWidth*BeamThickness;
AreaMoment = (BeamWidth*BeamThickness^3)/12;
AreaMoment = (BeamWidth*pow(BeamThickness,3))/12;

EI = Youngs*AreaMoment;
rhoA = Density*Area;

dz = BeamLength/N;		# Incremental length
dm = rhoA*dz;			# Incremental mass
dk = EI/dz;			# Incremental stiffness
dr = 0;				# Damping

K = sqrt(EI/rhoA)/BeamLength^2;	# Normalising factor
K = sqrt(EI/rhoA)/pow(BeamLength,2);	# Normalising factor


# EI= 58.6957			# from Reza
# rhoA= 0.7989			# from Reza
 


︙
1
2
3
4
5
6
7
8
9



10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

26
27
28
29
30
31
32
33
34
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

28
29
30
31
32
33
34
35
36
37









+
+
+















-
+









# -*-octave-*- Put Emacs into octave-mode
# Numerical parameter file (pPinnedBeam_numpar.txt)
# Generated by MTT at Mon Apr 19 06:24:08 BST 1999

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# %% Version control history
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# %% $Id$
# %% $Log$
# %% Revision 1.2  2003/06/11 16:03:06  gawthrop
# %% Updated examples for latest MTT.
# %%
# %% Revision 1.1  2000/12/28 17:59:05  peterg
# %% To RCS
# %%
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

## Number of lumps
Lumps = 20;			# Number of lumps

## Beam physical parameters
BeamLength = 0.60;
BeamWidth  = 0.05;
BeamThickness = 0.003;
Youngs = 68.94e9;
Density =  2712.8;
Area = BeamWidth*BeamThickness;
AreaMoment = (BeamWidth*BeamThickness^3)/12;
AreaMoment = (BeamWidth*pow(BeamThickness,3))/12;
EI = Youngs*AreaMoment;
rhoA = Density*Area;

## Segments
dz = BeamLength/Lumps;	        # Incremental length
dm = rhoA*dz;			# Incremental mass
dk = EI/dz;			# Incremental stiffness
dr = 0;				# Damping

1
2
3
4
5
6
7
8
9
10
11

12
1
2
3
4
5
6
7
8

9

10
11








-

-
+

## -*-octave-*- Put Emacs into octave-mode ##
 
## 
## System NonlinearMSD, representation input, language txt; 
## File NonlinearMSD_input.txt; 
## Generated by MTT on Thu Mar  7 10:50:46 GMT 2002; 

## First term is the equilibrium input; last term is the perturbation input.
## Removed by MTT on Tue Jun 10 16:50:53 BST 2003: NonlinearMSD_yu	= k*(l^2)*cos(eta/2)*2*(sin(eta/2)-sin(alpha/2)) + 1e-2;

nonlinearmsd__yu = k*(l^2)*cos(eta/2)*2*(sin(eta/2)-sin(alpha/2)) + 1e-2;
nonlinearmsd__yu = k*(pow(l,2))*cos(eta/2)*2*(sin(eta/2)-sin(alpha/2)) + 1e-2;

1
2
3
4
5
6
7
8



9
10
11
12
13
14
15
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18








+
+
+







# Numerical parameter file (SimpleGasTurbine_numpar.txt)
# Generated by MTT at Tue Mar 31 12:15:00 BST 1998

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# %% Version control history
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# %% $Id$
# %% $Log$
# %% Revision 1.1  2000/12/28 18:08:28  peterg
# %% To RCS
# %%
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

# Parameters
c_p = 	1005.0;
c_v = 	718.0; 
gamma_0 =  c_p/c_v;
alpha = (gamma_0-1)/gamma_0;
︙
30
31
32
33
34
35
36
37
38


39
40
41
42
43
44
45
46
47
48
49
50

51
52
53
54
55
56
33
34
35
36
37
38
39


40
41
42
43
44
45
46
47
48
49
50
51
52

53
54
55
56
57
58
59







-
-
+
+











-
+






m_c = (p_3*v_c)/(t_3*r);

%Equate pressures
p_4 = p_1;
p_2 = p_3;

%Compute ss temperatures (isentropic)
t_2 = t_1*(p_2/p_1)^alpha;
t_4 = t_3*(p_4/p_3)^alpha;
t_2 = t_1*pow((p_2/p_1),alpha);
t_4 = t_3*pow((p_4/p_3),alpha);

%Find the steady-state work output
w_0 = c_p*(t_3-t_4) - c_p*(t_2-t_1);

%Unit mass flow
mdot = 1;

%Corresponding shaft speed
omega_0 = mdot/k;

%Compute the corresponding load resistance (to absorb that work)
r_l = w_0/(omega_0)^2;
r_l = w_0/pow((omega_0),2);

%Compute shaft inertia to give unit time constant (j_s*r_l)
j_s = r_l;

%Find angular momentum to give shaft speed omega_0
mom_0 =  omega_0*j_s;
1
2
3
4
5
6
7
8



9
10
11
12
13
14
15
16
17
18
19
20
21


22
23
24
25
26
27

28
29
30
31
32
33
34
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22


23
24
25
26
27
28
29

30
31
32
33
34
35
36
37








+
+
+











-
-
+
+





-
+







% Steady-state parameter file (SimpleGasTurbine_sspar.r)
% Generated by MTT at Thu Mar 26 16:28:59 GMT 1998

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % Version control history
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % $Id$
% % $Log$
% % Revision 1.1  2000/12/28 18:08:28  peterg
% % To RCS
% %
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Find stored mass to give combustion chamber pressure p_3 (at
% temperature t_3
m_c := (p_3*v_c)/(t_3*r);

%Equate pressures
p_4 := p_1;
p_2 := p_3;

%Compute ss temperatures (isentropic)
t_2 := t_1*(p_2/p_1)^alpha;
t_4 := t_3*(p_4/p_3)^alpha;
t_2 := t_1*pow((p_2/p_1),alpha);
t_4 := t_3*pow((p_4/p_3),alpha);

%Find the steady-state work output
w_0 := c_p*(t_3-t_4) - c_p*(t_2-t_1);

%Compute the corresponding load resistance (to absorb that work)
r_l := w_0/(omega_0)^2;
r_l := w_0/pow((omega_0),2);

%Unit mass flow
mdot := 1;

%Corresponding shaft speed
omega_0 := mdot/k;

︙

MTT: Model Transformation Tools
[ GitHub | SourceHut | Sourceforge | Fossil RSS ]