''' ************************************************************************************ File : LAS_KN_1A.py File Purpose : Define the Class "LAS_KN" with hybrid technique to compute heeled section parameters Class Purpose : -> Compute and store KN values over a gien range of heel angles for each given draft -> Function to Compute Intersection of Water Line with Section Curves at various floating conditions -> Computations are done for Floating Conditions Defined by -> Draft -> Heel -> The range of Drafts and Heel are read from the File "Set_LAS_KN.dat" -> computations are performed for each draft and Heel combination -> Functions to write the Large angle Stability Results to File "Res_LAS_KN.dat" -> Plots and Saves KN curves and as image files of format ".png" App : PyNAT - Python Based Naval Architecture Tools Version : 2021.0 Paradigm : Object Based Framework Author : Praveen Kumar Ch *************************************************************************************** ''' import numpy as np from scipy.interpolate import interp1d from scipy.integrate import quad import matplotlib.pyplot as plt from HullGeometry import HullGeometry import math as math import os class LAS_KN_1A: # -------------------------------------------------------------------------------------- # The Constructor # PyNAT 2021 a : PraveenKumarCH # -------------------------------------------------------------------------------------- def __init__(self,fs,GivenShip,ep): self.techq="FAST" self.readFromFile(fs) self.compute(GivenShip,ep) self.writeToFile(fs) # -------------------------------------------------------------------------------------- # -------------------------------------------------------------------------------------- # The Method to Read the Data from File # PyNAT 2021 a : PraveenKumarCH # -------------------------------------------------------------------------------------- def readFromFile(self,fs): print("\n\t\tPyNAT: Reading File:\t...\\",fs.ikn_path) # Reading the Solver properties fp=open(fs.ikn_path,"rt") lines=fp.readlines() fp.close() nums=lines[0].split(",") self.nHeel=len(nums) self.rangeHeels=np.array(range(self.nHeel),dtype=np.float64) for i in range(self.nHeel): self.rangeHeels[i] =float(nums[i])*math.pi/180 # -------------------------------------------------------------------------------------- # -------------------------------------------------------------------------------------- # The Method to Compute KN Data # PyNAT 2021 a : PraveenKumarCH # -------------------------------------------------------------------------------------- def compute(self,GivenShip,ep): print("\n\t\tPyNAT: Computing KN Curves Data") # Finding the CL-WL intercept range for each heel angle self.nDraft=11 self.rangeDrafts=np.zeros((self.nHeel,self.nDraft)) for iha in range(self.nHeel): if(abs(self.rangeHeels[iha])%(math.pi/2)<=0.0001 and abs(self.rangeHeels[iha])>=0.0001): self.rangeHeels[iha]-=0.017453292 # Finding the Limits of Draft at CL dclmin=0 dclmax=0 for ist in range(GivenShip.nst): for ipt in range(GivenShip.npts[ist]): dcl=abs(GivenShip.y[ist][ipt]*math.tan(self.rangeHeels[iha])) db=-dcl+GivenShip.z[ist][ipt] da=dcl+GivenShip.z[ist][ipt] if(dbdclmax): dclmax=da delDraft=(dclmax-dclmin)/(self.nDraft-1) dclmin=dclmin+delDraft/2 dclmax=dclmax-delDraft/2 # creating a range of drafts self.rangeDrafts[iha,:]=np.linspace(dclmin,dclmax,self.nDraft,endpoint=True) # print(iha,self.rangeDrafts[iha,:]) self.VOL=np.zeros([self.nHeel,self.nDraft],dtype=np.float64) self.MASS=np.zeros([self.nHeel,self.nDraft],dtype=np.float64) self.LCB=np.zeros([self.nHeel,self.nDraft],dtype=np.float64) self.TCB=np.zeros([self.nHeel,self.nDraft],dtype=np.float64) self.VCB=np.zeros([self.nHeel,self.nDraft],dtype=np.float64) self.WPA=np.zeros([self.nHeel,self.nDraft],dtype=np.float64) self.WSA=np.zeros([self.nHeel,self.nDraft],dtype=np.float64) self.KN=np.zeros([self.nHeel,self.nDraft],dtype=np.float64) self.m=np.zeros([self.nHeel,self.nDraft],dtype=np.float64) self.c=np.zeros([self.nHeel,self.nDraft],dtype=np.float64) self.zN=np.zeros([self.nHeel,self.nDraft],dtype=np.float64) self.yN=np.zeros([self.nHeel,self.nDraft],dtype=np.float64) # Compute for each heel angle for iha in range(self.nHeel): # For each Draft in the draft range for idr in range(self.nDraft): zd=self.rangeDrafts[iha][idr] phi=self.rangeHeels[iha] BR=np.zeros(GivenShip.nst,dtype=np.float64) Ar=np.zeros(GivenShip.nst,dtype=np.float64) LM=np.zeros(GivenShip.nst,dtype=np.float64) TM=np.zeros(GivenShip.nst,dtype=np.float64) VM=np.zeros(GivenShip.nst,dtype=np.float64) CL=np.zeros(GivenShip.nst,dtype=np.float64) for ist in range(GivenShip.nst): pwli=self.computeIntersections(GivenShip.secCrv[ist],zd,phi) swli=self.computeIntersections(GivenShip.secCrv[ist],zd,-phi) # print("Station no =%5d Heel =%6.3f draft =%6.3f"%(ist,phi,zd),":\n") # print("Port Side Intersection Points:",pwli) # print("Starboard Side Intersection Points:",swli) # print("Praveen Kumar CH PyNAT : HAULT") # When WL intersects closed section at only two locations if (len(pwli)==1 and len(swli)==1): secData=GivenShip.secCrv[ist].sectionHydroStatics(pwli[0]) #print("The Section Prop at Draft Zp",pwli[0],"\n",secData) ht=zd-pwli[0] bt=abs(ht*math.tan(math.pi/2-(phi))) Atr=0.5*ht*bt hbw=math.sqrt(ht**2+bt**2) tl=Atr*bt/3 vl=Atr*(zd-ht*2/3) secData[0]=hbw secData[1]+=Atr secData[2]+=Atr*GivenShip.x[ist] secData[3]+=tl secData[4]+=vl #print("ht=%10.5f\tbt=%10.5f\tAtr=%10.5f\tHb=%10.5f\ttl=%10.5f\tvl=%10.5f\t"%(ht,bt,Atr,hbw,tl,vl)) secData2=GivenShip.secCrv[ist].sectionHydroStatics(swli[0]) #print("The Section Prop at Draft Zp",swli[0],"\n",secData2) ht=zd-swli[0] bt=abs(ht*math.tan(math.pi/2-(phi))) Atr=0.5*ht*bt hbw=math.sqrt(ht**2+bt**2) tl=Atr*bt/3 vl=Atr*(zd-ht*2/3) secData2[0]=hbw secData2[1]+=Atr secData2[2]+=Atr*GivenShip.x[ist] secData2[3]+=tl secData2[4]+=vl #print("ht=%10.5f\tbt=%10.5f\tAtr=%10.5f\tHb=%10.5f\ttl=%10.5f\tvl=%10.5f\t"%(ht,bt,Atr,hbw,tl,vl)) secData[0]+=secData2[0] # Breadth on WL secData[1]+=secData2[1] # Area secData[2]+=secData2[2] # LM secData[3]-=secData2[3] # TM secData[4]+=secData2[4] # VM secData[5]+=secData2[5] # CL if(abs(phi)>=math.pi/2): # In case if the ship is at a capsizing heel angle... fullSecData=GivenShip.secCrv[ist].sectionHydroStatics(GivenShip.secCrv[ist].Z[-1]) secData[1]=2*fullSecData[1]-secData[1] secData[2]=2*fullSecData[2]-secData[2] secData[3]=-secData[3] secData[4]=2*fullSecData[4]-secData[4] secData[5]=2*fullSecData[5]-secData[5] if secData[1]!=0: secData[6]=secData[3]/secData[1] # YC secData[7]=secData[4]/secData[1] # ZC else: secData[6]=0 # YC secData[7]=0 # ZC #print("Section Data Computed using Plane Geometry:\n",secData) else: if(phi<=0): side=1 else: side=-1 secData=self.compImmersedSectionProperties(GivenShip.secCrv[ist],zd,phi,side) # This portion of code for computation of Girth length of immersed section # This is improper and needs to be improved... # Since not a part of KN calculations so ignored for now secData[5]=0 if(len(pwli)!=0): secData[5]+=GivenShip.secCrv[ist].getGirthLength(pwli[-1]) if(len(swli)!=0): secData[5]+=GivenShip.secCrv[ist].getGirthLength(swli[-1]) # This portion of code for computation of section Breadth on water plane # This is improper and needs to be improved... # Since not a part of KN calculations so ignored for now secData[0]=0 pwli.extend(swli) for i in range(0,len(pwli)-1,2): dy=(pwli[i]-pwli[i+1])/abs(math.cos(phi-math.pi/2)) secData[0]+=dy #print("Section Data Computed using Numerical Integration:\n",secData) BR[ist]+=secData[0] Ar[ist]+=secData[1] LM[ist]+=secData[2] TM[ist]+=secData[3] VM[ist]+=secData[4] CL[ist]+=secData[5] #print("iSta=%5d, BRW=%10.5f, Ar=%10.5f, LM=%10.5f, TM=%10.5f, VM=%10.5f, CL=%10.5f"%(ist,BR[ist],Ar[ist],LM[ist],TM[ist],VM[ist],CL[ist])) # Integrating the Areas and Moments over the length Ar_fX=interp1d(GivenShip.x,Ar,kind='cubic') sq=quad(Ar_fX,GivenShip.x[0],GivenShip.x[-1]) self.VOL[iha,idr]=sq[0] self.MASS[iha,idr]=self.VOL[iha,idr]*ep.waterDensity Sa_fX=interp1d(GivenShip.x,CL,kind='cubic') sq=quad(Sa_fX,GivenShip.x[0],GivenShip.x[-1]) self.WSA[iha,idr]=sq[0] if self.VOL[iha,idr]!=0: lm_fX=interp1d(GivenShip.x,LM,kind='cubic') sq=quad(lm_fX,GivenShip.x[0],GivenShip.x[-1]) self.LCB[iha,idr]=sq[0]/self.VOL[iha,idr] tm_fX=interp1d(GivenShip.x,TM,kind='cubic') sq=quad(tm_fX,GivenShip.x[0],GivenShip.x[-1]) self.TCB[iha,idr]=sq[0]/self.VOL[iha,idr] vm_fX=interp1d(GivenShip.x,VM,kind='cubic') sq=quad(vm_fX,GivenShip.x[0],GivenShip.x[-1]) self.VCB[iha,idr]=sq[0]/self.VOL[iha,idr] else: self.LCB[iha,idr]=0 self.TCB[iha,idr]=0 self.VCB[iha,idr]=0 wpa_fX=interp1d(GivenShip.x,BR,kind='cubic') sq=quad(wpa_fX,GivenShip.x[0],GivenShip.x[-1]) self.WPA[iha,idr]=sq[0] # Finding KN self.m[iha,idr]=math.tan(math.pi/2-phi) self.c[iha,idr]=self.VCB[iha,idr]-self.m[iha,idr]*self.TCB[iha,idr] self.zN[iha,idr]=self.c[iha,idr]/(1+self.m[iha,idr]**2) self.yN[iha,idr]=-self.m[iha,idr]*self.zN[iha,idr] self.KN[iha,idr]=math.sqrt(self.zN[iha,idr]**2+self.yN[iha,idr]**2) print("\n\t\t\tVol. Disp=%12.6f\tHeel=%12.6f\tKN=%12.6f"%(self.VOL[iha,idr],phi*180/math.pi,self.KN[iha,idr])) print("\n\t\t\t\tyn=%10.5f \t zn=%10.5f"%(self.yN[iha,idr],self.zN[iha,idr])) # -------------------------------------------------------------------------------------- # -------------------------------------------------------------------------------------- # Function to compute heeled and immersed half section properties # PyNAT 2021 a : PraveenKumarCH # -------------------------------------------------------------------------------------- def compImmersedSectionProperties(self,sec,zd,ang,side): # Initialise return Variables with Default Values hBw=0 Ar=0 LM=0 VM=0 TM=0 CL=0 YC=0 ZC=0 dz=0.25 z=np.arange(sec.Z[0],sec.Z[-1],dz) da=np.arange(z.size,dtype=np.float64) dtm=np.arange(z.size,dtype=np.float64) dvm=np.arange(z.size,dtype=np.float64) # print("\nThe Heeled Section Calculations for %f"%(sec.X)) # print("Zd=",zd) # print("phi=",ang*180/math.pi) # print("side=",side) # PyNAT 2021 a : PraveenKumarCH if ang!=0: for i in range(z.size): yw=(z[i]-zd)*math.tan(ang-math.pi/2) ycp=sec.getHalfBreadth(z[i]) ycs=-ycp if (side>=0): # Measure strip ends Left Side of WL ya=max(ycp,yw) yb=max(ycs,yw) else: # Measure strip ends Right Side of WL ya=min(ycp,yw) yb=min(ycs,yw) da[i]=max(ya-yb,0) dtm[i]=da[i]*(ya+yb)/2 dvm[i]=da[i]*z[i] #print("yw=",yw,"yc=",yc,"ya=",ya,"yb=",yb,"da=",da[i]) a_fz=interp1d(z,da,kind='cubic') sq=quad(a_fz,z[0],z[-1]) Ar=sq[0] vm_fz=interp1d(z,dvm,kind='cubic') sq=quad(vm_fz,z[0],z[-1]) VM=sq[0] tm_fz=interp1d(z,dtm,kind='cubic') sq=quad(tm_fz,z[0],z[-1]) TM=sq[0] LM=Ar*sec.X if Ar!=0: YC=TM/Ar ZC=VM/Ar secData=list([hBw,Ar,LM,TM,VM,CL,YC,ZC]) else : secData=sec.sectionHydroStatics(zd) secData[0]=2*secData[0] secData[1]=2*secData[1] secData[2]=2*secData[2] secData[3]=0 secData[4]=2*secData[4] secData[5]=2*secData[5] secData[6]=0 if secData[0]!=0: secData[7]=secData[4]/secData[1] else: secData[7]=0 # print("SecData=",secData) return secData # -------------------------------------------------------------------------------------- # -------------------------------------------------------------------------------------- # The Method to predict water line intersection of a heeled section # PyNAT 2021 a : PraveenKumarCH # -------------------------------------------------------------------------------------- def computeIntersections(self,sec,zd,ang): # Tolerance to check intersection tol=0.001 # List of Intersection Z coordinates inPts=list() # If angle is zero i.e., zero heel angle if (math.tan(ang)==0): inPts.append(zd) return inPts # pre calculation of slope of water line tanang=math.tan(ang-math.pi/2) # Check for Water line cutting the bottom base line of section z1=sec.Z[0] yC=sec.getHalfBreadth(z1) yW=(z1-zd)*tanang if(yW<(yC-tol) and yW>=0): inPts.append(z1) # Check for water line cutting the section curve for i in range(sec.N-1): # take each pair of consecutive points z1=sec.Z[i] z2=sec.Z[i+1] dy1=sec.getHalfBreadth(z1)-(z1-zd)*tanang dy2=sec.getHalfBreadth(z2)-(z2-zd)*tanang #print("z1=%10.5f,yc1=%10.5f,yw1=%10.5f,dy1=%10.5f,z2=%10.5f,yc2=%10.5f,yw2=%10.5f,dy2=%10.5f"%(z1,yc1,yw1,dy1,z2,yc2,yw2,dy2)) if(dy1*dy2<0): for _ in range(25): if abs(dy1)<=tol: inPts.append(z1) #print("Appended z1",z1) break if abs(dy2)<=tol: inPts.append(z2) #print("Appended z2",z2) break zm=(z1+z2)/2 dym=sec.getHalfBreadth(zm)-(zm-zd)*tanang if(dy1*dym<0): z2=zm dy2=dym if(dy2*dym<0): z1=zm dy1=dym #print("z1=%10.5f,dy1=%10.5f,z2=%10.5f,dy2=%10.5f"%(z1,dy1,z2,dy2)) # Check for Water line cutting the top deck line z1=sec.Z[-1] yC=sec.getHalfBreadth(z1) yW=(z1-zd)*tanang if(yW<(yC-tol) and yW>=0): inPts.append(z1) # Return List of Points return inPts # -------------------------------------------------------------------------------------- # -------------------------------------------------------------------------------------- # The Method to Read the Data from File # PyNAT 2021 a : PraveenKumarCH # -------------------------------------------------------------------------------------- def writeToFile(self,fs): print("\n\t\tPyNAT: Writing KN Curves Data to File:\t ...\\",fs.okn_fast_path) if not os.path.exists(fs.las_dir): os.makedirs(fs.las_dir) fp=open(fs.okn_fast_path,"wt") fp.write('\nPy-NAT: Large angle Stability - KN Curves :\t') fp.write('\nDraft \tHeel \tKN \t') fp.write( 'zN \tyN \t') fp.write( 'm \tc \tWPA \t') fp.write( 'LCB \tTCB \tVCB \t') fp.write( 'Volume \tMass \tWet Area \t') for i in range(self.nHeel): for k in range(self.nDraft): # PyNAT 2021 a : PraveenKumarCH fp.write("\n%10.3E\t%10.3E\t%10.3E\t"%(self.rangeDrafts[i][k],self.rangeHeels[i],self.KN[i][k])) fp.write( "%10.3E\t%10.3E\t"%(self.zN[i][k],self.yN[i][k])) fp.write( "%10.3E\t%10.3E\t%10.3E\t"%(self.m[i][k],self.c[i][k],self.WPA[i][k])) fp.write( "%10.3E\t%10.3E\t%10.3E\t"%(self.LCB[i][k],self.TCB[i][k],self.VCB[i][k])) fp.write( "%10.3E\t%10.3E\t%10.3E\t"%(self.VOL[i][k],self.MASS[i][k],self.WSA[i][k])) fp.close() plt.figure("KN CURVES") for i in range(self.nHeel): hstr = format(self.rangeHeels[i]*180/math.pi,"0.2f") plt.plot(self.MASS[i,:],self.KN[i,:],label='KN@Heel='+hstr) plt.title("Cross Curves of Stability") plt.xlabel("Displacement") plt.ylabel("KN") lg=plt.legend(bbox_to_anchor=(1.05, 1)) plt.savefig(os.path.join(fs.las_dir,"LAS_KN_FAST.png"), dpi=300, format='png', bbox_extra_artists=(lg,), bbox_inches='tight') plt.close('all') # ----------------------------------------------------------------------------------------