C********************************************************************** C* Estimation of Roll Damping for conventional cargo ship & barge * C* * C* SA_RDPM.BAS * C* since 2007.09, developed & coded by Toru KATAYAMA, Dr.Eng * C* 2021.11, correction of eddy making componen * C* 2024.09, correction of YA in bilge^keel roll damping * C* Osaka Metroplitan University * C********************************************************************** C PROGRAM MAIN IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION NUE,L,KAI,L0,L1,L2,NABLA,KW,KAI1,LTHETA,MGA,KN & ,LO,LR,KG,LB1,LB2,LB3,LB2S0,M1,M2,M3,M4,M5,M6,M7,M8,LOMEGA INTEGER BWOSM,TYPE_BK,RDP,CSDMIN,SSBK1,SSBK2,P C SC:number of SS, PT:number of data in a SS, NNRA:number of wave length for culculation C NNFN:number of Fn for culculation, NNKAI:number of wave direction for culculation C NNRD:number of roll angles for culculation PARAMETER SC=200,PT=200,NNRA=70,NNFN=20,NNKAI=35,NNRD=30 C DIMENSION CSD(SC),NUMPT(SC),HT(0:SC,0:PT),HB(0:SC,0:PT) DIMENSION THETA1(NNRD),RAM(NNRA),FUN(NNFN),KAI(NNKAI) DIMENSION XX(SC,PT*2),YY(SC,PT*2) DIMENSION X(SC),HO(SC),SIG(SC),BX(SC),DX(SC) DIMENSION PERI(NNRA,NNFN,NNKAI),BE2(SC) DIMENSION RMAX1(2),V(2),CR(SC),CR1(SC) DIMENSION X1(30),YI(30),XBK(0:SC),RATIO(SC) DIMENSION X11(SC),Y11(SC),WX(SC),WY(SC) DIMENSION BKT(SC),BKB(SC),NUBK(SC),BBKHAT(SC) DIMENSION BBYD(SC),SGM(SC),BB0(SC),TT0(SC),WLHB(SC),WLAR(SC) DIMENSION GLEN(SC) C C C C****** C if BWOSM=0 then wave making component is not culculated BWOSM=0 C kinematic viscous coefficient NUE=1.14E-06 C circular constant PAI=3.14159 C Acceleration of gravity G=9.8 RO=102.0 C C OPEN(1,FILE='OSM2.RDP') C****** data of hull form C Lpp, Breadth, Draft READ(1,*)L,B,D C number of square stations, the number of mid-ship READ(1,*)NUMCS,IMID DO I=1 , NUMCS C SS READ(1,*)CSD(I) C number of datum in a square station (ss) READ(1,*)NUMPT(I) DO J=1 , NUMPT(I) C hight from Base Line, breadth from Center Line(please input from BL) READ(1,*)HT(I,J),HB(I,J) END DO END DO C****** data for bilgekeel C breadth of bilgekeel, afr-end(ss), forward-end(ss) (in Japan, AP=ss0, FP=ss10) READ(1,*)BBK,XBK1,XBK2 C Selection of the calculation method (1: simplified, 2:integration) READ(1,*)TYPE_BK IF( TYPE_BK.NE.1)then C the aft-end and the forward-end of the ss number of the section where the bilge keel exists on hull's ss. READ(1,*)SSBK1,SSBK2 DO I=SSBK1 , SSBK2 C The attached transverse position of the bilge-keel from BL READ(1,*)BKB(I) END DO END IF C****** eddy making component for two type ships (1: conventuonal ship, 2:barge) READ(1,*)RDP C****** OGD=(Draft-KG)/Draft READ(1,*)OGD C****** number of roll amplitude for culculation, roll amplitude(degree) READ(1,*)NTHETA DO I=1 , NTHETA READ(1,*)THETA1(I) END DO C****** Fn for culculation READ(1,*)NFN DO I=1 , NFN READ(1,*)FUN(I) END DO C****** wavelength/Lpp for culculation READ(1,*)NRA DO I=1 , NRA READ(1,*)RAM(I) END DO C****** wave direction for culculation (degree) READ(1,*)NKAI DO I=1 , NKAI READ(1,*)KAI(I) END DO CLOSE(1) C===== end of loading the data file C C=====in order to culculate bilge keel component by integration method, following culculation is required DO I=1 , NUMCS HT(I,NUMPT(I)+1)=999 HB(I,NUMPT(I)+1)=HB(I,NUMPT(I)) END DO DO I=1 , NUMCS DO J=1 , NUMPT(I) XX(I,NUMPT(I)-J+1)= HB(I,J) XX(I,NUMPT(I)+J-1)=-HB(I,J) YY(I,NUMPT(I)-J+1)= HT(I,J) YY(I,NUMPT(I)+J-1)= HT(I,J) END DO END DO C C====== culculation of data for hull form ===== C======half breadth, water line area, and girth length WLHT=D DO I=1 , NUMCS SUM=0.0 10 CONTINUE IF( HT(I,1).GT.WLHT)THEN WLHB(I)=0. WLAR(I)=0. GLEN(I)=0. GOTO 10 END IF IF( HT(I,1).EQ.WLHT)THEN WLHB(I)=HB(I,1) WLAR(I)=0. GLEN(I)=HB(I,1) END IF DO J=1 , NUMPT(I) 11 CONTINUE IF( (HT(I,J).LT.WLHT) .AND. (WLHT.LE.HT(I,J+1)))THEN GOTO 12 END IF SUM=SUM+(HT(I,J+1)-HT(I,J))*(HB(I,J+1)+HB(I,J))/2 GOTO 14 C ====== calculation on each water line 12 CONTINUE AA=(HB(I,J+1)-HB(I,J))/(HT(I,J+1)-HT(I,J)) BB=(-HB(I,J+1)*HT(I,J)+HB(I,J)*HT(I,J+1))/(HT(I,J+1)-HT(I,J)) WLHB(I)=AA*WLHT+BB WLAR(I)=SUM+(WLHT-HT(I,J))*(WLHB(I)+HB(I,J))/2 GOTO 15 14 CONTINUE END DO 15 CONTINUE END DO C ===== Calculation of a section modulus DO I=1 , NUMCS IF( HT(1,1) .LT.D)THEN CSDMIN=1 END IF IF( (HT(I-1,1) .GT.D) .AND. (HT(I,1).GE.D))THEN GOTO21 END IF IF( (HT(I-1,1).GE.D) .AND. (HT(I,1).LT.D))THEN CSDMIN=I END IF IF( (HT(I,1).GE.D) .AND. (HT(I-1,1).LT.D))THEN NUMH0=I-1 GOTO 22 END IF HB2=HB(I,1) DO J=2 , NUMPT(I) IF(HT(I,J).GT.WLHT)THEN GOTO 18 END IF IF(HB(I,J).GT.HB2)THEN HB2=HB(I,J) END IF END DO 18 CONTINUE IF( WLHB(I).GT.HB2)THEN HB2=WLHB(I) END IF IF( WLHB(I).EQ.0)THEN BBYD(I)=HB2/(WLHT-HT(I,1)) GOTO 19 END IF BBYD(I)=WLHB(I)/(WLHT-HT(I,1)) 19 CONTINUE IF( WLHB(I).EQ.0)THEN SGM(I)=WLAR(I)/(HB2*(WLHT-HT(I,1))) GOTO 20 END IF SGM(I)=WLAR(I)/(WLHB(I)*(WLHT-HT(I,1))) 20 CONTINUE TT0(I)=(WLHT-HT(I,1))/D IF( TT0(I).LT.0)THEN TT0(I)=0 END IF IF( WLHB(I).EQ.0)THEN BB0(I)=HB2/B*2 GOTO 21 END IF BB0(I)=WLHB(I)/B*2 21 CONTINUE END DO NUMH0=NUMCS 22 CONTINUE C N=NUMCS-CSDMIN+1 DO I=1 , N X(I) =CSD(CSDMIN-1+I) HO(I) =BBYD(CSDMIN-1+I) SIG(I)=SGM(CSDMIN-1+I) BX(I) =BB0(CSDMIN-1+I)*B DX(I) =TT0(CSDMIN-1+I)*D END DO C C====== culculation of sidplacement ODD=0.0 if (mod(NUMCS,2).LT.0.5) then ODD=1.0 end if DO J=1 , NUMCS-2-ODD , 2 Q0=CSD(J) Q1=CSD(J+1) Q2=CSD(J+2) QQ=(Q0+Q2)/2 HH=(Q2-Q0)/2*L/10/3 L0=4*(QQ-Q1)*(QQ-Q2)/(Q0-Q1)/(Q0-Q2)+1 L1=4*(QQ-Q0)*(QQ-Q2)/(Q1-Q0)/(Q1-Q2) L2=4*(QQ-Q0)*(QQ-Q1)/(Q2-Q0)/(Q2-Q1)+1 SUM3=SUM3+(L0*WLAR(J)+L1*WLAR(J+1)+L2*WLAR(J+2))*HH END DO IF( ODD.LT.0.5)THEN GOTO 1404 END IF J=NUMCS-1 HH=(CSD(J+1)-CSD(J))/2*L/10 SUM3=SUM3+(WLAR(J)+WLAR(J+1))*HH 1404 NABLA=2*SUM3 C C====== Calculation of a block coefficient and an inside crossection coefficient CB=NABLA/L/B/D CM=2*WLAR(IMID)/B/D C C======Calculation of the encounter period which calculates a roll damping (an incident-wave period, Fn, angle of wave direction) DO K=1 , NKAI DO J=1 , NFN DO I=1 , NRA OMG=SQRT(2*PAI*G/RAM(I)/L) KW=2*PAI/RAM(I)/L V2=FUN(J)*SQRT(L*G) KAI1=KAI(K)*PAI/180 MGA=OMG-KW*V2*COS(KAI1) MGA=ABS(MGA) PERI(I,J,K)=2*PAI/MGA END DO END DO END DO C C===== Calculation of the attachment height of a bilge keel, and an index of a data number IF( TYPE_BK.NE.1)THEN DO I=SSBK1 , SSBK2 J=1 3 CONTINUE NUBK(I)=NUMPT(I)-J+1 IF( HB(I,J).LE.BKB(I))THEN J=J+1 GOTO 3 END IF BKT(I)=(HT(I,J)-HT(I,J-1))*(BKB(I)-HB(I,J-1)) & /(HB(I,J)-HB(I,J-1))+HT(I,J-1) END DO END IF C C********************************************************************* C OPEN(2,FILE='OSMB44.CSV') NO=0 DO ITHETA=1 , NTHETA THTA=THETA1(ITHETA)*3.14159/180 DO IKAI=1 , NKAI DO IFN=1 , NFN DO IRAM=1 , NRA NO=NO+1 C OMEGA=6.28318/PERI(IRAM,IFN,IKAI) C C *********** Frictional Component ************ C Ref. H.Kato (JZK. No.102) and S.Tamiya et al (JZK. No.132) C SF=L*(1.7*D+CB*B) RF=((0.887+0.145*CB)*(1.7*D+CB*B)-2*OGD*D)/3.1415 BFHAT=0.787*SF*RF**2.*SQRT(OMEGA*NUE*B/19.6133) & /(NABLA*B**2.)*(1.+4.1*FUN(IFN)/OMEGA*SQRT(9.80665/L)) C C C ********** Wave Making Component ********** C Ref. Y.Ikeda et al (JZK. No.143) C GUZAID=OMEGA**2.*D/9.80665 A1=1+GUZAID**(-1.2)*EXP(-2*GUZAID) A2=0.5+GUZAID**(-1)*EXP(-2*GUZAID) LOMEGA=OMEGA*FUN(IFN)*SQRT(L/9.80665) SIT=20.*(LOMEGA-0.3) TAH=(EXP(SIT)-EXP(-SIT))/(EXP(SIT)+EXP(-SIT)) BWHAT=BWOSM*.5*(((A2+1)+(A2-1)*TAH) & +(2*A1-A2-1)*EXP(-150*(LOMEGA-.25)**2)) C C C ********** Lift Component ********** C Ref. Y.Ikeda et al (JZK. No.143) C IF( CM.LE.0.92 )THEN KAPA=0.0 END IF IF( (CM.LE.0.97).AND.(CM.GT.0.92) )THEN KAPA=0.1 END IF IF( CM.GT.0.97 )THEN KAPA=0.3 END IF KN=6.28319*D/L+KAPA*(4.1*B/L-0.045) OG=OGD*D LO=0.3*D LR=0.5*D BLHAT=L*D*KN*LO*LR*FUN(IFN)*0.5/(NABLA*B**2) & *SQRT(0.5*L*B)*(1-1.4*OG/LR+0.7*OG**2/(LO*LR)) C C IF( RDP.EQ.2)THEN C C ******* Eddy Making Component for type barge ships ******* C Ref. Y.Ikeda et al(Proc. Ship Motion, Wave Loads and C Propaulsive Performance in a Seaway,1984) C Y.Ikeda, T.Fujiwara, T.Katayama (Proc ISOPE Vol.3, 1993) C OG=OGD*D DO I=1,N BE2(I)=(2/3.141592)*DX(I)**4/NABLA/BX(I)**2*(B/2/9.8)**.5 & *(HO(I)**2+1-OG/DX(I))*(HO(I)**2 & +(1-OG/DX(I))**2)*THTA*OMEGA END DO BEHAT=0 DO I=1,N-1 XHA=(X(I+1)-X(I)) SBE=XHA*(BE2(I)+BE2(I+1))/2 BEHAT=BEHAT+SBE END DO BEHAT=BEHAT*L/(X(N)-X(1)) IF (FUN(IFN).EQ.0) then BEHAT=BEHAT GOTO 5215 END IF AK=OMEGA/FUN(IFN)*SQRT(L/9.8) BEHAT=BEHAT*(.04*AK)**2/((.04*AK)**2+1) 5215 CONTINUE C C GOTO 5064 END IF C C *********** Eddy Making Component ********** C Ref. Y.Ikeda et al (JZK. No.142 , JZK. No.143) C Ref. T.Katayama et al (PJASNAOE No.33) C DO J=1 , N AHO=HO(J)/(1-OGD) SGUMA=(SIG(J)-OGD)/(1-OGD) E=(AHO-1)/(AHO+1) E2=E**2 A=4*SGUMA*(1-E2)/3.1415+E2 O=-A/(A+3) O2=SQRT(O**2-(A-1)/(A+3)) A3=O+O2 A1=E*(1+A3) AM=BX(J)/(1+A1+A3)*0.5 AA1=A1*(1+A3)/A3*0.25 IF( AA1.GT.1 )THEN AA1=1 END IF IF( AA1.LE.-1 )THEN AA1=-1 END IF DO I=1 , 2 CALL ARC(AA1,ACS) LTHETA=.5*ACS IF( I.EQ.1 )THEN LTHETA=0 END IF AH=1+A1**2+9*A3**2+2*A1*(1-3*A3)*COS(2*LTHETA) & -6*A3*COS(4*LTHETA) AA=-2*A3*COS(5*LTHETA)+A1*(1-A3)*COS(3*LTHETA) & +((6-3*A1)*A3**2+(A1**2-3*A1)*A3+A1**2)*COS(LTHETA) BB=-2*A3*SIN(5*LTHETA)+A1*(1-A3)*SIN(3*LTHETA) & +((6+3*A1)*A3**2+(3*A1+A1**2)*A3+A1**2)*SIN(LTHETA) V(I)=2*AM*SQRT(AA**2+BB**2)/AH RMAX1(I)=AM*SQRT(((1+A1)*SIN(LTHETA)-A3*SIN(3*LTHETA))**2 & +((1-A1)*COS(LTHETA)+A3*COS(3*LTHETA))**2) END DO RMAX=RMAX1(1) VMAX=V(1) IF( RMAX1(1).LE.RMAX1(2) )THEN GOTO 5159 END IF GOTO 5161 5159 RMAX=RMAX1(2) VMAX=V(2) 5161 RMEAN=2*DX(J)*(1-OGD)*SQRT(AHO*SGUMA/3.1415) P1=VMAX/RMEAN P2=RMAX/RMEAN PP3=P1+P2 GAMMA=(1+4*EXP(-165000*(1-SGUMA)**2))*PP3 CP=0.5*(0.87*EXP(-GAMMA)-4*EXP(-0.187*GAMMA)+3) SIT=20*(SIG(J)-0.7) CALL TANH(TAH,SIT) F1=0.5*(1+TAH) F2=0.5*(1-COS(3.1415*SIG(J)))-1.5*(1-EXP(-5*(1-SIG(J)))) & *SIN(3.1415*SIG(J))**2 IF( SIG(J).GT.1 )THEN SIG(J)=1. END IF R=2*DX(J)*SQRT(HO(J)*(SIG(J)-1)/(-0.8584)) RD=R/DX(J) IF( (HO(J).LE.1).AND.(RD.GE.AHO) )THEN R=.5*BX(J) END IF IF( (HO(J).GE.1).AND.(RD.GE.1) )THEN R=DX(J) END IF RD=R/DX(J) CR1(J)=RMAX**2/DX(J)**2*CP*((1-F1*RD)*(1-3/2*OGD-F1*RD) & +F2*(HO(J)-F1*RD)**2) IF( CR1(J).LT.0 )THEN CR1(J)=0 END IF END DO DO K=1 , 21 X1(K)=0+.5*(K-1) MAX=N MAXX=SC NN=MAX M1=1 M2=0 XXX=X1(K) DO JJ=1 , MAXX X11(JJ)=X(JJ) Y11(JJ)=CR1(JJ) END DO CALL HOKAN1(NN,IL,XXX,YX,Y2,M1,M2,X11,Y11,WX,WY) CR(K)=Y2 DAM=YX END DO CR(1)=1.5*(1-3/2*OGD) CR(21)=1.5*(1-3/2*OGD) SAM=0 DO K=1,10 K2=2*K K1=K2-1 K3=K2+1 SAM1=CR(K1)+4*CR(K2)+CR(K3) SAM=SAM+SAM1 END DO CRT=SAM/60 BEHAT=4*L*D**4/3/3.1415*OMEGA*SQRT(B/19.6)/NABLA/B**2*CRT*THTA DO ISSN=1,N BE2(ISSN)=CR1(ISSN)*4*D**4/3/3.1415*OMEGA & *SQRT(B/19.6)/NABLA/B**2*THTA END DO C IF (FUN(IFN).EQ.0) then BEHAT=BEHAT GOTO 5211 END IF C AK=OMEGA/FUN(IFN)*SQRT(L/9.8) BEHAT=BEHAT*(.04*AK)**2/((.04*AK)**2+1) DO ISSN=1,N BE2(ISSN)=BE2(ISSN)*(.04*AK)**2/((.04*AK)**2+1) END DO 5211 CONTINUE C C 5064 IF( BBK.LE.1E-10)THEN GOTO 5067 END IF IF( TYPE_BK.EQ.1)THEN C C ********** Damping due to Bilge Keels ********** C Ref. Y.Ikeda et al (JZK. No.161 , JZK. No.165) C XBK(1)=XBK1 XBK(11)=XBK2 DO I=2,10 XBK(I)=XBK(I-1)+(XBK2-XBK1)*.1 END DO MAX=N DO I=1,11 MAXX=SC NN=MAX M1=1 M2=0 XXX=XBK(I) DO JJ=1,MAXX X11(JJ)=X(JJ) Y11(JJ)=HO(JJ) END DO CALL HOKAN1(NN,IL,XXX,YX,Y2,M1,M2,X11,Y11,WX,WY) HO1=Y2 DAM=YX MAXX=SC NN=MAX M1=1 M2=0 XXX=XBK(I) DO JJ=1,MAXX X11(JJ)=X(JJ) Y11(JJ)=SIG(JJ) END DO CALL HOKAN1(NN,IL,XXX,YX,Y2,M1,M2,X11,Y11,WX,WY) SG1=Y2 DAM=YX MAXX=SC NN=MAX M1=1 M2=0 XXX=XBK(I) DO JJ=1,MAXX X11(JJ)=X(JJ) Y11(JJ)=DX(JJ) END DO CALL HOKAN1(NN,IL,XXX,YX,Y2,M1,M2,X11,Y11,WX,WY) DX1=Y2 DAM=YX MAXX=SC NN=MAX XXX=XBK(I) M1=1 M2=0 XXX=XBK(I) DO JJ=1,MAXX X11(JJ)=X(JJ) Y11(JJ)=BX(JJ) END DO CALL HOKAN1(NN,IL,XXX,YX,Y2,M1,M2,X11,Y11,WX,WY) BX1=Y2 DAM=YX R=2*DX1*sqrt(HO1*(SG1-1)/(-.8585)) RD=R/DX1 IF((HO1.LE.1).AND.(RD.GE.HO1)) then R=.5*BX1 END IF IF((HO1.GT.1).AND.(RD.GE.1)) then R=DX1 END IF RD=R/DX1 F=1+.3*EXP(-160*(1-SG1)) RBK=DX1*SQRT((HO1-.2929*RD)**2+(1-OGD-.2929*RD)**2) M1=RD M2=OGD M3=1-M1-M2 M4=HO1-M1 M5=(.414*HO1+.0651*M1**2-(.382*HO1+.0106)*M1) & /((HO1-.215*M1)*(1-.215*M1)) M6=(.414*HO1+.0651*M1**2-(.382+.0106*HO1)*M1) & /((HO1-.215*M1)*(1-.215*M1)) SO=.3*(3.1415*F*RBK*THTA)+1.95*BBK M7=SO/DX1-.25*3.1415*M1 R1=.25*3.1415*R IF(SO.LT.R1) then M7=0 END IF M8=M7+.414*M1 IF(SO.LT.R1) then M8=M7+1.414*(1-cos(SO/R))*M1 END IF A=(M3+M4)*M8-M7**2 BB=M4**3/3/(HO1-.215*M1)+(1-M1)**2*(2*M3-M2) & /6/(1-.215*M1)+M1*(M3*M5+M4*M6) CPPLAS=1.2 CPMINS=-22.5*BBK/(3.1415*RBK*F*THTA)-1.2 CD=CPPLAS-CPMINS C C *** BBKhat for unit length *** C RATIO(I)=RBK*BBK*CD/(RBK*BBK*CD+.5*DX1**2 & *(-A*CPMINS+BB*CPPLAS)) BBKHAT(I)=8*RBK**2*OMEGA*SQRT(B/19.6) & *THTA*F**2/(3*3.1415*NABLA*B**2) & *(RBK*BBK*CD+.5*DX1**2*(-A*CPMINS+BB*CPPLAS)) END DO C C *** BBKhat for Three Dimentional Ship Form *** C SAM=0 DO I=1,5 I2=2*I SAM1=BBKHAT(I2-1)+4*BBKHAT(I2)+BBKHAT(I2+1) SAM=SAM+SAM1 END DO BKHAT=SAM*(XBK2-XBK1)*.1/3*L*.1 C C ELSE C C ********** Damping due to Bilge Keels ********** C Ref. Y.Ikeda, T.Katayama et al (JKSNAJ. No.222) C BBK0=0 BBK1=0 BBK2=0 BBK3=0 RO=102 KG=(1-OGD)*D DO P=SSBK1,SSBK2 SGUMA=(SIG(P)-OGD)/(1-OGD) C ******* normal pressure component ******* C ******* F1,R,CD----->FBN K=NUBK(P) XX0=XX(P,K) XX1=BKB(P) XX2=XX(P,K+1) YY0=YY(P,K) YY1=BKT(P) YY2=YY(P,K+1) IF(XX2.EQ.XX1) THEN XX2=XX(P,K+2) YY2=YY(P,K+2) END IF LOW=102 F1=1+.3*EXP((-1)*160*(1-SGUMA)) R=SQRT((KG-YY1)**2+XX1**2) CD=22.5*BBK/(3.14159*F1*R*THTA)+2.4 FBN=8/3*RO*R**2*BBK*THTA*F1**2*CD*OMEGA/6.28 C ******* XA,YA ***** Quadric approximation of a hull surface inclination IF(XX0.EQ.XX2) THEN XA=0 YA=YY1 GOTO 3550 END IF IF(YY0.EQ.YY2) THEN XA=XX1 YA=KG GOTO 3550 END IF FD0=YY0/((XX0-XX1)*(XX0-XX2)) FD1=YY1/((XX1-XX0)*(XX1-XX2)) FD2=YY2/((XX2-XX0)*(XX2-XX1)) FD=2*XX1*(FD0+FD1+FD2)-(XX1+XX2) & *FD0-(XX0+XX2)*FD1-(XX0+XX1)*FD2 XA=(XX1+FD*(YY1-KG))/(1+FD**2) YA=(FD*XX1+FD**2*(YY1-KG))/(1+FD**2)+KG C ******* LB1,BN 3550 LB1=SQRT((XA-XX1)**2+(YA-YY1)**2) BN=FBN*LB1*2 C C ******* negative pressure component ******* C ******* S0,CPN----->FBSN S0=.3*3.14159*F1*R*THTA+1.95*BBK CPN=1.2-CD FBSN=8/3*RO*R**2*F1**2*THTA*CPN*OMEGA/6.28 C ******* starboard side C ******* XX0,YY0,XX1,YY1 MPN=0 LB2S0=0 OWARI=0 BSN=0 BSNS=0 BSNP=0 K=NUBK(P) XX0=BKB(P) YY0=BKT(P) XX1=XX(P,K+1) YY1=YY(P,K+1) GOTO 3950 3900 K=K+1 XX0=XX(P,K) YY0=YY(P,K) XX1=XX(P,K+1) YY1=YY(P,K+1) GOTO 3950 3920 XX1=(XX1-XX0)*(S0+LB2-LB2S0)/LB2+XX0 YY1=(YY1-YY0)*(S0+LB2-LB2S0)/LB2+YY0 C ******* XM,YM,XA,YA 3950 XM=(XX0+XX1)/2 YM=(YY0+YY1)/2 IF(XX0.EQ.XX1) THEN XA=0 YA=YM GOTO 4010 END IF IF(YY0.EQ.YY1) THEN XA=XM YA=KG GOTO 4010 END IF XA=(YM+(XX0-XX1)/(YY0-YY1)*XM-KG)*(XX0-XX1) & *(YY0-YY1)/((YY0-YY1)**2+(XX0-XX1)**2) YA=(YY0-YY1)/(XX0-XX1)*XA+KG C ******* LB2,LB3----->BSNS 4010 LB2=SQRT((XX0-XX1)**2+(YY0-YY1)**2) LB2S0=LB2S0+LB2 IF((OWARI.EQ.0).AND.(S0.LT.LB2S0)) THEN OWARI=1 GOTO 3920 END IF LB3=SQRT((XM-XA)**2+(YM-YA)**2) BSNS=BSNS+FBSN*LB2/LB3*(XA*(YM-YA)-(YA-KG)*(XM-XA)) IF(OWARI.EQ.1) THEN OWARI=0 GOTO 4110 END IF GOTO 3900 C ******* port side C ******* XX0,YY0,XX1,YY1 4110 LB2S0=0 K=2*NUMPT(P)-1-NUBK(P) XX0=-BKB(P) YY0=BKT(P) XX1=XX(P,K+1) YY1=YY(P,K+1) GOTO 4160 4140 K=K+1 XX0=XX(P,K) YY0=YY(P,K) XX1=XX(P,K+1) YY1=YY(P,K+1) 4160 IF(YY1.GT.D) THEN OWARI=1 XX1=(XX1-XX0)*(D-YY0)/(YY1-YY0)+XX0 YY1=D END IF GOTO 4210 4180 XX1=(XX1-XX0)*(S0+LB2-LB2S0)/LB2+XX0 YY1=(YY1-YY0)*(S0+LB2-LB2S0)/LB2+YY0 C ******* XM,YM,XA,YA 4210 XM=(XX0+XX1)/2 YM=(YY0+YY1)/2 IF(XX0.EQ.XX1) THEN XA=0 YA=YM GOTO 4270 END IF IF(YY0.EQ.YY1) THEN XA=XM YA=KG GOTO 4270 END IF XA=(YM+(XX0-XX1)/(YY0-YY1)*XM-KG)*(XX0-XX1) & *(YY0-YY1)/((YY0-YY1)**2+(XX0-XX1)**2) YA=(YY0-YY1)/(XX0-XX1)*XA+KG C ******* LB2,LB3----->BSNP 4270 LB2=SQRT((XX0-XX1)**2+(YY0-YY1)**2) LB2S0=LB2S0+LB2 IF((OWARI.EQ.0).AND.(S0.LT.LB2S0)) THEN OWARI=1 GOTO 4180 END IF LB3=SQRT((XM-XA)**2+(YM-YA)**2) BSNP=BSNP+FBSN*LB2/LB3*(XA*(YM-YA)-(YA-KG)*(XM-XA)) IF(OWARI.EQ.1) THEN OWARI=0 GOTO 4350 END IF GOTO 4140 4350 BSN=BSNS+BSNP C C ******* positive pressure component ******* C ******* CPP CPP=1.2 C ******* starboard side C ******* GLWBK GLWBK=0 GLBKK0=0 GLBKK1=0 BSP=0 BSPS=0 BSPP=0 DO K=1,NUBK(P) XX0=XX(P,K) YY0=YY(P,K) XX1=XX(P,K+1) YY1=YY(P,K+1) IF(YY1.GE.D) THEN goto 4520 END IF IF((YY1.LT.D).AND.(YY0.GE.D)) then XX0=(XX1-XX0)*(D-YY0)/(YY1-YY0)+XX0 YY0=D END IF C IF(K.EQ.NUBK(P)) THEN XX1=BKB(P) YY1=BKT(P) END IF GLWBK=GLWBK+SQRT((XX0-XX1)**2+(YY0-YY1)**2) 4520 END DO C ******* GLBKK0,GLBKK1,FBSP0,FBSP1 DO K=1,NUBK(P) XX0=XX(P,K) YY0=YY(P,K) XX1=XX(P,K+1) YY1=YY(P,K+1) IF(YY1.GE.D) then goto 4740 END IF IF ((YY1.LT.D).AND.(YY0.GE.D)) then XX0=(XX1-XX0)*(D-YY0)/(YY1-YY0)+XX0 YY0=D END IF IF(K.EQ.NUBK(P)) THEN XX1=BKB(P) YY1=BKT(P) END IF GLBKK0=GLBKK1 GLBKK1=GLBKK0+SQRT((XX0-XX1)**2+(YY0-YY1)**2) FBSP0=8/3*RO*R**2*THTA*F1**2*CPP*OMEGA/6.28*GLBKK0/GLWBK FBSP1=8/3*RO*R**2*THTA*F1**2*CPP*OMEGA/6.28*GLBKK1/GLWBK C ******* XM,YM,XA,YA,LB2,LB3----->BSP XM=(FBSP0+2*FBSP1)/(3*FBSP0+3*FBSP1)*(XX1-XX0)+XX0 YM=(FBSP0+2*FBSP1)/(3*FBSP0+3*FBSP1)*(YY1-YY0)+YY0 IF (XX0.EQ.XX1) THEN XA=0 YA=YM GOTO 4700 END IF IF(YY0.EQ.YY1) THEN XA=XM YA=KG GOTO 4700 END IF XA=(YM+(XX0-XX1)/(YY0-YY1)*XM-KG)*(XX0-XX1) & *(YY0-YY1)/((YY0-YY1)**2+(XX0-XX1)**2) YA=(YY0-YY1)/(XX0-XX1)*XA+KG 4700 LB2=SQRT((XX0-XX1)**2+(YY0-YY1)**2) LB3=SQRT((XM-XA)**2+(YM-YA)**2) BSPS=BSPS+(FBSP0+FBSP1)/2*LB2/LB3 & *(XA*(YM-YA)-(YA-KG)*(XM-XA)) 4740 END DO C ******* port side C ******* GLWBK GLWBK=0 GLBKK0=0 GLBKK1=0 OWARI=0 KSTA=NUMPT(P) KEND=2*NUMPT(P)-1-NUBK(P) DO K=KSTA,KEND XX0=XX(P,K) YY0=YY(P,K) XX1=XX(P,K+1) YY1=YY(P,K+1) IF(K.EQ.KEND) THEN XX1=-BKB(P) YY1=BKT(P) END IF GLWBK=GLWBK+SQRT((XX0-XX1)**2+(YY0-YY1)**2) END DO C ******* GLBKK0,GLBKK1,FBSP0,FBSP1 DO K=KSTA,KEND XX0=XX(P,K) YY0=YY(P,K) XX1=XX(P,K+1) YY1=YY(P,K+1) IF(K.EQ.KEND) THEN XX1=-BKB(P) YY1=BKT(P) END IF GLBKK0=GLBKK1 GLBKK1=GLBKK0+SQRT((XX0-XX1)**2+(YY0-YY1)**2) FBSP0=8/3*RO*R**2*THTA*F1**2*CPP*OMEGA/6.28*GLBKK0/GLWBK FBSP1=8/3*RO*R**2*THTA*F1**2*CPP*OMEGA/6.28*GLBKK1/GLWBK C ******* XM,YM,XA,YA,LB2,LB3----->BSP XM=(FBSP0+2*FBSP1)/(3*FBSP0+3*FBSP1)*(XX1-XX0)+XX0 YM=(FBSP0+2*FBSP1)/(3*FBSP0+3*FBSP1)*(YY1-YY0)+YY0 IF (XX0.EQ.XX1) THEN XA=0 YA=YM GOTO 4990 END IF IF (YY0.EQ.YY1) THEN XA=XM YA=KG GOTO 4990 END IF XA=(YM+(XX0-XX1)/(YY0-YY1)*XM-KG)*(XX0-XX1) & *(YY0-YY1)/((YY0-YY1)**2+(XX0-XX1)**2) YA=(YY0-YY1)/(XX0-XX1)*XA+KG 4990 LB2=SQRT((XX0-XX1)**2+(YY0-YY1)**2) LB3=SQRT((XM-XA)**2+(YM-YA)**2) BSPP=BSPP+(FBSP0+FBSP1)/2*LB2/LB3 & *(XA*(YM-YA)-(YA-KG)*(XM-XA)) END DO BSP=BSPS+BSPP C C******************************************************************** BBK0=BBK1 BBK1=BN+BSN+BSP BBK2=BBK0+BBK1 IF (P.EQ.SSBK1) then BBK3=BBK3+L*(CSD(P)-XBK1)/10*BBK1 GOTO 5140 END IF BBK3=BBK3+L*(CSD(P)-CSD(P-1))/10*BBK2/2 5140 IF (P.EQ.SSBK2) THEN BBKDAMP=BBK3+L*(XBK2-CSD(P))/10*BBK1 END IF C ******* END DO BKHAT=BBKDAMP/(RO*NABLA*B**2)*SQRT(B/19.6) C ******* C======================================= C C END IF B2KHAT=BKHAT 5067 IF( BBK.LE.1E-10)THEN B2KHAT=0. END IF BB44H=BFHAT+BWHAT+BLHAT+BEHAT+B2KHAT B44 =BB44H*(RO*NABLA*B**2)/SQRT(B/19.6) BF =BFHAT *(RO*NABLA*B**2)/SQRT(B/19.6) BW =BWHAT *(RO*NABLA*B**2)/SQRT(B/19.6) BL =BLHAT *(RO*NABLA*B**2)/SQRT(B/19.6) BE =BEHAT *(RO*NABLA*B**2)/SQRT(B/19.6) BBBK=B2KHAT*(RO*NABLA*B**2)/SQRT(B/19.6) C WRITE(*,100)KAI(IKAI),BB44H,BFHAT WRITE(2,100)THETA1(ITHETA),FUN(IFN),RAM(IRAM) & ,KAI(IKAI),BB44H,BFHAT,BWHAT,BLHAT,BEHAT & ,B2KHAT,B44,BF,BW,BL,BE,BBBK 100 FORMAT( 16(F20.9,',')) c ANSWER=0 c DO ISSN=1 , N-1 c ANSWER=ANSWER+(X(ISSN+1)-X(ISSN)) c & *(BE2(ISSN)+BE2(ISSN+1))/2./10.*L c END DO END DO END DO END DO END DO CLOSE(2) STOP END C C C C ************ Subroutine ************* C C ********** Lagrange 3 Points Interpolation ********** SUBROUTINE HOKAN1(NN,IL,XXX,YX,Y2,M1,M2,X11,Y11,WX,WY) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION NUE,L,KAI,L0,L1,L2,NABLA,KW,KAI1,LTHETA,MGA,KN & ,LO,LR,KG,LB1,LB2,LB3,LB2S0,M1,M2,M3,M4,M5,M6,M7,M8 INTEGER BWOSM,TYPE_BK,RDP,CSDMIN,SSBK1,SSBK2 PARAMETER SC=200,PT=200,NNRA=70,NNFN=20,NNKAI=35,NNRD=30 C DIMENSION X11(SC),Y11(SC),WX(SC),WY(SC) C C N1=NN-1 DO IL=2,N1 IF(XXX.LE.X11(IL)) THEN GOTO 5324 END IF END DO 5324 I1=IL-1 IF (XXX.GT.X11(N1)) THEN I1=NN-2 END IF I2=I1+2 DO IL=I1,I2 II=IL+1-I1 WX(II)=X11(IL) WY(II)=Y11(IL) END DO IF (M1.NE.1) THEN GOTO 5334 END IF CALL LAG3(XXX,Y2,X11,Y11,WX,WY) 5334 IF (M2.NE.1) THEN GOTO 5336 END IF YX=0 5336 RETURN END C C C SUBROUTINE LAG3(XXX,Y2,X11,Y11,WX,WY) IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER SC=200,PT=200,NNRA=70,NNFN=20,NNKAI=35,NNRD=30 DOUBLE PRECISION NUE,L,KAI,L0,L1,L2,NABLA,KW,KAI1,LTHETA,MGA,KN & ,LO,LR,KG,LB1,LB2,LB3,LB2S0,M1,M2,M3,M4,M5,M6,M7,M8 INTEGER BWOSM,TYPE_BK,RDP,CSDMIN,SSBK1,SSBK2 DIMENSION X11(SC),Y11(SC),WX(SC),WY(SC) C Y2=0 DO I3=1,3 W=1 Z=1 DO J3=1,3 IF(J3.EQ.I3) THEN GOTO 5350 END IF W=W*(XXX-WX(J3)) Z=Z*(WX(I3)-WX(J3)) 5350 END DO Y2=Y2+WY(I3)*W/Z END DO RETURN END C C C ********** arc cos ********** SUBROUTINE ARC(AA1,ACS) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C IF(ABS(AA1).EQ.1) THEN ACS=0 GOTO 5361 END IF ACS=-ATAN(AA1/SQRT(-AA1*AA1+1))+3.14159/2 5361 RETURN END C C C ********** hyperbolic tangent ********** SUBROUTINE TANH(TAH,SIT) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C TAH=(EXP(SIT)-EXP(-SIT))/(EXP(SIT)+EXP(-SIT)) RETURN END C C