forked from GTSL-UC/T-Blade3
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathfuncNsubs.f90
1530 lines (1374 loc) · 58.1 KB
/
funcNsubs.f90
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
41
42
43
44
45
46
47
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
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
!List of functions and subroutines used in 3DBGB code.
!--------------------------Kiran Siddappaji
!---------------------------------------------------------
! Subroutine to display the welcome message
!---------------------------------------------------------
subroutine displayMessage
write(*,*)
write(*,*)'************************************************************'
write(*,*)'************************************************************'
write(*,*)'**** T-BLADE3:Turbomachinery BLADE 3D Geometry Builder ****'
write(*,*)'**** ****'
write(*,*)'**** Version 0.8 ****'
write(*,*)'**** ****'
write(*,*)'**** ...was also called as below till Aug 2016... ****'
write(*,*)'**** 3DBGB: 3 Dimensional Blade Geometry Builder ****'
write(*,*)'**** ****'
write(*,*)'**** Version 1.3 ****'
write(*,*)'**** ****'
write(*,*)'**** This software comes with ABSOLUTELY NO WARRANTY ****'
write(*,*)'**** ****'
write(*,*)'**** This is a program which generates a 3D blade... ****'
write(*,*)'**** ...shape and outputs 3D blade section files. ****'
write(*,*)'**** ****'
write(*,*)'**** Inputs: LE and TE curve(x,r), inlet angle, ****'
write(*,*)'**** exit angle, chord, tm/c, incidence, ****'
write(*,*)'**** deviation, secondary flow angles, ****'
write(*,*)'**** streamline coordinates:(x,r) ****'
write(*,*)'**** control points for sweep, lean, ****'
write(*,*)'**** blade scaling factor. ****'
write(*,*)'**** ****'
write(*,*)'**** Outputs: 3D blade sections (x,y,z), ****'
write(*,*)'**** 2D airfoils (mprime,theta). ****'
write(*,*)'**** ****'
write(*,*)'**** ---------------by Kiran Siddappaji ---- ****'
write(*,*)'**** ---------------by Mark G. Turner ---- ****'
write(*,*)'**** ------------------- [email protected] ---- ****'
write(*,*)'**** ---------------by Karthik Balasubramanian ---- ****'
write(*,*)'**** ---------------by Syed Moez Hussain Mahmood---- ****'
write(*,*)'**** ---------------by Ahmed Nemnem ---- ****'
write(*,*)'**** ---------------by Marshall C. Galbraith ---- ****'
write(*,*)'************************************************************'
write(*,*)'************************************************************'
write(*,*)
write(*,*) 'T-Blade3 Copyright (C) 2017 University of Cincinnati, developed by Dr. Mark Turner,'
write(*,*) 'Karthik Balasubramanian, Syed Moez Hussain, Ahmed Farid Nemnem, '
write(*,*) 'Kiran Siddappaji and Marshall C. Galbraith.'
write(*,*)
write(*,*) 'This program is free software; you can redistribute it and/or modify it under the '
write(*,*) 'terms of the GNU General Public License as published by the Free Software Foundation; '
write(*,*) 'either version 2 of the License, or (at your option) any later version.'
write(*,*)
write(*,*) 'This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;'
write(*,*) 'without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. '
write(*,*) 'See the GNU General Public License for more details. For the complete terms of the '
write(*,*) 'GNU General Public License, please see this URL: http://www.gnu.org/licenses/gpl-2.0.html'
write(*,*)
write(*,*)'************************************************************'
write(*,*)
return
end subroutine displayMessage
!---------------------------------------------------------
!---------------------------------------------------------
!function to calculate the inlet angle including incidence
!----------------------------------------------------------
real*8 function inBetaInci(inbeta,inci)
implicit none
real*8 inbeta(1),inci(1)
!Positive inlet angle
if(inbeta(1).ge.0.)then
inBetaInci = inbeta(1) - inci(1)
else
!Negative inlet angle
inBetaInci = inbeta(1) + inci(1)
endif
!print*,inci(1),inBetaInci(1)
end function
!----------------------------------------------------------
!----------------------------------------------------------
!function to calculate the exit angle including deviation
!----------------------------------------------------------
real*8 function outBetaDevn(inbeta,outbeta,devn)
implicit none
real*8 inbeta(1),outbeta(1),camber(1), devn(1)
!calculating camber
camber(1) = outbeta(1) - inbeta(1)
!negative camber
if(camber(1).le.0.)then
outBetaDevn = outbeta(1) - devn(1)
else
!positive camber
outBetaDevn = outbeta(1) + devn(1)
endif
!print*,devn(1),outbeta(1)
end function
!----------------------------------------------------------
!--------------------------------------------------------------
!subroutine to write the dimensional hub and casing streamlines
!--------------------------------------------------------------
subroutine hubTipStreamline(xhub,rhub,nphub,xtip,rtip,nptip,nsl,scf,casename)
implicit none
integer i,nphub,nptip,nsl
real*8 xhub(nphub,1),rhub(nphub,1),xtip(nptip,1),rtip(nptip,1)
real*8 scf
character*80 fname1,fname2
character*32 casename
print*,' Writing dimensional hub and casing streamlines...'
fname1 = 'hub.'//trim(casename)//'.sldcrv'
open(1,file = fname1,status ='unknown')
do i = 1,nphub
write(1,*)scf*xhub(i,1),0.0,scf*rhub(i,1)
enddo
close(1)
! Casing streamline
fname2 = 'casing.'//trim(casename)//'.sldcrv'
open(2,file = fname2,status ='unknown')
do i = 1,nptip
write(2,*)scf*xtip(i,1),0.0,scf*rtip(i,1)
enddo
close(2)
return
end subroutine
!--------------------------------------------------------------
!--------------------------------------------------------------
!subroutine to write the dimensional streamlines (without hub or casing streamlines)nemnem
!--------------------------------------------------------------
subroutine streamlines(xml,rml,np,scf,casename,ia)
implicit none
integer i,ia,np
real*8 xml(np),rml(np)
real*8 scf
character*80 fname
character*32 casename,temp
write(temp,*)ia
print*,' Writing dimensional streamline'//trim(adjustl(temp))//'...'
fname = 'streamlines'//trim(adjustl(temp))//'.'//trim(casename)//'.sldcrv'
open(1,file = fname,status ='unknown')
do i = 1,np
write(1,*)scf*xml(i),0.0,scf*rml(i)
enddo
close(1)
return
end subroutine
!--------------------------------------------------------------
!--------------------------------------------------------------
! subroutine calculating hub offset
!--------------------------------------------------------------
subroutine huboffset(mphub,x,r,dxds,drds,hub,nphub,scf,casename)
implicit none
integer i,nphub
real*8, intent(out):: mphub(nphub,1)
real*8 xhub(nphub,1),rhub(nphub,1),xms_hub(nphub,1),rms_hub(nphub,1)
real*8 x(nphub,1),r(nphub,1),dxds(nphub,1),drds(nphub,1),hub,scf,deltan
real*8 b,xnorm(nphub,1),rnorm(nphub,1),dxn(nphub,1),drn(nphub,1)
character*80 fname1, fname2
character*32 casename
!Calcultating the normal and offset coordinates
do i = 1, nphub
!print*,dxds(i,1)
xnorm(i,1) = drds(i,1)
rnorm(i,1) = -dxds(i,1)
! xnorm(1,1) = -drds(i,1)
! rnorm(1,1) = dxds(i,1)
!print*,'xnorm(i,1) rnorm(i,1)', xnorm(i,1), rnorm(i,1)
deltan = hub*1.
b = deltan/((xnorm(i,1)**2 + rnorm(i,1)**2)**0.5)
dxn(i,1) = b*(xnorm(i,1))
drn(i,1) = b*(rnorm(i,1))
!offset hub coordinates
xhub(i,1) = x(i,1) + dxn(i,1)
rhub(i,1) = r(i,1) + drn(i,1)
enddo
write(*,*)
!writing the offset coordinates to a file
fname1 = 'hub-offset.'//trim(casename)//'.dat'
open(1,file = fname1, status = 'unknown')
write(*,*)'Calculating hub-offset streamline coordinates...'
write(*,*)'Writing the coordinates to hub-offset.dat file'
!writing the dimensionless offset coordinates to a file
fname1 = 'hub-offset-dimlss.'//trim(casename)//'.dat'
open(2,file = fname1, status = 'unknown')
!Calculating the mprime coordinates of the offset streamline/construction line.
mphub(1,1)= 0.
write(2,*)xhub(1,1),rhub(1,1)
!mp(1,1) = mphub(1,1) ! replacing the old streamline m' to the offset ones
do i = 2, nphub
!print*,'abs(rhub(i,1)+rhub(i-1,1))',abs(rhub(i,1)+rhub(i-1,1))
! bug: when radius of streamline is ~= 0 ... Nemnem 6/3/2014
! if(rhub(i,1) < radius_tolerance.and.(xhub(i,1)-xhub(i-1,1)) /= 0.0) then
! print*,'Small streamline radius exception used ... (huboffset subroutine)'
! mphub(i,1) = mphub(i-1,1) + 2.*(xhub(i,1) - xhub(i-1,1))/(rhub(i,1)+rhub(i-1,1))
! else if ((xhub(i,1)-xhub(i-1,1)) == 0.0) then
! mphub(i,1) = mphub(i-1,1)
! else
! mphub(i,1) = mphub(i-1,1) + 2.*sqrt((rhub(i,1)-rhub(i-1,1))**2+(xhub(i,1)-xhub(i-1,1))**2)/(rhub(i,1)+rhub(i-1,1))
! endif
mphub(i,1) = mphub(i-1,1) + 2.*sqrt((rhub(i,1)-rhub(i-1,1))**2+(xhub(i,1)-xhub(i-1,1))**2)/(rhub(i,1)+rhub(i-1,1))
write(1,*)scf*xhub(i,1),0.0,scf*rhub(i,1)
write(2,*)xhub(i,1),rhub(i,1)
!mp(i,1) = mphub(i,1) ! replacing the old streamline m' to the offset ones
enddo
! Splining the offset xhub, rhub:
call spline(xhub(1,1),xms_hub(1,1),mphub(1,1),nphub, 999.0, -999.0)
call spline(rhub(1,1),rms_hub(1,1),mphub(1,1),nphub, 999.0, -999.0)
! over writing the xm, rm values with hub spline coefficients:
x = xhub
r = rhub
dxds = xms_hub
drds = rms_hub
! print*,'xms_hub',x
! print*,'rms_hub',r
close(2)
close(1)
return
end subroutine
!--------------------------------------------------------------
!--------------------------------------------------------------
! subroutine calculating tip offset
!--------------------------------------------------------------
subroutine tipoffset(mptip,x,r,dxds,drds,tip,nptip,scf,nsl,casename)
implicit none
integer i,nptip,nsl
real*8,intent(out) :: mptip(nptip,1)
real*8 xtip(nptip,1),rtip(nptip,1),xms_tip(nptip,1),rms_tip(nptip,1)
real*8 x(nptip,1),r(nptip,1),dxds(nptip,1),drds(nptip,1),tip,scf
real*8 b,xnorm(nptip,1),rnorm(nptip,1),dxn(nptip,1),drn(nptip,1),deltan
character*80 fname1, fname2
character*32 casename
!Calcultating the normal and offset coordinates
do i = 1, nptip
!print*,dxds(i,1)
xnorm(i,1) = drds(i,1)
rnorm(i,1) = -dxds(i,1)
!print*, xnorm(i,1),rnorm(i,1)
deltan = tip*1.
b = deltan/((xnorm(i,1)**2 + rnorm(i,1)**2)**0.5)
dxn(i,1) = b*(xnorm(i,1))
drn(i,1) = b*(rnorm(i,1))
!offset tip coordinates
xtip(i,1) = x(i,1) - dxn(i,1)
rtip(i,1) = r(i,1) - drn(i,1)
enddo
write(*,*)
!writing the offset coordinates to a file
fname1 = 'tip-offset.'//trim(casename)//'.dat'
open(1,file = fname1, status = 'unknown')
write(*,*)'Calculating tip-offset streamline coordinates...'
write(*,*)'Writing the coordinates to tip-offset.dat file'
!writing the dimenstionless offset coordinates to a file
fname1 = 'tip-offset-dimlss.'//trim(casename)//'.dat'
open(2,file = fname1, status = 'unknown')
write(2,*)xtip(1,1),rtip(1,1)
!Calculating the mprime coordinates of the offset streamline/construction line.
mptip(1,1)= 0.
!mp(1,nsl) = mptip(1,nsl) ! replacing the old streamline m' to the offset ones
do i = 2, nptip
! if(rtip(i,1) < radius_tolerance.and.((xtip(i,1)-xtip(i-1,1)) /= 0)) then
! mptip(i,1) = mptip(i-1,1) + 2.*(xtip(i,1) - xtip(i-1,1))/(rtip(i,1)+rtip(i-1,1))
! else if ((xtip(i,1)-xtip(i-1,1)) == 0.0) then
! mptip(i,1) = mptip(i-1,1)
! else
! mptip(i,1) = mptip(i-1,1) + 2.*sqrt((rtip(i,1)-rtip(i-1,1))**2+(xtip(i,1)-xtip(i-1,1))**2)/(rtip(i,1)+rtip(i-1,1))
! endif
mptip(i,1) = mptip(i-1,1) + 2.*sqrt((rtip(i,1)-rtip(i-1,1))**2+(xtip(i,1)-xtip(i-1,1))**2)/(rtip(i,1)+rtip(i-1,1))
write(1,*)scf*xtip(i,1),0.0,scf*rtip(i,1)
write(2,*)xtip(i,1),rtip(i,1)
!mp(i,nsl) = mptip(i,nsl) ! replacing the old streamline m' to the offset ones
!print*,mptip(i,1)
enddo
! Splining the offset xtip, rtip:
call spline(xtip(1,1),xms_tip(1,1),mptip(1,1),nptip, 999.0, -999.0)
call spline(rtip(1,1),rms_tip(1,1),mptip(1,1),nptip, 999.0, -999.0)
! overwriting xm, rm with new tip spline coefficients:
x = xtip
r = rtip
dxds = xms_tip
drds = rms_tip
! print*,'xms_tip',x
! print*,'rms_tip',r
close(2)
close(1)
return
end subroutine
!--------------------------------------------------------------
!--------------------------------------------------------------
!-----INTERSECTION_POINTS--------------------------------------
!--------------------------------------------------------------
SUBROUTINE INTERP(XA,YA,XB,YB,XC,YC,XD,YD,XINT,YINT)
real*8 XA,YA,XB,YB,XC,YC,XD,YD,XINT,YINT
XINT=(XA*(YB-YA)+YC*(XB-XA)-((YD-YC)/(XD-XC))*XC*(XB-XA)-YA*(XB-XA))/((YB-YA)-(((YD-YC)*(XB-XA))/(XD-XC)))
YINT = YC+((XINT-XC)/(XD-XC))*(YD-YC)
RETURN
END subroutine
!*******************************************************
!*******************************************************
!*******************************************************
subroutine outputfiledata(bladedata,nsl,amount_data,throat_pos,casename,units)
integer nsl,amount_data,js
real*8 bladedata(amount_data,nsl)
real*8 in_beta(nsl),out_beta(nsl)
character*20 throat_pos(nsl)
character*32 casename
character*80 file1
character(len=2) :: units
file1 = 'blade_section_data.'//trim(casename)//'.dat'
open(unit= 100,file=file1, form="formatted")
write(100,*)trim(casename)
write(100,*),'Blade sections Data:'
write(100,*),'---------------------'
write(100,*)
do js = 1,nsl
if (js == 1) then
if (units == 'mm') then
if (throat_pos(nsl) == 'le') then
write(100,202),'section','span','betaZ*le(deg)','betaZ*te(deg)','betaM*le(deg)','betaM*te(deg)','chord [mm]','sweep','lean','area[mm^2]','lethk [mm]','tethk [mm]','throat [mm]','(r*delta_theta)LE','Geom Zweifel','le/te/btween/none'
elseif (throat_pos(nsl) == 'te') then
write(100,202),'section','span','betaZ*le(deg)','betaZ*te(deg)','betaM*le(deg)','betaM*te(deg)','chord [mm]','sweep','lean','area[mm^2]','lethk [mm]','tethk [mm]','throat [mm]','(r*delta_theta)TE','Geom Zweifel','le/te/btween/none'
else
write(100,202),'section','span','betaZ*le(deg)','betaZ*te(deg)','betaM*le(deg)','betaM*te(deg)','chord [mm]','sweep','lean','area[mm^2]','lethk [mm]','tethk [mm]','throat [mm]','(r*delta_theta)WT','Geom Zweifel','le/te/btween/none'
endif
elseif (units == 'cm') then
if (throat_pos(nsl) == 'le') then
write(100,202),'section','span','betaZ*le(deg)','betaZ*te(deg)','betaM*le(deg)','betaM*te(deg)','chord [cm]','sweep','lean','area[cm^2]','lethk [cm]','tethk [cm]','throat [cm]','(r*delta_theta)LE','Geom Zweifel','le/te/btween/none'
elseif (throat_pos(nsl) == 'te') then
write(100,202),'section','span','betaZ*le(deg)','betaZ*te(deg)','betaM*le(deg)','betaM*te(deg)','chord [cm]','sweep','lean','area[cm^2]','lethk [cm]','tethk [cm]','throat [cm]','(r*delta_theta)TE','Geom Zweifel','le/te/btween/none'
else
write(100,202),'section','span','betaZ*le(deg)','betaZ*te(deg)','betaM*le(deg)','betaM*te(deg)','chord [cm]','sweep','lean','area[cm^2]','lethk [cm]','tethk [cm]','throat [cm]','(r*delta_theta)WT','Geom Zweifel','le/te/btween/none'
endif
elseif ((units == 'm ').or.(units == 'm)')) then
if (throat_pos(nsl) == 'le') then
write(100,202),'section','span','betaZ*le(deg)','betaZ*te(deg)','betaM*le(deg)','betaM*te(deg)','chord [ m]','sweep','lean','area[ m^2]','lethk [ m]','tethk [ m]','throat [ m]','(r*delta_theta)LE','Geom Zweifel','le/te/btween/none'
elseif (throat_pos(nsl) == 'te') then
write(100,202),'section','span','betaZ*le(deg)','betaZ*te(deg)','betaM*le(deg)','betaM*te(deg)','chord [ m]','sweep','lean','area[ m^2]','lethk [ m]','tethk [ m]','throat [ m]','(r*delta_theta)TE','Geom Zweifel','le/te/btween/none'
else
write(100,202),'section','span','betaZ*le(deg)','betaZ*te(deg)','betaM*le(deg)','betaM*te(deg)','chord [ m]','sweep','lean','area[ m^2]','lethk [ m]','tethk [ m]','throat [ m]','(r*delta_theta)WT','Geom Zweifel','le/te/btween/none'
endif
else
print*,'Error in Blade scaling factor ...'
print*,'Check the Scaling factor in the input file ...'
stop
endif
endif
write(100,201),js,bladedata(1:(amount_data-1),js),throat_pos(js)
enddo
write(100,*),'Total Volume of the blade =',bladedata(amount_data,nsl),'[m^3]'
close(100)
202 format((A7,2x),(A4,11x),4(A14,3x),(A13,1x),2(A6,5x),(A13,2x),(A10,2x),(A10,1x),(A11),2x,A17,2x,A12,2x,A17)
201 format(1x,I3,4x,(f11.8,4x),4(sf14.8,3x),3x,8(f11.8,1x),7x,f11.8,4x,A5)
return
end subroutine
!***************************************************************
!*******************************************************
!***************************************************************
subroutine throatindex(throat_pos,throat_index,n_normal_distance,js,nsl)
implicit none
integer, intent(inout):: js
integer, intent(in):: nsl,throat_index(nsl),n_normal_distance
character*20, intent(out):: throat_pos(nsl)
if(throat_index(js) == 0) then
throat_pos(js) = 'none'
print*,'No Throat found'
!exit
elseif(throat_index(js) < 0.25*n_normal_distance) then
throat_pos(js) = 'le'
print*,'throat_index',js,throat_index(js)
elseif(throat_index(js) > 0.75*n_normal_distance) then
throat_pos(js) = 'te'
print*,'throat_index',js,throat_index(js)
else
throat_pos(js) = 'btween'
print*,'throat_index',js,throat_index(js)
endif! end if for throat options
return
end subroutine
!***************************************************************
!*******************************************************
!***************************************************************
subroutine cascade_nondim_file(msle,mste,mprime_ble,mprime_bte,chordm,pitch,nsl,ibrow,casename)
implicit none
integer, intent(in) :: nsl,ibrow
integer ia
character*80 file1
character*32 casename
real*8, intent(in) :: msle(nsl),mste(nsl),pitch
real*8, intent(in) :: mprime_ble(nsl),mprime_bte(nsl),chordm(nsl)
file1 = 'cascade_nondim.'//trim(casename)//'.dat'
open(13,file = file1,status ='unknown')
write(13,*)trim(casename)
write(13,*)'Blade row:',ibrow
write(13,*)"Non-dimensional quantities:"
write(13,*)"section m'sLE m'sTE m'blade_LE m'blade_TE chord pitch "
do ia = 1, nsl
write(13,101)ia,msle(ia),mste(ia),mprime_ble(ia),mprime_bte(ia),chordm(ia),pitch
enddo
close(13)
101 format(I2,2X,6(f19.16,1x))
return
end subroutine
!***************************************************************
!**********************************************************************************************
! BLADEGEN SUBROUTINES ------------
!**********************************************************************************************
!
subroutine aecalc(n,x,y,t, itype, area,xcen,ycen,ei11,ei22,apx1,apx2)
dimension x(n),y(n),t(n)
!
!---------------------------------------------------------------
! calculates geometric properties of shape x,y
!
! input:
! n number of points
! x(.) shape coordinate point arrays
! y(.)
! t(.) skin-thickness array, used only if itype = 2
! itype = 1 ... integration is over whole area dx dy
! = 2 ... integration is over skin area t ds
!
! output:
! xcen,ycen centroid location
! ei11,ei22 principal moments of inertia
! apx1,apx2 principal-axis angles
!---------------------------------------------------------------
!
sint = 0.0
aint = 0.0
xint = 0.0
yint = 0.0
xxint = 0.0
xyint = 0.0
yyint = 0.0
!
do 10 io = 1, n
if(io.eq.n) then
ip = 1
else
ip = io + 1
endif
!
dx = x(io) - x(ip)
dy = y(io) - y(ip)
xa = (x(io) + x(ip))*0.50
ya = (y(io) + y(ip))*0.50
ta = (t(io) + t(ip))*0.50
!
ds = sqrt(dx*dx + dy*dy)
sint = sint + ds
if(itype.eq.1) then
!-------- integrate over airfoil cross-section
da = ya*dx
aint = aint + da
xint = xint + xa *da
yint = yint + ya *da/2.0
xxint = xxint + xa*xa*da
xyint = xyint + xa*ya*da/2.0
yyint = yyint + ya*ya*da/3.0
else
!-------- integrate over skin thickness
da = ta*ds
aint = aint + da
xint = xint + xa *da
yint = yint + ya *da
xxint = xxint + xa*xa*da
xyint = xyint + xa*ya*da
yyint = yyint + ya*ya*da
endif
!
10 continue
!
area = aint
!
if(aint .eq. 0.0) then
xcen = 0.0
ycen = 0.0
ei11 = 0.0
ei22 = 0.0
apx1 = 0.0
apx2 = atan2(1.0,0.0)
return
endif
!
!
!---- calculate centroid location
xcen = xint/aint
ycen = yint/aint
!
!---- calculate inertias
eixx = yyint - ycen*ycen*aint
eixy = xyint - xcen*ycen*aint
eiyy = xxint - xcen*xcen*aint
!
!---- set principal-axis inertias, ei11 is closest to "up-down" bending inertia
eisq = 0.25*(eixx - eiyy )**2 + eixy**2
sgn = sign( 1.0 , eiyy-eixx )
ei11 = 0.5*(eixx + eiyy) - sgn*sqrt(eisq)
ei22 = 0.5*(eixx + eiyy) + sgn*sqrt(eisq)
!
if(eisq/(ei11*ei22) .lt. (0.001*sint)**4) then
!----- rotationally-invariant section (circle, square, etc.)
apx1 = 0.0
apx2 = atan2(1.0,0.0)
!
else
c1 = eixy
s1 = eixx - ei11
c2 = eixy
s2 = eixx - ei22
!
sgn1 = sign( 1.0 , c1 )
sgn2 = sign( 1.0 , c2 )
!
apx1 = atan2(sgn1*s1,sgn1*c1)
apx2 = atan2(sgn2*s2,sgn2*c2)
!
endif
return
end subroutine
!*************************************************************
!*******************************************************************
subroutine stacking(xb,yb,xbot,ybot,xtop,ytop,js,np,stack_switch,stack,stk_u,stk_v,area,LE)
!---------------------------------------------------------------------------------------------------
! Airfoil Stacking options-----------------------Nov 2012
! Sample Input value : -053075 means stacked at 53% chord and 75% below meanline.
! Sample Input value : +100035 means stacked at 100% chord (TE) and 35% above meanline.
! Sample Input value : +200100 means stacked at area centroid and 100% above meanline.
! stack = above stacking number
! stack_switch = switch to activate variable radial stacking.
! xb = u, yb = v
! xb = xb - u_stack
! yb = yb - v_stack
! u_stack is the offset along the chord.
! v_stack is the offset above or below meanline.
! Stacking on meanline: v_stack = 0
! stku = stack position as %age of the chord, Range: [0,100].
! stkv = stack position from the meanline (stk_v=0). +ve is above the meanline and -ve is belove the meanline. Range:[-100,100].
! Stacking at the area centroid: value is 200.
! print*,'stack:',stack
implicit none
integer i,j,k,nx,np_side,np,nphalf,ny,LE
integer, intent(in)::js,stack,stack_switch
parameter (nx=500,ny=300)
!real, intent (inout):: xb(nx),yb(nx),xbot(ny),ybot(ny),xtop(ny),ytop(ny)
real, intent (inout):: xb(np),yb(np),xbot((np+1)/2),ybot((np+1)/2),xtop((np+1)/2),ytop((np+1)/2)
real, intent(in):: stk_u(1), stk_v(1)
real umin,umax,sb(nx)
real u_stack,v_stack,vtop_stack,vbot_stack,v_zero_stack
real const_stk_u,const_stk_v,stku,stkv
real area,ucen,vcen,ei11,ei22,apx1,apx2
! Calculating the area centroid cordinates of the airfoil ucen,vcen.
print*,np
call aecalc(np,xb,yb,sb,1,area,ucen,vcen,ei11,ei22,apx1,apx2)
!
if(mod(np,2).eq.0)then
np_side = np/2 ! For even no. of points LE = np/2. So, the current formula (np+1)/2 ==> np/2 with this change.
else
! points on bottom and top curve
np_side = (np+1)*0.5
endif
if(stack_switch.eq.0)then !constant radial stacking.
const_stk_u = floor(abs(real(stack))/1000)
stku = const_stk_u
if(stack.ge.0)then
const_stk_v = real(stack)-(const_stk_u*1000)
stkv = const_stk_v
else
const_stk_v = real(stack)+(const_stk_u*1000)
stkv = const_stk_v
endif
elseif(stack_switch.eq.1)then ! variable radial stacking.Kiran 7/52/13
stku = stk_u(js)
stkv = stk_v(js)
endif
write(*,*)
print*,'Stacking position on the chord as % of chord:',stku
print*,'Stacking position in % above(+ve) or below(-ve) meanline :',stkv
! stacking on the chord...........
if(stku.eq.200.)then ! 200 is for stacking at the area centroid.
u_stack = ucen
v_stack = vcen
print*,'Stacking at the area centroid of the airfoil...'
goto 16
else
! Calculating u_stack using the % chord stacking information
!nphalf = (np+1)/2
if(LE.eq.2)then ! sting LE option
print*,'LE with sting.'
! nphalf = nphalf+2
!xb(nphalf) = 0.5*(xb(nphalf-2) + xb(nphalf+2))
!yb(nphalf) = 0.5*(yb(nphalf-2) + yb(nphalf+2))
endif
! u_stack = xb(nphalf) + (real(stku)/100)*abs(sqrt((xb(nphalf) - xb(np))**2 + (yb(nphalf) - yb(np))**2))
! u_stack = xb(nphalf) + (real(stku)/100)*abs(xb(nphalf) - xb(np))! + (yb(nphalf) - yb(np))**2))
u_stack = xb(1)*(real(stku)/100)
! write(*,*)
! print*,'ustack',u_stack
endif
!write(*,*)
!
! stacking on/above/below meanline.......
! Calculating vtop_stack and vbot_stack for v_stack.
j = 1
do i = 1, np_side-1 ! from 1 to 99 (top curve) for 199 points [moving from right to left]
umin = xbot(i)
umax = xbot(i+1)
if(u_stack.ge.umin.and.u_stack.le.umax) then
j = i
exit
endif
enddo
print*,'top curve points index range for % chord stacking:',j,j+1
vbot_stack = ybot(j) + ((u_stack - xbot(j))*(ybot(j+1) - ybot(j))/(xbot(j+1) - xbot(j))) !*sqrt((xbot(j+1) - xbot(j))**2 + (ybot(j+1) - ybot(j))**2)
k = 1
do i = 1, np_side-1 ! from 101 to 198 for 199 points (bottom curve)
umin = xtop(i)
umax = xtop(i+1)
if(u_stack.ge.umin.and.u_stack.le.umax) then
k = i
exit
endif
enddo
print*,'bottom curve points index range for % chord stacking:',k,k+1
vtop_stack = ytop(k) + ((u_stack - xtop(k))*(ytop(k+1) - ytop(k))/(xtop(k+1) - xtop(k))) !*sqrt((xtop(k+1) - xtop(k))**2 + (ytop(k+1) - ytop(k))**2)
! Stacking on meanline: v_stack = 0
! v_zero_stack = average of yb(istack) and yb(1+np - istack)
! Above meanline stack: vstack= v_zero_stack - %age(stack)*(distance between meanline and 100% ABOVE meanline @ istack)
! Below meanline stack: vstack= v_zero_stack - %age(stack)*(distance between meanline and 100% BELOW meanline @ istack)
v_zero_stack = (vtop_stack + vbot_stack)/2
print*, "vtop_stack vbot_stack", vtop_stack, vbot_stack
if(stku.ne.200.and.stkv.gt.0)then! above meanline stacking
v_stack = v_zero_stack + (real(stkv)/100)*(abs(vtop_stack - v_zero_stack))
print*,'+ve stack v_stack',v_stack,(real(stkv)/100)
elseif(stku.ne.200.and.stkv.eq.0)then
v_stack = 0.!v_zero_stack
elseif(stku.ne.200.and.stkv.lt.0)then !below meanline stacking
v_stack = v_zero_stack + (real(stkv)/100)*(abs(vbot_stack - v_zero_stack))
print*,'v_stack',v_stack,(real(stkv)/100)
endif
write(*,*)
print*,'u_stack v_stack', u_stack, v_stack
write(*,*)
!-----stacked coordinates
16 do i = 1,np
xb(i) = xb(i) - u_stack
yb(i) = yb(i) - v_stack
!print*,xb(i),yb(i)
enddo
return
end subroutine
!***************************************************
!*******************************************************
!*******************************************************
subroutine rotate(xb,yb,x,y,angle)
implicit none
real x,y
real, intent(out) :: xb(1),yb(1)
real angle
xb(1) = x*cos(-angle) + y*sin(-angle)
yb(1) = y*cos(-angle) - x*sin(-angle)
!print*,xb(1),yb(1)
return
end
!***************************************************
!*******************************************************
!*******************************************************
subroutine rotate2(xb,yb,x,y,angle)
implicit none
real x,y
real, intent(out) :: xb,yb
real angle
xb = x*cos(-angle) + y*sin(-angle)
yb = y*cos(-angle) - x*sin(-angle)
!print*,xb(1),yb(1)
return
end
!
!**************************************************
!*******************************************************
!*******************************************************
real function scaled(x,scalefactor)
implicit none
real, intent(in):: x,scalefactor
scaled = x*scalefactor
return
end
!
!*******************************************************
!*******************************************************
!*******************************************************
subroutine bladesection(xb,yb,np,nbls,TE_del,sinls,sexts,chrdd,fext,js,pitch,mble,mbte,airfoil)
implicit none
integer TE_del,nx,nbls,ii,i,js,np,np_side
parameter(nx=500)
character*32 fname,fext,blext(100)
character*20 airfoil
real*8 sexts,sinls,chrdd,pitch,pi,x1,y1,x2,y2
real*8, intent(inout) :: xb(nx),yb(nx)
real*8 mble,mbte
if(mod(np,2).eq.0)then
np_side = np/2 ! For even no. of points LE = np/2. So, the current formula (np+1)/2 ==> np/2 with this change.
else
! points on bottom and top curve
np_side = (np+1)*0.5
endif
!If the TE points do not coincide.
x1 = xb(1)
y1 = yb(1)
x2 = xb(np)
y2 = yb(np)
if(len_trim(airfoil).eq.4)then ! for NACA cases to have a closed straight TE.
xb(1) = 0.5*(x1 + x2)
yb(1) = 0.5*(y1 + y2)
xb(np) = xb(1)
yb(np) = yb(1)
endif
if((x1.ne.x2).and.(y1.ne.y2))then
xb(1) = 0.5*(xb(np) + xb(1))
yb(1) = 0.5*(yb(np) + yb(1))
xb(np) = xb(1)
yb(np) = yb(1)
endif
pi = 4.*atan(1.0)
! Calculation of the pitch between adjacent blades:
pitch = 2*pi/nbls
! Switch for deleting the trailing edge coordinates for each blade:
!TE_del = 1 --- I moved it in the parameters nemnem
if(TE_del == 1) then
ii = 18 ! ii is the number of points deleted from the Trailing edge.
!print*,'The switch to delete Trailing edge from blade files is ON'
elseif (TE_del == 0) then
ii = 0
!print*,'The switch to delete Trailing edge from blade files is OFF'
endif
!Assiging the LE and TE to variables
mble = xb(np_side)
mbte = xb(np)
!---- output
fname = 'blade.'//trim(fext)
open(2,file=fname,status='unknown')
write(2,100) trim(fext)
write(2,101) sinls,sexts,0.5*chrdd,0.5*chrdd,pitch
! print*,'np before blade file',np
do i = ii+1,np-ii
write(2,*) xb(i),yb(i)
enddo
close(2)
100 format(a)
101 format(5(f19.16,1x))
102 format(2(f35.16,1x))
! 105 format(' Blade Area : ',f8.5)
!106 format(' Blade Centroid : ',f8.5,1x,f8.5)
return
end subroutine
!*****************************************************
!*****************************************************
!*****************************************************
SUBROUTINE st_line_intersection(XA,YA,XB,YB,XC,YC,XD,YD,XINT,YINT)
integer i,j
real*8 XA,YA,XB,YB,XC,YC,XD,YD,xmin,ymin,xmax,ymax
real*8 , intent(out):: XINT,YINT
! xmin = min(XA,XB,XC,XD)
! ymin = min(YA,YB,YC,YD)
! xmax = max(XA,XB,XC,XD)
! ymax = max(YA,YB,YC,YD)
!
xint =(XA*(YB-YA)+YC*(XB-XA)-((YD-YC)/(XD-XC))*XC*(XB-XA)-YA*(XB-XA))/((YB-YA)-(((YD-YC)*(XB-XA))/(XD-XC)))
yint = YC+((XINT-XC)/(XD-XC))*(YD-YC)
RETURN
END subroutine
!*******************************************************
!********************************************************
!********************************************************
! This subroutine is to calculate the 2D throat between nondimensional blade sections.
! The approach here is to find the perpendicular distance on the pitch line(camber at half pitch)
! between 2 adj blades.
! calculating the throat from bottom point on the blade ---(for +ve ainl and -ve aext)---
subroutine throat_calc_pitch_line(xb,yb,np,camber,angle,sang,u,pi,pitch,throat_coord, mouth_coord,exit_coord, &
min_throat_2D,throat_index,n_normal_distance,casename,js,nsl,develop,isdev)
! np = 2*np_sidee-1
! All data used is nondimensional
implicit none
integer i,j,k,np,np_sidee,throat_index,n_normal_distance,js,nsl
integer i_interup1,i_interup2,i_interdwn1,i_interdwn2
real v1_top((np+1)/2),u1_top((np+1)/2),v2_bot((np+1)/2),u2_bot((np+1)/2)
real xb(np),yb(np),yb_upper(np)
real xint_up,yint_up,xint_dwn,yint_dwn
real pitch, throat_coord(4),mouth_coord(4),exit_coord(4)
real inter_coord(4,((np+1)/2))
real pitch_line((np+1)/2), camber((np+1)/2),angle((np+1)/2),u((np+1)/2)
real camber_upper((np+1)/2)
real min_throat_2D,pi,sang,angle_up
real x_interup,y_interup,x_interdwn,y_interdwn
real, allocatable, dimension(:) :: throat
character*80 file4
character*32 casename,develop
logical isdev
if(mod(np,2).eq.0)then
np_sidee = np/2 ! For even no. of points LE = np/2. So, the current formula (np+1)/2 ==> np/2 with this change.
else
! points on bottom and top curve
np_sidee = (np+1)/2
endif
!np_sidee = (np+1)/2
print*,'np_sidee',np_sidee
! intializing the values to 0...
n_normal_distance = 0
throat_index = 0
inter_coord = 0.
throat_coord = 0.; mouth_coord = 0. ;exit_coord = 0.
! Calculating the pitch line in mid way between 2 blades:
pitch_line = camber + (0.5*pitch)
! Create the upper blade:
yb_upper = yb + pitch
camber_upper = camber + pitch
! Finding the upper point using the point and angle:
! after rotation by stagger, 'sang' should be taking into consideration,
! write(85,*),'u1_top v1_top u2_bot v2_bot'
do i = 1, np_sidee
v1_top(i) = pitch_line(i)+ 0.5*pitch
u1_top(i) = 1/tan(angle(i)+sang+pi/2.)*(v1_top(i)-pitch_line(i))+u(i)
v2_bot(i) = pitch_line(i)- 0.5*pitch
u2_bot(i) = 1/tan(angle(i)+sang-pi/2.)*(v2_bot(i)-pitch_line(i))+u(i)
!write(85,*),u1_top(i),v1_top(i)
!write(85,*),u2_bot(i),v2_bot(i)
enddo
! Identify the throat by intersection point with upper surface of the blade:
! intials:
k = 1
!write(85,*)'intersection points'
do j = 1, np_sidee
i_interup1 = 0
i_interdwn1 = 0
do i = np_sidee,np-1 ! intersection with upper airfoil
call st_line_intersection(u(j),pitch_line(j),u1_top(j),v1_top(j),xb(i),yb_upper(i),xb(i+1),yb_upper(i+1),xint_up,yint_up)
!print*,'xint_up,yint_up',xint_up,yint_up
if((xint_up >= xb(i) ).and.(xint_up <= xb(i+1)))then
!print*,"correct point of intersection ..."
!print*,'xint_up,yint_up ----',xint_up,yint_up
x_interup = xint_up ;y_interup = yint_up
!write(85,*),'xint_up,yint_up',xint_up,yint_up
i_interup1 = i ;i_interup2 = i+1
!print*,'i_interup',i_interup1,i_interup2
endif
enddo
do i = 1,np_sidee-1 ! intersection with lower airfoil
call st_line_intersection(u(j),pitch_line(j),u2_bot(j),v2_bot(j),xb(i),yb(i),xb(i+1),yb(i+1),xint_dwn,yint_dwn)
!print*,'xint_dwn,yint_dwn',xint_dwn,yint_dwn
if((xint_dwn <= xb(i) ).and.(xint_dwn >= xb(i+1)))then
!print*,"correct point of intersection ..."
!print*,'xint_dwn,yint_dwn ----',xint_dwn,yint_dwn
x_interdwn = xint_dwn ;y_interdwn = yint_dwn
!write(85,*),'xint_dwn,yint_dwn',xint_dwn,yint_dwn
i_interdwn1 = i ;i_interdwn2 = i+1
!print*,'i_interdwn',i_interdwn1,i_interdwn2
endif
enddo
if((i_interup1.ne.0).and.(i_interdwn1.ne.0)) then
inter_coord(1,k) = x_interup;inter_coord(2,k) = y_interup ! first intersection point
inter_coord(3,k) = x_interdwn ;inter_coord(4,k) = y_interdwn ! second intersection point
!print*,'inter_coord(:,k)',inter_coord(:,k),k
! write(85,*),inter_coord(1:2,k)
! write(85,*),inter_coord(3:4,k),k
n_normal_distance = k ! number of throats for each section
k = k+1
endif
enddo
print*, 'n_normal_distance =',n_normal_distance
if(n_normal_distance == 0) then
print*, 'No throats found because of Low number of blades.'
return
endif
! Writing the Mouth and Exit areas:
mouth_coord = inter_coord(:,1)
exit_coord = inter_coord(:,n_normal_distance)
print*,'number of intersection points (k) =',k-1,'from np_side of',np_sidee
if (allocated(throat)) deallocate(throat)
Allocate (throat(n_normal_distance))
! Calculation of the throat:
throat(1) = sqrt((inter_coord(1,1)-inter_coord(3,1))**2+(inter_coord(2,1)-inter_coord(4,1))**2) ! Nondimensional
min_throat_2D = throat(1)
throat_index = 1
do k = 2,n_normal_distance
throat(k) = sqrt((inter_coord(1,k)-inter_coord(3,k))**2+(inter_coord(2,k)-inter_coord(4,k))**2) ! Nondimensional
if(throat(k) < min_throat_2D) then
min_throat_2D = throat(k)
throat_index = k
throat_coord = inter_coord(:,k)
endif
enddo
! Writing to a file for debugging
!-----------------------------------------------
if(isdev)then
print*,""
print*,'Writing non-dimensional throat points to a file for debugging...'
print*,""
if(js==1) then
file4 = 'throat_points.'//trim(casename)//'.txt'
open(unit= 85,file=file4, form="formatted")
print*,file4
write(85,*)'pitch',pitch
endif
write(85,*)'section',js
write(85,*)'u camber upper_camber'
do i = 1, np_sidee