Ahmed Siddique
$6/hr
Data Scraping and Analysis,scientific code conversion, Chemical engineering
- Reply rate:
- -
- Availability:
- Full-time (40 hrs/wk)
- Age:
- 32 years old
- Location:
- Chakwal, Punjab, Pakistan
- Experience:
- 4 years
Fortran-MATLAB_conversion
%TRACING OF THE THERMODYNAMIC VARIABLES STATION-BY-STATION
%%REQUIRED INPUTS
%==========================================================
% INPUT REQUIRED IN MATRIX FORM OF (NSTATN,NSTRM) DIMENSIONS
% PSI,CZ,CR,RCU,DENSITY,RADIUS,HTOTAL,KOUNT,
% PTOTAL,RHS,ERRRHS,ERRDENS,OMEGLOS,ENTROPY
NSTATN = 2;
NSTRM = 2;
% ?????????????????????????????????????????????
% INPUTS FOR BLOCK 1
RCU = [1, 1; 1, 1];
RADIUS = [1, 1; 1, 1];
CZ = [1, 1; 1, 1];
CR = [1, 1; 1, 1];
HTOTAL = [1, 1; 1, 1];
PTOTAL = [1, 1; 1, 1];
DENSITY = [1,1; 1, 1]; % NOT REQUIRED BUT INTILIZED FOR SPEED AND TO AVOID SYNTEX ERROR
%???????????????????????????????????????????????
% ?????????????????????????????????????????????
% INPUTS FOR BLOCK 2
% BLOCK 2 SECTION 1
PSI = [1, 1; 1, 1];
% BLOCK 2 SECTION 1
ENTROPY = [1, 1; 1, 1];
%???????????????????????????????????????????????
% ?????????????????????????????????????????????
% INPUTS FOR BLOCK 3
% BLOCK 3 SECTION 2
RHS = [1, 1; 1, 1];
%
%???????????????????????????????????????????????
%==========================================================
% INPUT REQRIED AS SINGLE VALUE
% NSTATN, NSTRM, ERRDENS, H,H2,TWOH,ZERO,ONE,TWO,XMASS,TWOPI,NSWEEP,NMIN,
% ERRLIMP,ERRLIMR,ERRLIMD,CP,TWOCP,RGAS,OMEGA,GAMA,GAMAM,NLE,NTE
% ??????????????????????????????????????????????
% INPUTS FOR BLOCK 1
TWO = 2;
GAMAM = 2;
RGAS = 2;
HSTATIC = 2;
CP = 2;
% ???????????????????????????????????????????????
% ??????????????????????????????????????????????
% INPUTS FOR BLOCK 2
ZERO = 0;
% BLOCK 2 SECTION 1
NLE = 2;
NTE = 2;
% BLOCK 2 SECTION 3
OMEGA = 2;
HSTATIC = 2;
CP = 2;
% ???????????????????????????????????????????????
% ??????????????????????????????????????????????
% INPUTS FOR BLOCK 3
% BLOCK 3 SECTION 1
SPECIFIED_PRESSURE_RATIO = 1;
% BLOCK 3 SECTION 2
XMASS = 2;
TWOPI = 2;
TWOH = 2;
% ???????????????????????????????????????????????
%%START OF THE EXECUTION PART
%-----------------------------------------------------------
%!!!!!!!!!!!! BLOCK 1 !!!!!!!!!!!!!!!!
% UPDATING DENSITY TO THE INLET CONDITIONS
for j = 1:NSTRM;
CSQ= (RCU(1,J)/RADIUS(1,J))^2 + CZ(1,J)^2 + CR(1,J)^2;
HSTATIC = HTOTAL(1,J) - CSQ/TWO;
PSTATIC = PTOTAL(1,J) * (HSTATIC/HTOTAL(1,J))^GAMAM;
DENSITY(1,J) = PSTATIC/(RGAS*HSTATIC/CP);
end
%!!!!!!!!!!!! BLOCK 2 !!!!!!!!!!!!!!!!
% SWEEP THE PLANES FROM PLANE 2 TO NSTATN
ERRDENS = ZERO;
CHANGE = ZERO;
%__________________SECTION 1_________________
for I = 2:NSTATN;
% TRACING EACH STREAMLINE BACK TO THE PREVIOUS STATION
% FOR A DUCT (REFERENCE PLANE)
% OR TO THE LEADING EDGE IN CASE OF A BLADE ROW
LEFT = I-1;
if (I > NLE && I <= NTE)
LEFT = NLE;
end
NSTART = 1;
PSIDN = PSI(LEFT,NSTART);
PSIUP = PSI(LEFT,NSTART+1);
% __________________SECTION 2_________________
% SWEEP FROM HUB TO SHROUD AT EACH REFERENCE PLANE
% TO LOCATE ORIGIN OF STREAMLINE
for J = 1, NSTRM;
DESIRED = PSI(I,J);
while(DESIRED > PSIUP | DESIRED < PSIDN)
% IN FORTRAN CODE IT WAS STATED THAT IF PRESSURE IS WITHIN LIMITS SKIP THE
% FOLLOWING 13 LINES TILL DELTA CALCULATIONS, HERE LOGIC IS MODIFIFED
% I.E. EXECUTE THE LINES IF ANY OF PSIUP OR PSIDN IS NOT WITHIN LIMITS
NSTART = NSTART + 1;
PSIDN = PSI(LEFT,NSTART);
PSIUP = PSI(LEFT,NSTART+1);
if(NSTART == NSTRM)
fprintf('PSI VALUE = %f STATION I = %f STREAMLINE J = %f PSI UP = %f PSI DOWN = %f \n',...
DESIRED,I,J,PSIUP,PSIDN)
fprintf('CANNOT TRACE STREAMLINE')
% ERROR TERMINATES THE PROGRAM
return
end
end
%________________SECTION 3_________________________
% STREAMLINE ORIGIN HAS BEEN LOCATED
% ----------------------------------
DELTA = (DESIRED - PSIDN) / (PSIUP - PSIDN);
% DETERMINE WHETHER THIS A ROTATING BLADE
% THEN CALCULATE ROTHALPY AT PREVIOUS STATION
% IF NOT CALCULATE TOTAL ENTHALPY
% -------------------------------------------
ROTATE = ZERO;
if(I > NLE && I <= NTE)
ROTATE = OMEGA;
end
% QUANTITIES AT "1" NEEDED FOR CONSERVATION PRINCIPLE
% ---------------------------------------------------
RCU1 = DELTA * (RCU(LEFT,NSTART+1) - RCU(LEFT,NSTART))+...
RCU(LEFT,NSTART);
HTOTAL1 = DELTA * (HTOTAL(LEFT,NSTART+1) - HTOTAL(LEFT,NSTART))...
+ HTOTAL(LEFT,NSTART);
PTOTAL1 = DELTA * (PTOTAL(LEFT,NSTART+1) - PTOTAL(LEFT,NSTART))...
+ PTOTAL(LEFT,NSTART);
RAD1 = DELTA * (RADIUS(LEFT,NSTART+1) - RADIUS(LEFT,NSTART))...
+ RADIUS(LEFT,NSTART);
% QUANTITIES AT "1" NEEDED FOR LOSS CALCULATION
% ---------------------------------------------
DENS1 = DELTA*(DENSITY(LEFT,NSTART+1)-DENSITY(LEFT,NSTART))...
+ DENSITY(LEFT,NSTART);
CZ1 = DELTA*(CZ(LEFT,NSTART+1) - CZ(LEFT,NSTART))...
+ CZ(LEFT,NSTART);
CR1 = DELTA*(CR(LEFT,NSTART+1) - CR(LEFT,NSTART))...
+ CR(LEFT,NSTART);
ENTROP1 = DELTA*(ENTROPY(LEFT,NSTART+1)-ENTROPY(LEFT,NSTART))...
+ ENTROPY(LEFT,NSTART);
C1SQ = CZ1^2 + CR1^2 + (RCU1/RAD1)^2;
HSTAT1 = HTOTAL1 - C1SQ / TWO;
PSTAT1 = PTOTAL1 * (HSTAT1/HTOTAL1)^GAMAM;
% ROTATING AND NON-ROTATING QUANTITIES AT REFERENCE STATION
% ---------------------------------------------------------
ROTALP1 = HTOTAL1 - ROTATE * RCU1;
HOR2 = ROTALP1 + (ROTATE * RADIUS(I,J))^2 / TWO;
HOR1 = ROTALP1 + (ROTATE * RAD1)^2 / TWO;
POR1 = PTOTAL1 * (HOR1/HTOTAL1)^GAMAM;
POR2IDL = POR1 *(HOR2/HOR1 )^GAMAM;
if(I > NLE && I <= NTE)
OMEGLOS = 0.03*(I- NLE)/(NTE-NLE);
PLOSS = OMEGLOS * ( POR1 - PSTAT1);
else
OMEGLOS = 0.0;
PLOSS = 0.0;
end
POR2 = POR2IDL - PLOSS;
%!!!!!!!!!!!! BLOCK 3 !!!!!!!!!!!!!!!!
% ________________SECTION 1____________________________
% PROJECTS 1 AND 2 : ROTOR WITH "RCU" SPECIFIED
% -------------------------------------------------------
RCU(I,J) = RCU1;
if(I == NLE+1)
RCU(I,J) = 117.85;
end
if(I == NLE+2)
RCU(I,J) = 157.1;
end
if(I == NLE+3)
RCU(I,J) = 196.35;
end
if(I == NTE)
RCU(I,J) = 235.6;
end
HTOTAL(I,J) = HTOTAL1 + ROTATE * (RCU(I,J) - RCU1);
PTOTAL(I,J) = POR2 * (HTOTAL(I,J)/HOR2)^GAMAM;
% PROJECT 3 : ROTOR WITH "PO2/PO1" SPECIFIED
% -------------------------------------------------------
PTOTAL(I,J) = PTOTAL1 * SPECIFIED_PRESSURE_RATIO;
HTOTAL(I,J) = HOR2 * (PTOTAL(I,J)/POR2)^GAMAM;
RCU(I,J) = RCU1 + (HTOTAL(I,J) - HTOTAL1)/ROTATE;
% ________________SECTION 2____________________________
% COMMON CALCULATION BLOCK
% ------------------------
CU = RCU(I,J) / RADIUS(I,J);
VU = CU - ROTATE * RADIUS(I,J);
V2SQ = VU^2 + CZ(I,J)^2 + CR(I,J)^2;
C2SQ = CU^2 + CZ(I,J)^2 + CR(I,J)^2;
HSTATIC = HOR2 - V2SQ / TWO;
HTOTAL(I,J) = HSTATIC + C2SQ / TWO;
PTOTAL(I,J) = POR2 * (HTOTAL(I,J)/HOR2)^GAMAM;
PSTATIC = PTOTAL(I,J) * (HSTATIC/HTOTAL(I,J))GAMAM;
DENSOLD = DENSITY(I,J);
DENSITY(I,J) = PSTATIC/(RGAS*HSTATIC/CP);
CHANGE = abs(DENSOLD-DENSITY(I,J));
ERRDENS = max(CHANGE,ERRDENS);
ENTROPY(I,J) = CP * log ( HTOTAL(I,J) / HTOTAL1) -...
RGAS * log( PTOTAL(I,J) / PTOTAL1) + ENTROP1
end
end
%______________________________________________________________
% ASSEMBLE THE RHS FOR UNKNOWN NODES
% I.E. ALL NODES EXCEPT HUB/SHROUD
% ----------------------------------
ERRRHS = ZERO;
CHANGE = ZERO;
for I = 1,NSTATN
for J = 2,NSTRM-1
RHSOLD = RHS(I,J);
CU = RCU(I,J) / RADIUS(I,J);
TSTATIC = (HTOTAL(I,J)-(CZ(I,J)^2+CR(I,J)^2+CU^2)/TWO)/CP;
RHS(I,J) = - ONE/CZ(I,J) * TWOPI / XMASS *...
( CU/RADIUS(I,J) *(RCU(I,J+1) - RCU(I,J-1)) / TWOH...
- (HTOTAL(I,J+1) - HTOTAL(I,J-1))/ TWOH ...
+ TSTATIC*(ENTROPY(I,J+1)- ENTROPY(I,J-1)) / TWOH );
CHANGE = abs(RHS(I,J)-RHSOLD);
ERRRHS = max (ERRRHS, CHANGE);
end
end
