diff --git a/Source_Script/ENV_PROP.py b/Source_Script/ENV_PROP.py new file mode 100644 index 0000000..0d1089f --- /dev/null +++ b/Source_Script/ENV_PROP.py @@ -0,0 +1,33 @@ +""" +************************************************************************************ + File : ENV_PROP.py + + File Purpose : Define the Class "ENV_PROP" + + Class Purpose : + -> Reads and stores Environment properties from file "Set_EP.dat" + -> Environment Properties are defined in this class + + App : PyNAT - Python Based Naval Architecture Tools + Version : 2021.0 + Paradigm : Object Based Framework + Author : Praveen Kumar Ch +************************************************************************************* +""" + +class ENV_PROP: + + # -------------------------------------------------------------------------------------- + # The Constructor + # -------------------------------------------------------------------------------------- + def __init__(self,fs): + + print("\n\t\tPyNAT: Reading File:\t...\\",fs.iep_path) + + # Reading the Solver properties + fp=open(fs.iep_path,"rt") + lines=fp.readlines() + fp.close() + + nums=lines[0].split('\t') + self.waterDensity = float(nums[0]) \ No newline at end of file diff --git a/Source_Script/FileSystem.py b/Source_Script/FileSystem.py new file mode 100644 index 0000000..4078561 --- /dev/null +++ b/Source_Script/FileSystem.py @@ -0,0 +1,85 @@ +""" +************************************************************************************ + File : HullGeometry.py + + File Purpose : Define the Class "FileSystems" + + Class Purpose : + -> Read the Hull Geometry definition data from the File "SetHullGeom.dat" + -> Create the data variables defining the ship Hull Geometry + -> Create and prepare a "SectionCurves" object for each station + -> Print the read Data to console + + App : PyNAT - Python Based Naval Architecture Tools + Version : 2021.0 + Paradigm : Object Based Framework + Author : Praveen Kumar Ch +************************************************************************************* +""" +import shutil +import os + +class FileSystem: + + # -------------------------------------------------------------------------------------- + # The Constructor + # -------------------------------------------------------------------------------------- + def __init__(self): + + # The Directories for input and output files + self.in_dir = "in" + self.out_dir = "out" + self.sec_dnm = "sec" + self.sec_dir =os.path.join(self.out_dir,self.sec_dnm) + self.hs_dnm="hs" + self.hs_dir =os.path.join(self.out_dir,self.hs_dnm) + self.las_dnm="las" + self.las_dir =os.path.join(self.out_dir,self.las_dnm) + + + # Check if Output Directory exists and delete if exist + # if(os.path.isdir(self.out_dir)): + # shutil.rmtree(self.out_dir) + + # Create the out directory + if not os.path.exists(self.out_dir): + os.makedirs(self.out_dir) + if not os.path.exists(self.sec_dir): + os.makedirs(self.sec_dir) + if not os.path.exists(self.hs_dir): + os.makedirs(self.hs_dir) + if not os.path.exists(self.las_dir): + os.makedirs(self.las_dir) + + self.iep_file = "Set_EP.dat" + self.iep_path = os.path.join(self.in_dir,self.iep_file) + + self.ihg_file = "Set_HG.dat" + self.ihg_path = os.path.join(self.in_dir,self.ihg_file) + + self.ihs_file = "Set_HS.dat" + self.ihs_path = os.path.join(self.in_dir,self.ihs_file) + + self.ikn_file = "Set_KN.dat" + self.ikn_path = os.path.join(self.in_dir,self.ikn_file) + + self.ilc_file = "Set_LC.dat" + self.ilc_path = os.path.join(self.in_dir,self.ilc_file) + + self.ohg_file = "Res_HG.dat" + self.ohg_path = os.path.join(self.sec_dir,self.ohg_file) + + self.ohs_file = "Res_HS.dat" + self.ohs_path = os.path.join(self.hs_dir,self.ohs_file) + + self.okn_file = "Res_KN.dat" + self.okn_path = os.path.join(self.las_dir,self.okn_file) + + self.ogz_file = "Res_GZ.dat" + self.ogz_path = os.path.join(self.las_dir,self.ogz_file) + + self.okn_fast_file = "Res_KN_FAST.dat" + self.okn_fast_path = os.path.join(self.las_dir,self.okn_fast_file) + + self.ogz_fast_file = "Res_GZ_FAST.dat" + self.ogz_fast_path = os.path.join(self.las_dir,self.ogz_fast_file) \ No newline at end of file diff --git a/Source_Script/HullGeometry.py b/Source_Script/HullGeometry.py new file mode 100644 index 0000000..fe35ada --- /dev/null +++ b/Source_Script/HullGeometry.py @@ -0,0 +1,363 @@ + +""" +************************************************************************************ + File : HullGeometry.py + + File Purpose : Define the Class "HullGeometry" + + Class Purpose : + -> Read the Hull Geometry definition data from the File "Set_HG.dat" + -> Create the data variables defining the ship Hull Geometry + -> Create and prepare a "SectionCurves" object for each station + -> Print the read Data to console + -> Write the prepared data to "Res_HG.dat" + + App : PyNAT - Python Based Naval Architecture Tools + Version : 2021a + Paradigm : Object Oriented Model + Author : Praveen Kumar Ch +************************************************************************************* +""" +import numpy as np +import scipy.interpolate as spi +import matplotlib.pyplot as plt +from SectionCurves import SectionCurves +import os + +class HullGeometry: + + # -------------------------------------------------------------------------------------- + # The Constructor + # Py-NAT 2021a : PRAVEEN KUMAR CH + # -------------------------------------------------------------------------------------- + def __init__(self,fs): + self.tol=0.001 + self.readFromFile(fs) + self.writeToFile(fs) + #self.printToConsol(fs) + # -------------------------------------------------------------------------------------- + + # -------------------------------------------------------------------------------------- + # Function to Read Hull Geometry Data from file + # Creates and Prepares SectionCurve Objects + # Py-NAT 2021a : PRAVEEN KUMAR CH + # -------------------------------------------------------------------------------------- + def readFromFile(self,fs): + + print("\n\t\t PyNAT: Reading File : ...\\",fs.ihg_path) + + fp=open(fs.ihg_path,"rt") + lines=fp.readlines() + fp.close() + + # Reading the Main Dimensions + nums=lines[1].split("\t") + self.L=float(nums[0]) + + nums=lines[2].split("\t") + self.B=float(nums[0]) + + nums=lines[3].split("\t") + self.D=float(nums[0]) + + nums=lines[4].split("\t") + self.AP=float(nums[0]) + self.FP=float(nums[1]) + self.MS=(self.AP+self.FP)/2 + + x=list() + y=list() + z=list() + c=list() + + nums=lines[5].rstrip("\n").split("\t") + tx=float(nums[0]) + ty=float(nums[1]) + tz=float(nums[2]) + tc=1 + if(len(nums)>=4): + if nums[3].count("k"): + tc=0 + if nums[3].count("b"): + tc=-1 + + x.append(tx) + y.append(list([ty])) + z.append(list([tz])) + c.append(list([tc])) + + # Reading the points and segregating into stations + # Py-NAT 2021a : PRAVEEN KUMAR CH + for i in range(6,len(lines),1): + + nums=lines[i].rstrip("\n").split("\t") + tx=float(nums[0]) + ty=float(nums[1]) + tz=float(nums[2]) + tc=1 + if(len(nums)>=4): + if nums[3].count("k"): + tc=0 + if nums[3].count("b"): + tc=-1 + + found=-1 + loc=0 + # check which station to add this point + for j in range(len(x)): + delx=tx-x[j] + if delx>0: + loc=j+1 + if abs(delx)<=self.tol: + found=j + loc=j + break + + if found<=-1: + # If the point dosent belong to exisitng station then add new station + x.insert(loc,tx) + y.insert(loc,list([ty])) + z.insert(loc,list([tz])) + c.insert(loc,list([tc])) + else: + # If the point belongs to existing station check were to place it + for j in range(len(z[loc]),0,-1): + delz=tz-z[loc][j-1] + if(delz>self.tol): + y[loc].insert(j,ty) + z[loc].insert(j,tz) + c[loc].insert(j,tc) + break + elif(abs(delz)<=self.tol): + dely=ty-y[loc][j-1] + if(dely>self.tol): + y[loc].insert(j,ty) + z[loc].insert(j,tz) + c[loc].insert(j,tc) + break + elif(abs(dely)<=self.tol): + break + + # Dealing with Horizontal Lines due to points with same z values + # Py-NAT 2021a : PRAVEEN KUMAR CH + for i in range(len(x)): + # In each station + for j in range(len(z[i])-1): + # At each point + if abs(z[i][j]-z[i][j+1])<=self.tol: + # If two adjacent points have same z + z[i][j+1]=z[i][j]+self.tol + + self.nst=len(x) + self.x=np.array(x) + self.y=list() + self.z=list() + self.c=list() + self.xms=np.array(range(self.nst)) + self.npts=np.array(range(self.nst)) + + for i in range(self.nst): + self.y.append(np.array(y[i])) + self.z.append(np.array(z[i])) + self.c.append(np.array(c[i])) + self.npts[i]=self.y[i].size + # Py-NAT 2021a : PRAVEEN KUMAR CH + print("\n\t\t PyNAT : Preparing Section Curve Spline Fit Data") + + self.secCrv=list() + for i in range(self.nst): + ind=np.argsort(self.z[i],axis=0) + self.z[i]=np.take_along_axis(self.z[i],ind,axis=0) + self.y[i]=np.take_along_axis(self.y[i],ind,axis=0) + self.c[i]=np.take_along_axis(self.c[i],ind,axis=0) + self.c[i][0]=-1 + self.c[i][-1]=-1 + self.secCrv.append(SectionCurves(self.x[i],self.y[i],self.z[i],self.c[i])) + + # Compute X locations form Mid Ship + self.xms=self.x-self.MS + + # --------------------------------------------------------------------------------------- + + # -------------------------------------------------------------------------------------- + # Function to write prepared hull geometry data to files + # Py-NAT 2021a : PRAVEEN KUMAR CH + # -------------------------------------------------------------------------------------- + def writeToFile(self, fs): + + print("\n\t\tPy-NAT: Writing Prepared Hull Geometry Data to File : ...\\",fs.ohg_path) + + if not os.path.exists(fs.sec_dir): + os.makedirs(fs.sec_dir) + + fp=open(fs.ohg_path,'wt') + fp.write("\nPy-NAT: Hull Geometry Data") + # Printing the read data to Console + fp.write("\nLength =\t%12.4E"%(self.L)) + fp.write("\nBreadth =\t%12.4E"%(self.B)) + fp.write("\nDepth =\t%12.4E"%(self.D)) + fp.write("\nX location of AP =\t%12.4E"%(self.AP)) + fp.write("\nX location of FP =\t%12.4E"%(self.FP)) + fp.write("\nX location of MS =\t%12.4E"%(self.MS)) + fp.write("\nNo of Stations =\t%12.4E"%(self.nst)) + fp.write("\n\n Printing Station Wise Data... \n\n") + for i in range(self.nst): + fp.write("\n-----------------------------------------------------------") + fp.write("\nStation no =\t%12.4E"%(i)) + fp.write("\nX =\t%12.4E"%(self.x[i])) + fp.write("\nX frm MS =\t%12.4E"%(self.xms[i])) + fp.write("\nNo of pts =\t%12.4E"%(self.npts[i])) + fp.write("\n i \t Z \t Y \t C") + for j in range(self.npts[i]): + fp.write("\n%4d\t%12.4E\t%12.4E\t%12d"%(j,self.z[i][j],self.y[i][j],self.c[i][j])) + + fp.write("\n\n Total No of Curvelets =%d"%(self.secCrv[i].nCrv)) + for k in range(self.secCrv[i].nCrv): + + fp.write("\n\n The Curve No: %d \t Start index= %10.5f \t End Index= %10.5f \n"%(k,self.secCrv[i].sCrv[k],self.secCrv[i].eCrv[k])) + fp.write("\nThe Spline curve Fit Y=F(z)") + tck=self.secCrv[i].y_fz[k] + fp.write("\nKnots, t=") + for l in range(len(tck[0])): + fp.write("%10.5f,"%(tck[0][l])) + fp.write("\nB-Spline Coeff, c=") + for l in range(len(tck[1])): + fp.write("%10.5f,"%(tck[1][l])) + fp.write("\nDegree, k=") + fp.write("%10.5f,"%(tck[2])) + + fp.write("\nThe Spline curve Fit ZY=F(z)") + tck=self.secCrv[i].zy_fz[k] + fp.write("\nKnots, t=") + for l in range(len(tck[0])): + fp.write("%10.5f,"%(tck[0][l])) + fp.write("\nB-Spline Coeff, c=") + for l in range(len(tck[1])): + fp.write("%10.5f,"%(tck[1][l])) + fp.write("\nDegree, k=") + fp.write("%10.5f,"%(tck[2])) + + fp.write("\nThe Spline curve Fit YY=F(z)") + tck=self.secCrv[i].yy_fz[k] + fp.write("\nKnots, t=") + for l in range(len(tck[0])): + fp.write("%10.5f,"%(tck[0][l])) + fp.write("\nB-Spline Coeff, c=") + for l in range(len(tck[1])): + fp.write("%10.5f,"%(tck[1][l])) + fp.write("\nDegree, k=") + fp.write("%10.5f,"%(tck[2])) + + fp.write("\nThe Spline curve Fit L=F(z)") + tck=self.secCrv[i].L_fz[k] + fp.write("\nKnots, t=") + for l in range(len(tck[0])): + fp.write("%10.5f,"%(tck[0][l])) + fp.write("\nB-Spline Coeff, c=") + for l in range(len(tck[1])): + fp.write("%10.5f,"%(tck[1][l])) + fp.write("\nDegree, k=") + fp.write("%10.5f,"%(tck[2])) + fp.close() + + #------------------------------------------------------------------------------------------------- + # Activate this Section to Plot Each individual Section + # Py-NAT 2021a : PRAVEEN KUMAR CH + # for i in range(self.nst): + # plt.figure(i) + # zs,ys=self.secCrv[i].getPoints(0.1) + # plt.plot(ys,zs) + # plt.plot(self.y[i],self.z[i],'o') + # xstr = format(self.x[i],"6.3f") + # istr = format(i,"3d") + # plt.title("Section @ X="+xstr) + # plt.xlabel("Y") + # plt.ylabel("Z") + # plt.grid(True) + # plt.axis('equal') + # #lg=plt.legend(bbox_to_anchor=(1.05, 1)) + # plt.savefig(os.path.join(fs.sec_dir,"Sec_Curv_"+istr+".png"), + # dpi=300, + # format='png', + # bbox_extra_artists=(), + # bbox_inches='tight') + # plt.close() + #------------------------------------------------------------------------------------------------- + + # Plot Section Curves from the Read Points + for i in range(self.nst): + plt.figure("Body Plan") + if(i Contains Function to Compute Hydrostatics + -> Computations are done for Floating Conditions Defined by + -> Draft + -> Trim + -> The range of Drafts and Trim are read from the File "Set_HS.dat" + -> Computations are performed over a range of drafts at all Individual Sections + -> Section wise drafts are computed for each Floating Condition + -> Computations are performed for each trim and draft combination Entire Hull Volume + -> Functions to write the sectional and Hull Hydrostatics to File "Res_HS.dat" + -> Plots and Saves Hydrostatics 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 +# ------------------------------------------------------------------------------------------ +# The Class to Compute and Store Hydrostatic Data +# Py-NAT 2021a : PRAVEEN KUMAR CH +# ------------------------------------------------------------------------------------------ + +class HydroStatics: + + # -------------------------------------------------------------------------------------- + # The Constructor + # Py-NAT 2021a : PRAVEEN KUMAR CH + # -------------------------------------------------------------------------------------- + def __init__(self,fs,shipGrid,ep): + self.readFromFile(fs,shipGrid) + self.computeSectionalHydroStatics(shipGrid) + self.computeHullHydroStatics(shipGrid,ep) + self.writeSectionalHydrostaticsToFile(fs,shipGrid) + self.writeHullHydroStaticsToFile(fs,shipGrid) + # -------------------------------------------------------------------------------------- + + # -------------------------------------------------------------------------------------- + # Read The Hydrostatic Solver Settings File + # create list of arrays of sectional drafts for each floating condition + # Py-NAT 2021a : PRAVEEN KUMAR CH + # -------------------------------------------------------------------------------------- + def readFromFile(self,fs,shipGrid): + + print("\n\t\tPyNAT : Reading File:\t...\\",fs.ihs_path) + # Reading the Solver properties + fp=open(fs.ihs_path,"rt") + lines=fp.readlines() + fp.close() + + nums=lines[0].split(",") + self.nDraft=len(nums) + self.rangeDrafts = np.array(range(self.nDraft),dtype=np.float64) + for i in range(self.nDraft): + self.rangeDrafts[i] =float(nums[i]) + + nums=lines[1].split(",") + self.nTrims=len(nums) + self.rangeTrims = np.array(range(self.nTrims),dtype=np.float64) + for i in range(self.nTrims): + self.rangeTrims[i] =float(nums[i])*math.pi/180 + + # Compute Sectional Drafts for each FLoating condition + # For each floating condition the section drafts are stored in an array + # A list of all such arrays is created + self.FC=list() + for j in range(self.nTrims): + temp=list() + for i in range(self.nDraft): + drafts=np.array(range(shipGrid.nst),dtype=np.float64) + for k in range(shipGrid.nst): + drafts[k]=self.rangeDrafts[i]-shipGrid.xms[k]*math.tan(self.rangeTrims[j]) + temp.append(drafts) + self.FC.append(temp) + self.nFC=self.nTrims*self.nDraft + + # -------------------------------------------------------------------------------------- + + # -------------------------------------------------------------------------------------- + # Function to compute sectional hydrostatics over a range of drafts + # Py-NAT 2021a : PRAVEEN KUMAR CH + # -------------------------------------------------------------------------------------- + def computeSectionalHydroStatics(self,shipGrid): + + print("\n\t\tPyNAT : Computing Section Wise Hydrostatic Properties") + + # Initializing the Arrays to Hold the Sectional Data + # Each Row is for a particular station + # Each Column is for a particular draft + self.Bw=np.zeros((shipGrid.nst, self.nDraft)) + self.Area=np.zeros((shipGrid.nst, self.nDraft)) + self.VM=np.zeros((shipGrid.nst, self.nDraft)) + self.TM=np.zeros((shipGrid.nst, self.nDraft)) + self.LM=np.zeros((shipGrid.nst, self.nDraft)) + self.CL=np.zeros((shipGrid.nst, self.nDraft)) + self.YC=np.zeros((shipGrid.nst, self.nDraft)) + self.ZC=np.zeros((shipGrid.nst, self.nDraft)) + + # Loop to Compute for Every section in the shipGrid + for i in range(shipGrid.nst): + print("\t\t\t\tStation:",i) + # Loop to compute for each draft in the range + for k in range(self.nDraft): + secData=shipGrid.secCrv[i].sectionHydroStatics(self.rangeDrafts[k]) + self.Bw[i][k]=2*secData[0] + self.Area[i][k]=2*secData[1] + self.LM[i][k]=2*secData[2] + self.TM[i][k]=secData[3] + self.VM[i][k]=2*secData[4] + self.CL[i][k]=2*secData[5] + self.YC[i][k]=secData[5] + self.ZC[i][k]=secData[7] + + # -------------------------------------------------------------------------------------- + + # -------------------------------------------------------------------------------------- + # Function to Compute Hull Volume Hydrostatics + # Py-NAT 2021a : PRAVEEN KUMAR CH + # -------------------------------------------------------------------------------------- + def computeHullHydroStatics(self,shipGrid,ep): + + print("\n\t\tPyNAT : Computing the Hull Hydrostatic Properties") + # Initializing the arrays to hold hydrostatics + self.VOL=np.zeros((self.nTrims,self.nDraft)) + self.WSA=np.zeros((self.nTrims,self.nDraft)) + self.WPA=np.zeros((self.nTrims,self.nDraft)) + self.LCF=np.zeros((self.nTrims,self.nDraft)) + self.MASS=np.zeros((self.nTrims,self.nDraft)) + self.VCB=np.zeros((self.nTrims,self.nDraft)) + self.TCB=np.zeros((self.nTrims,self.nDraft)) + self.LCB=np.zeros((self.nTrims,self.nDraft)) + self.IXX0=np.zeros((self.nTrims,self.nDraft)) + self.IXXG=np.zeros((self.nTrims,self.nDraft)) + self.IYY0=np.zeros((self.nTrims,self.nDraft)) + self.BMT=np.zeros((self.nTrims,self.nDraft)) + self.BML=np.zeros((self.nTrims,self.nDraft)) + self.KMT=np.zeros((self.nTrims,self.nDraft)) + self.KML=np.zeros((self.nTrims,self.nDraft)) + self.TPC=np.zeros((self.nTrims,self.nDraft)) + self.MCT=np.zeros((self.nTrims,self.nDraft)) + self.draftAP=np.zeros((self.nTrims,self.nDraft)) + self.draftFP=np.zeros((self.nTrims,self.nDraft)) + + # Loop for each Floating Condition [List of Array of Sectional Drafts] + for itr in range(self.nTrims): + for idr in range(self.nDraft): + # Py-NAT 2021a : PRAVEEN KUMAR CH + Area=np.zeros(shipGrid.nst) + VM=np.zeros(shipGrid.nst) + TM=np.zeros(shipGrid.nst) + LM=np.zeros(shipGrid.nst) + CL=np.zeros(shipGrid.nst) + YC=np.zeros(shipGrid.nst) + ZC=np.zeros(shipGrid.nst) + + y=np.arange(shipGrid.nst) + yx=np.arange(shipGrid.nst) + yxx=np.arange(shipGrid.nst) + yyy=np.arange(shipGrid.nst) + + for j in range(shipGrid.nst): + secData=shipGrid.secCrv[j].sectionHydroStatics(self.FC[itr][idr][j]) + y[j]=secData[0] + Area[j]=secData[1] + LM[j]=secData[2] + TM[j]=secData[3] + VM[j]=secData[4] + CL[j]=secData[5] + YC[j]=secData[6] + ZC[j]=secData[7] + yx[j]=y[j]*shipGrid.x[j] + yxx[j]=2*yx[j]*shipGrid.x[j] + yyy[j]=8*(y[j]**3)/12 + + Ar_fX=interp1d(shipGrid.x,Area,kind='cubic') + sq=quad(Ar_fX,shipGrid.x[0],shipGrid.x[-1]) + self.VOL[itr][idr]=2*sq[0] + self.MASS[itr][idr]=self.VOL[itr][idr]*ep.waterDensity + + Sa_fX=interp1d(shipGrid.x,CL,kind='cubic') + sq=quad(Sa_fX,shipGrid.x[0],shipGrid.x[-1]) + self.WSA[itr][idr]=2*sq[0] + + lm_fX=interp1d(shipGrid.x,LM,kind='cubic') + sq=quad(lm_fX,shipGrid.x[0],shipGrid.x[-1]) + self.LCB[itr][idr]=2*sq[0]/self.VOL[itr][idr] + + tm_fX=interp1d(shipGrid.x,TM-TM,kind='cubic') + sq=quad(tm_fX,shipGrid.x[0],shipGrid.x[-1]) + self.TCB[itr][idr]=sq[0]/self.VOL[itr][idr] + + vm_fX=interp1d(shipGrid.x,VM,kind='cubic') + sq=quad(vm_fX,shipGrid.x[0],shipGrid.x[-1]) + self.VCB[itr][idr]=2*sq[0]/self.VOL[itr][idr] + + wpa_fX=interp1d(shipGrid.x,y,kind='cubic') + sq=quad(wpa_fX,shipGrid.x[0],shipGrid.x[-1]) + self.WPA[itr][idr]=2*sq[0] + + wpm_fX=interp1d(shipGrid.x,yx,kind='cubic') + sq=quad(wpm_fX,shipGrid.x[0],shipGrid.x[-1]) + self.LCF[itr][idr]=2*sq[0]/self.WPA[itr][idr] + + IXX0_fx=interp1d(shipGrid.x,yxx,kind='cubic') + sq=quad(IXX0_fx,shipGrid.x[0],shipGrid.x[-1]) + self.IXX0[itr][idr]=sq[0] + + IYY0_fx=interp1d(shipGrid.x,yyy,kind='cubic') + sq=quad(IYY0_fx,shipGrid.x[0],shipGrid.x[-1]) + self.IYY0[itr][idr]=sq[0] + + self.IXXG[itr][idr]=self.IXX0[itr][idr]-self.WPA[itr][idr]*(self.LCF[itr][idr]**2) + + self.draftAP[itr][idr]=self.rangeDrafts[idr]+(shipGrid.MS-shipGrid.AP)*math.tan(self.rangeTrims[itr]) + self.draftFP[itr][idr]=self.rangeDrafts[idr]-(shipGrid.FP-shipGrid.MS)*math.tan(self.rangeTrims[itr]) + # Py-NAT 2021a : PRAVEEN KUMAR CH + print("\n\t\t\tTrim=%f\tDraft=%f\tVolume=%f"%(self.rangeTrims[itr],self.rangeDrafts[idr],self.VOL[itr][idr])) + + self.BML=self.IXXG/self.VOL + self.BMT=self.IYY0/self.VOL + self.KML=self.BML+self.VCB + self.KMT=self.BMT+self.VCB + self.TPC=self.WPA*ep.waterDensity/100000 + self.MCT=self.MASS*self.BML/(shipGrid.L*100) + # -------------------------------------------------------------------------------------- + + # -------------------------------------------------------------------------------------- + # Function to Write Section wise Hydrostatics to file + # Py-NAT 2021a : PRAVEEN KUMAR CH + # -------------------------------------------------------------------------------------- + def writeSectionalHydrostaticsToFile(self,fs,shipGrid): + + print("\n\t\tPyNAT : Writing Section wise Hydrostatic Data to File:\t...\\",fs.ohs_path) + + if not os.path.exists(fs.hs_dir): + os.makedirs(fs.hs_dir) + + fp=open(fs.ohs_path,"wt") + fp.write("\nPy-NAT: Hydrostatic Properties") + for i in range(shipGrid.nst): + fp.write('\n\n\nSection :\t'+str(i)) + fp.write('\nDraft \tB/2 \tArea \tLong. Mom \tVert. Mom \tTrans. Mom\tCurv Leng ') + for k in range(self.nDraft): + fp.write("\n%10.3E\t%10.3E\t%10.3E\t%10.3E\t%10.3E\t%10.3E\t%10.3E"%(self.rangeDrafts[k],self.Bw[i,k],self.Area[i,k],self.LM[i,k],self.VM[i,k],self.TM[i,k],self.CL[i,k])) + fp.close() + # Py-NAT 2021a : PRAVEEN KUMAR CH + for i in range(shipGrid.nst): + + plt.figure("Sec Ar vs Draft") + plt.plot(self.rangeDrafts,self.Area[i],label=str(i)) + plt.title("Sectional Areas") + plt.xlabel("Draft") + plt.ylabel("Area") + + plt.figure("Sec VM vs Draft") + plt.plot(self.rangeDrafts,self.VM[i],label=str(i)) + plt.title("Vertical Moment abt BL") + plt.xlabel("Draft") + plt.ylabel("Moment") + + plt.figure("Sec TM vs Draft") + plt.plot(self.rangeDrafts,self.TM[i],label=str(i)) + plt.title("Transverse Moment abt CL (for Half Section only)") + plt.xlabel("Draft") + plt.ylabel("Moment") + + plt.figure("Sec Girth vs Draft") + plt.plot(self.rangeDrafts,self.CL[i],label=str(i)) + plt.title("Girth Length") + plt.xlabel("Draft") + plt.ylabel("Girth") + + plt.figure("Sec Zc vs Draft") + plt.plot(self.rangeDrafts,self.ZC[i],label=str(i)) + plt.title("Z Centroid from BL ") + plt.xlabel("Draft") + plt.ylabel("Zc") + + plt.figure("Sec Yc vs Draft") + plt.plot(self.rangeDrafts,self.YC[i],label=str(i)) + plt.title("Y Centroid from CL (for Half Section only)") + plt.xlabel("Draft") + plt.ylabel("Yc") + + plt.figure("Sec Ar vs Draft") + lg=plt.legend(bbox_to_anchor=(1.05, 1)) + plt.savefig(os.path.join(fs.hs_dir,'Sect_Areas.png'), + dpi=300, + format='png', + bbox_extra_artists=(lg,), + bbox_inches='tight') + + plt.figure("Sec VM vs Draft") + lg=plt.legend(bbox_to_anchor=(1.05, 1)) + plt.savefig(os.path.join(fs.hs_dir,'Sect_VertMom.png'), + dpi=300, + format='png', + bbox_extra_artists=(lg,), + bbox_inches='tight') + + plt.figure("Sec TM vs Draft") + lg=plt.legend(bbox_to_anchor=(1.05, 1)) + plt.savefig(os.path.join(fs.hs_dir,'Sect_TransMom.png'), + dpi=300, + format='png', + bbox_extra_artists=(lg,), + bbox_inches='tight') + + plt.figure("Sec Girth vs Draft") + lg=plt.legend(bbox_to_anchor=(1.05, 1)) + plt.savefig(os.path.join(fs.hs_dir,'Sect_Girth.png'), + dpi=300, + format='png', + bbox_extra_artists=(lg,), + bbox_inches='tight') + + plt.figure("Sec Zc vs Draft") + lg=plt.legend(bbox_to_anchor=(1.05, 1)) + plt.savefig(os.path.join(fs.hs_dir,'Sect_ZCent.png'), + dpi=300, + format='png', + bbox_extra_artists=(lg,), + bbox_inches='tight') + + plt.figure("Sec Yc vs Draft") + lg=plt.legend(bbox_to_anchor=(1.05, 1)) + plt.savefig(os.path.join(fs.hs_dir,'Sect_YCent.png'), + dpi=300, + format='png', + bbox_extra_artists=(lg,), + bbox_inches='tight') + plt.close('all') + # -------------------------------------------------------------------------------------- + + # -------------------------------------------------------------------------------------- + # Function to Write HULL Hydrostatics to file + # Py-NAT 2021a : PRAVEEN KUMAR CH + # -------------------------------------------------------------------------------------- + def writeHullHydroStaticsToFile(self,fs,shipGrid): + # Py-NAT 2021a : PRAVEEN KUMAR CH + print("\n\t\tPyNAT : Writing Hull Hydrostatic Data to File:\t...\\",fs.ohs_path) + fp=open(fs.ohs_path,"at") + fp.write('\n\n\nHull Hydrostatics :\t') + fp.write('\nDraft MS \tTrim \t') + fp.write(' Draft AP \tDraft FP \t') + fp.write( 'Volume \tMass \tWet Hull Ar.\t') + fp.write( 'LCB \tTCB \tVCB \t') + fp.write( 'WPA \tLCF \t') + fp.write( 'IXX0 \tIXXG \tIYY0 \t') + fp.write( 'BM_L \tBM_T \tKM_L \tKM_T \t') + fp.write( 'TPC \tMCT \t') + fp.write( '') + + for i in range(self.nTrims): + for k in range(self.nDraft): + fp.write("\n%10.3E\t%10.3E\t"%(self.rangeDrafts[k],self.rangeTrims[i])) + fp.write( "%10.3E\t%10.3E\t"%(self.draftAP[i][k],self.draftFP[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.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"%(self.WPA[i][k],self.LCF[i][k])) + fp.write( "%10.3E\t%10.3E\t%10.3E\t"%(self.IXX0[i][k],self.IXXG[i][k],self.IYY0[i][k])) + fp.write( "%10.3E\t%10.3E\t%10.3E\t%10.3E\t"%(self.BML[i][k],self.BMT[i][k],self.KML[i][k],self.KMT[i][k])) + fp.write( "%10.3E\t%10.3E\t"%(self.TPC[i][k],self.MCT[i][k])) + + fp.close() + + plt.figure("Vol Disp vs Draft") + for i in range(self.nTrims): + hstr = format(self.rangeTrims[i],"6.4f") + plt.plot(self.rangeDrafts,self.VOL[i],label='Volume@Trim='+hstr) + plt.title("Vol Disp vs Draft") + plt.xlabel("Draft@MS") + plt.ylabel("Volume") + lg=plt.legend(bbox_to_anchor=(1.05, 1)) + plt.savefig(os.path.join(fs.hs_dir,"HS_VolDisp.png"), + dpi=300, + format='png', + bbox_extra_artists=(lg,), + bbox_inches='tight') + + plt.figure("Mass Disp vs Draft") + for i in range(self.nTrims): + hstr = format(self.rangeTrims[i],"6.4f") + plt.plot(self.rangeDrafts,self.MASS[i],label='Mass@Trim='+hstr) + plt.title("Mass Disp vs Draft") + plt.xlabel("Draft@MS") + plt.ylabel("Mass") + lg=plt.legend(bbox_to_anchor=(1.05, 1)) + plt.savefig(os.path.join(fs.hs_dir,"HS_MassDisp.png"), + dpi=300, + format='png', + bbox_extra_artists=(lg,), + bbox_inches='tight') + + plt.figure("Wet Area vs Draft") + for i in range(self.nTrims): + hstr = format(self.rangeTrims[i],"6.4f") + plt.plot(self.rangeDrafts,self.WSA[i],label='WSA@Trim='+hstr) + plt.title("Wetted Surface Area vs Draft") + plt.xlabel("Draft@MS") + plt.ylabel("Wet Surf Area") + lg=plt.legend(bbox_to_anchor=(1.05, 1)) + plt.savefig(os.path.join(fs.hs_dir,"HS_WetArea.png"), + dpi=300, + format='png', + bbox_extra_artists=(lg,), + bbox_inches='tight') + + plt.figure("WPA vs Draft") + for i in range(self.nTrims): + hstr = format(self.rangeTrims[i],"6.4f") + plt.plot(self.rangeDrafts,self.WPA[i],label='Aw@Trim='+hstr) + plt.title("Water Plane Area vs Draft") + plt.xlabel("Draft@MS") + plt.ylabel("WPA") + lg=plt.legend(bbox_to_anchor=(1.05, 1)) + plt.savefig(os.path.join(fs.hs_dir,"HS_WPA.png"), + dpi=300, + format='png', + bbox_extra_artists=(lg,), + bbox_inches='tight') + + plt.figure("LCB vs Draft") + for i in range(self.nTrims): + hstr = format(self.rangeTrims[i],"6.4f") + plt.plot(self.rangeDrafts,self.LCB[i]-shipGrid.MS,label='LCB@Trim='+hstr) + plt.title("LCB (Frm Mid Ship) vs Draft") + plt.xlabel("Draft@MS") + plt.ylabel("LCB") + lg=plt.legend(bbox_to_anchor=(1.05, 1)) + plt.savefig(os.path.join(fs.hs_dir,"HS_LCB.png"), + dpi=300, + format='png', + bbox_extra_artists=(lg,), + bbox_inches='tight') + + plt.figure("LCF vs Draft") + for i in range(self.nTrims): + hstr = format(self.rangeTrims[i],"6.4f") + plt.plot(self.rangeDrafts,self.LCF[i]-shipGrid.MS,label='LCF@Trim='+hstr) + plt.title("LCF (Frm Mid Ship) vs Draft") + plt.xlabel("Draft@MS") + plt.ylabel("LCF") + lg=plt.legend(bbox_to_anchor=(1.05, 1)) + plt.savefig(os.path.join(fs.hs_dir,"HS_LCF.png"), + dpi=300, + format='png', + bbox_extra_artists=(lg,), + bbox_inches='tight') + + plt.figure("VCB vs Draft") + for i in range(self.nTrims): + hstr = format(self.rangeTrims[i],"6.4f") + plt.plot(self.rangeDrafts,self.VCB[i],label='VCB@Trim='+hstr) + plt.title("VCB (Frm Base) vs Draft") + plt.xlabel("Draft@MS") + plt.ylabel("VCB") + lg=plt.legend(bbox_to_anchor=(1.05, 1)) + plt.savefig(os.path.join(fs.hs_dir,"HS_VCB.png"), + dpi=300, + format='png', + bbox_extra_artists=(lg,), + bbox_inches='tight') + + plt.figure("TCB vs Draft") + for i in range(self.nTrims): + hstr = format(self.rangeTrims[i],"6.4f") + plt.plot(self.rangeDrafts,self.TCB[i],label='TCB@Trim='+hstr) + plt.title("TCB From CL vs Draft") + plt.xlabel("Draft@MS") + plt.ylabel("TCB") + lg=plt.legend(bbox_to_anchor=(1.05, 1)) + plt.savefig(os.path.join(fs.hs_dir,"HS_TCB.png"), + dpi=300, + format='png', + bbox_extra_artists=(lg,), + bbox_inches='tight') + + plt.figure("TPC vs Draft") + for i in range(self.nTrims): + hstr = format(self.rangeTrims[i],"6.4f") + plt.plot(self.rangeDrafts,self.TPC[i],label='TPC@Trim='+hstr) + plt.title("Tons Per CM immersion vs Draft") + plt.xlabel("Draft@MS") + plt.ylabel("TPC") + lg=plt.legend(bbox_to_anchor=(1.05, 1)) + plt.savefig(os.path.join(fs.hs_dir,"HS_TPC.png"), + dpi=300, + format='png', + bbox_extra_artists=(lg,), + bbox_inches='tight') + + plt.figure("MCT vs Draft") + for i in range(self.nTrims): + hstr = format(self.rangeTrims[i],"6.4f") + plt.plot(self.rangeDrafts,self.MCT[i],label='MCT@Trim='+hstr) + plt.title("Mom. to Change Trim vs Draft") + plt.xlabel("Draft@MS") + plt.ylabel("MCT") + lg=plt.legend(bbox_to_anchor=(1.05, 1)) + plt.savefig(os.path.join(fs.hs_dir,"HS_MCT.png"), + dpi=300, + format='png', + bbox_extra_artists=(lg,), + bbox_inches='tight') + + plt.figure("KMT vs Draft") + for i in range(self.nTrims): + hstr = format(self.rangeTrims[i],"6.4f") + plt.plot(self.rangeDrafts,self.KMT[i],label='KMT@Trim='+hstr) + plt.title("KM Trans.") + plt.xlabel("Draft@MS") + plt.ylabel("KM_T") + lg=plt.legend(bbox_to_anchor=(1.05, 1)) + plt.savefig(os.path.join(fs.hs_dir,"HS_KMT.png"), + dpi=300, + format='png', + bbox_extra_artists=(lg,), + bbox_inches='tight') + + plt.figure("KML vs Draft") + for i in range(self.nTrims): + hstr = format(self.rangeTrims[i],"6.4f") + plt.plot(self.rangeDrafts,self.KML[i],label='KML@Trim='+hstr) + plt.title("KM Long.") + plt.xlabel("Draft@MS") + plt.ylabel("KM_L") + lg=plt.legend(bbox_to_anchor=(1.05, 1)) + plt.savefig(os.path.join(fs.hs_dir,"HS_KML.png"), + dpi=300, + format='png', + bbox_extra_artists=(lg,), + bbox_inches='tight') + + plt.figure("BMT vs Draft") + for i in range(self.nTrims): + hstr = format(self.rangeTrims[i],"6.4f") + plt.plot(self.rangeDrafts,self.BMT[i],label='BMT@Trim='+hstr) + plt.title("Trans. Metacentric Radius") + plt.xlabel("Draft@MS") + plt.ylabel("BM_T") + lg=plt.legend(bbox_to_anchor=(1.05, 1)) + plt.savefig(os.path.join(fs.hs_dir,"HS_BMT.png"), + dpi=300, + format='png', + bbox_extra_artists=(lg,), + bbox_inches='tight') + + plt.figure("BML vs Draft") + for i in range(self.nTrims): + hstr = format(self.rangeTrims[i],"6.4f") + plt.plot(self.rangeDrafts,self.BML[i],label='BML@Trim='+hstr) + plt.title("Long. Metacentric Radius") + plt.xlabel("Draft@MS") + plt.ylabel("BM_L") + lg=plt.legend(bbox_to_anchor=(1.05, 1)) + plt.savefig(os.path.join(fs.hs_dir,"HS_BML.png"), + dpi=300, + format='png', + bbox_extra_artists=(lg,), + bbox_inches='tight') + + plt.close('all') + # -------------------------------------------------------------------------------------- + # Py-NAT 2021a : PRAVEEN KUMAR CH +# ------------------------------------------------------------------------------------------ \ No newline at end of file diff --git a/Source_Script/LAS_GZ.py b/Source_Script/LAS_GZ.py new file mode 100644 index 0000000..cd8316c --- /dev/null +++ b/Source_Script/LAS_GZ.py @@ -0,0 +1,136 @@ +''' +************************************************************************************ + File : LAS_GZ.py + + File Purpose : Define the Class "LAS_GZ" + + Class Purpose : + -> Read the Load Case Data from the File set_LC.dat + -> Interpolate KN values at various heel angles for the given displacement + -> Compute GZ values + -> write GZ data to file + -> Plot GZ data in curves + + App : PyNAT - Python Based Naval Architecture Tools + Version : 2021.0 + Paradigm : Object Based Framework + Author : Praveen Kumar Ch +*************************************************************************************** +''' +import numpy as np +from scipy import interpolate +import matplotlib.pyplot as plt +from HullGeometry import HullGeometry +import math as math +import os +# Py-NAT 2021a : PRAVEEN KUMAR CH +class LAS_GZ: + + # -------------------------------------------------------------------------------------- + # The Constructor + # Py-NAT 2021a : PRAVEEN KUMAR CH + # -------------------------------------------------------------------------------------- + def __init__(self,fs,GivenShip,knData,ep): + self.readFromFile(fs) + self.compute(GivenShip,knData,ep) + self.writeToFile(fs,knData) + + # -------------------------------------------------------------------------------------- + # Function to Read Load Case Data from File + # Py-NAT 2021a : PRAVEEN KUMAR CH + # -------------------------------------------------------------------------------------- + def readFromFile(self,fs): + + print("\t\tReading Load Condition from file...") + # Reading the Solver properties + fp=open(fs.ilc_path,"rt") + lines=fp.readlines() + fp.close() + + self.nCases=len(lines) + self.Mass=list() + self.KG=list() + + for i in range(self.nCases): + nums=lines[i].split("\t") + self.Mass.append(float(nums[0].rstrip("\n"))) + self.KG.append(float(nums[1].rstrip("\n"))) + # -------------------------------------------------------------------------------------- + + + # -------------------------------------------------------------------------------------- + # Function to Compute GZ Data + # Py-NAT 2021a : PRAVEEN KUMAR CH + # -------------------------------------------------------------------------------------- + + def compute(self,GivenShip,knData,ep): + + print("\n\t\tComputing GZ Data") + + # Interpolate KN at each Heel + self.KN=np.ones((self.nCases,knData.nHeel),dtype=np.float64) + self.GZ=np.ones((self.nCases,knData.nHeel),dtype=np.float64) + # Py-NAT 2021a : PRAVEEN KUMAR CH + for iha in range(knData.nHeel): + x=np.copy(knData.MASS[iha,:]) + y=np.copy(knData.KN[iha,:]) + ind=np.argsort(x,axis=0) + x=np.take_along_axis(x,ind,axis=0) + y=np.take_along_axis(y,ind,axis=0) + kn_fMass=interpolate.splrep(x,y,k=3) + + for ics in range(self.nCases): + self.KN[ics][iha]=interpolate.splev(self.Mass[ics],kn_fMass,der=0) + phi=abs(knData.rangeHeels[iha]) + self.GZ[ics][iha]=self.KN[ics][iha]-self.KG[ics]*math.sin(phi) + #print("\n\t\t\tKN=",self.KN[ics][iha]) + print("\n\t\t\tKG= %10.4f Disp=%10.4f Heel=%6.4f GZ=%10.4f"%(self.KG[ics],self.Mass[ics],knData.rangeHeels[iha],self.GZ[ics][iha])) + # -------------------------------------------------------------------------------------- + + + # -------------------------------------------------------------------------------------- + # Function to Write GZ Data to File + # Py-NAT 2021a : PRAVEEN KUMAR CH + # -------------------------------------------------------------------------------------- + def writeToFile(self,fs,knData): + + print("\n\t\tPyNAT: Writing GZ Data to File:\t...\\",fs.ogz_path) + + if not os.path.exists(fs.las_dir): + os.makedirs(fs.las_dir) + # Py-NAT 2021a : PRAVEEN KUMAR CH + if(knData.techq=="DIRECT"): + filepath=fs.ogz_path + plotpath=os.path.join(fs.las_dir,"LAS_GZ.png") + else: + filepath=fs.ogz_fast_path + plotpath=os.path.join(fs.las_dir,"LAS_GZ_FAST.png") + + fp=open(filepath,'wt') + fp.write("\nPyNAT: Righting Lever (GZ) Curve Data") + + for ics in range(self.nCases): + fp.write("\nMass Displacement=%f"%(self.Mass[ics])) + fp.write("\nVCG (KG) = %f"%(self.KG[ics])) + fp.write('\nHeel Angle\tGZ \t') + + for i in range(knData.nHeel): + fp.write("\n%10.3E\t%10.3E\t"%(knData.rangeHeels[i]*180/math.pi,self.GZ[ics][i])) + + plt.figure("GZ Curve") + hst=format(ics,"2d") + plt.plot(knData.rangeHeels*180/math.pi,self.GZ[ics,:],label="Case No"+hst) + plt.title("Righting Lever (GZ) Curve") + plt.xlabel("Heel Angle") + plt.ylabel("GZ") + lg=plt.legend(bbox_to_anchor=(1.05, 1)) + plt.savefig(plotpath, + dpi=300, + format='png', + bbox_extra_artists=(lg,), + bbox_inches='tight') + + fp.close() + plt.close('all') + # Py-NAT 2021a : PRAVEEN KUMAR CH + # -------------------------------------------------------------------------------------- \ No newline at end of file diff --git a/Source_Script/LAS_KN_0B.py b/Source_Script/LAS_KN_0B.py new file mode 100644 index 0000000..8e6d0d5 --- /dev/null +++ b/Source_Script/LAS_KN_0B.py @@ -0,0 +1,332 @@ +''' +************************************************************************************ + File : LAS_KN.py + + File Purpose : Define the Class "LAS_KN" + with Direct Integration technique to compute heeled section prop. + + 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 +# Py-NAT 2021a : PRAVEEN KUMAR CH +class LAS_KN_0B: + + # -------------------------------------------------------------------------------------- + # The Constructor + # Py-NAT 2021a : PRAVEEN KUMAR CH + # -------------------------------------------------------------------------------------- + def __init__(self,fs,GivenShip,ep): + self.techq="DIRECT" + self.readFromFile(fs) + self.compute(GivenShip,ep) + self.writeToFile(fs) + # -------------------------------------------------------------------------------------- + + + # -------------------------------------------------------------------------------------- + # The Method to Read the Data from File + # Py-NAT 2021a : PRAVEEN KUMAR CH + # -------------------------------------------------------------------------------------- + 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 + # Py-NAT 2021a : PRAVEEN KUMAR CH + # -------------------------------------------------------------------------------------- + + 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(self.rangeDrafts) + + 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 the Water Line Intersections for each heel angle + for iha in range(self.nHeel): + + # For each Draft + 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) + + #print("------\t------\t------\t------\t------\t------\t------\t------\t------\t------\t") + #print("\t\tZd=%12.6f\tphi=%12.6f"%(zd,phi*180/math.pi)) + # Py-NAT 2 021a : PRAVEEN KUMAR CH + #print("------\t------\t------\t------\t------\t------\t------\t------\t------\t------\t") + + for ist in range(GivenShip.nst): + + if(phi<=0): + side=1 + else: + side=-1 + secData=self.compImmersedSectionProperties(GivenShip.secCrv[ist],zd,phi,side) + + 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] + + # Py-NAT 2021 a : PRAVEEN KUMAR CH + # 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])) + # -------------------------------------------------------------------------------------- + + # -------------------------------------------------------------------------------------- + # Function to compute heeled and immersed half section properties + # Py-NAT 2021a : PRAVEEN KUMAR CH + # -------------------------------------------------------------------------------------- + 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) + + 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] + # Py-NAT 2021a : PRAVEEN KUMAR CH + #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 Read the Data from File + # Py-NAT 2021a : PRAVEEN KUMAR CH + # -------------------------------------------------------------------------------------- + def writeToFile(self,fs): + + print("\n\t\tPyNAT: Writing KN Curves Data to File:\t ...\\",fs.okn_path) + + if not os.path.exists(fs.las_dir): + os.makedirs(fs.las_dir) + + fp=open(fs.okn_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): + # Py-NAT 2021a : PRAVEEN KUMAR CH + 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") + # PyNAT 2021 a : PraveenKumarCH + 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.png"), + dpi=300, + format='png', + bbox_extra_artists=(lg,), + bbox_inches='tight') + plt.close('all') + # ---------------------------------------------------------------------------------------- \ No newline at end of file diff --git a/Source_Script/LAS_KN_1A.py b/Source_Script/LAS_KN_1A.py new file mode 100644 index 0000000..620d3f0 --- /dev/null +++ b/Source_Script/LAS_KN_1A.py @@ -0,0 +1,499 @@ +''' +************************************************************************************ + 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') + # ---------------------------------------------------------------------------------------- diff --git a/Source_Script/PyNAT.py b/Source_Script/PyNAT.py new file mode 100644 index 0000000..400b1b7 --- /dev/null +++ b/Source_Script/PyNAT.py @@ -0,0 +1,137 @@ +''' +**************************************************************** + File : main.py + + Purpose : The Base Scipt from where the execution begins + Class Objects are Created to perform analysis + + App : PyNAT - Python Based Naval Architecture Tools + Version : 2021a + Paradigm : Object Oriented Model + Author : Praveen Kumar Ch +***************************************************************** +''' + +from HullGeometry import HullGeometry +from HydroStatics import HydroStatics +from FileSystem import FileSystem +from LAS_KN_1A import LAS_KN_1A +from LAS_KN_0B import LAS_KN_0B +from LAS_GZ import LAS_GZ +from ENV_PROP import ENV_PROP +import matplotlib.pyplot as plt +import os +import msvcrt as m +import time + +#----------------------------------------------------------------------------- +# Creation of 'FileSystem' Object +# Object Contains The File Names and paths to be utilized by program +# See "FileSystem.py" +#----------------------------------------------------------------------------- +t0=time.time() +print("\n\n*********************************************************************") +print(" Py-NAT 2021a") +print("*********************************************************************") +fs=FileSystem() +#----------------------------------------------------------------------------- + +#----------------------------------------------------------------------------- +# Creation of 'ENV_PROP' Object +# Contains Environment properties like water density +#----------------------------------------------------------------------------- +print("\n\nThe Environemnt Properties are being read...\t...\t...\t...\t...") +ep=ENV_PROP(fs) +print("\nEnvironment Properties File Read") +#----------------------------------------------------------------------------- + +# Creation of 'HullGeometry' Object +# Object Contains Data defining the ship hull geometry +# object Contains Script to read and prepare Hull Geometry +# Contians Script to fit Splines +# See the Script file "HullGeometry.py" +#----------------------------------------------------------------------------- +print("\n\nThe Hull Geometry is being read and prepared...\t...\t...\t...\t...") +givenShip=HullGeometry(fs) +tend=time.time() +print("\nThe Hull Geometry Data Preperation Complete.") +print("\nCheck the Prepared Hull Geometry Data in Directory:\t...\\",fs.sec_dir) +print("\nTotal Time Elapsed for Hull Geometry Preperations :\t",tend-t0) +print("\nClose Body Plan Plot to Continue.") +plt.show() +#---------------------------------------------------------------------------------- +iLoops=0 +maxLoops=50 +while(iLoops Create and Store Curve Fitting Data for a Section Curve + -> Assumes symmetry: Stores and Computes for half sections only. + -> Section Data Stored + -> X location of the stations + -> Curve point coordinates (Y,Z) + -> Continuity Type at each point C + -> No. of Curvelets and Ends/Knuckle points + -> Spline Curve Fit Data of each Curvelet + -> Functions to compute and Return the Section Properties at a given draft + -> Follwing Section Properties computed at a given draft + -> Half Breadth + -> Area + -> Moments : Long., Trans. , Vertical. + -> Girth Lengths + + App : PyNAT - Python Based Naval Architecture Tools + Version : 2021.0 + Paradigm : Object Based Framework + Author : Praveen Kumar Ch +************************************************************************************* +''' +import numpy as np +import matplotlib.pyplot as plt +from scipy import interpolate +from scipy.integrate import quad +from scipy.misc import derivative + +class SectionCurves: + + # -------------------------------------------------------------------------------------- + # The Constructor for a Section + # PyNAT 2021 a : PraveenKumarCH + # -------------------------------------------------------------------------------------- + def __init__(self,X,Y,Z,C): + + # Set Data from Arguments + self.N=Y.size + self.X=X + self.Y=Y + self.Z=Z + self.C=C + self.chined=0 + self.fittype=1 + self.sf=0.25 + self.maxOrd=2 + + # Compute no of Curvelets and their end points + # temp=np.where(self.C==0) + # self.iCE=temp[0] + # self.nCP=self.iCE.size-1 + + # + self.nCrv=0 + self.sCrv=list() + self.eCrv=list() + self.sCrv.append(0) + flag=1 + for i in range(1,self.N): + if(flag==1 and self.C[i]==0): + self.nCrv+=1 + self.eCrv.append(i) + self.sCrv.append(i) + flag=1 + elif(flag ==1 and self.C[i]==-1): + self.nCrv+=1 + self.eCrv.append(i) + flag=0 + elif(flag==0): + self.sCrv.append(i) + flag=1 + + # Set Type to Chined Hull Section is all points are knuckles + if self.nCrv==(self.N-1): + self.chined=1 + + # Fit Splines Function for each curvelet + self.y_fz=list() + self.zy_fz=list() + self.yy_fz=list() + self.L_fz=list() + for i in range(self.nCrv): + s=self.sCrv[i] + e=self.eCrv[i] + + if(self.fittype==0): + # Interpolation Function Cubic Spline Fit + e=e+1 + self.y_fz.append(interpolate.interp1d(self.Z[s:e],self.Y[s:e],kind='cubic')) + self.zy_fz.append(interpolate.interp1d(self.Z[s:e],self.Z[s:e]*self.Y[s:e],kind='cubic')) + self.yy_fz.append(interpolate.interp1d(self.Z[s:e],self.Y[s:e]*self.Y[s:e]/2,kind='cubic')) + df=np.arange(e-s,dtype=np.float64) + df[0]=(self.y_fz[i](self.Z[s+1])-self.y_fz[i](self.Z[s]))/(self.Z[s+1]-self.Z[s]) + df[1:-1]=(self.y_fz[i](self.Z[s+2:e])-self.y_fz[i](self.Z[s:e-2]))/(self.Z[s+2:e]-self.Z[s:e-2]) + df[-1]=(self.y_fz[i](self.Z[e-1])-self.y_fz[i](self.Z[e-2]))/(self.Z[e-1]-self.Z[e-2]) + self.L_fz.append(interpolate.interp1d(self.Z[s:e],np.sqrt(1+df**2),kind='cubic')) + + if(self.fittype==1): + # Smoothing Function Cubic Spline Fit + if(e-s>self.maxOrd): + order=self.maxOrd + else: + order=e-s + e=e+1 + self.y_fz.append(interpolate.splrep(self.Z[s:e],self.Y[s:e],k=order,s=self.sf)) + self.zy_fz.append(interpolate.splrep(self.Z[s:e],self.Z[s:e]*self.Y[s:e],k=order,s=self.sf)) + self.yy_fz.append(interpolate.splrep(self.Z[s:e],self.Y[s:e]*self.Y[s:e]/2,k=order,s=self.sf)) + df=interpolate.splev(self.Z[s:e],self.y_fz[i],der=1) + self.L_fz.append(interpolate.splrep(self.Z[s:e],np.sqrt(1+df**2),k=order,s=self.sf)) + + # -------------------------------------------------------------------------------------- + + # -------------------------------------------------------------------------------------- + # Function to Compute half breadths at a given draft + # PyNAT 2021 a : PraveenKumarCH + # -------------------------------------------------------------------------------------- + def getHalfBreadth(self,draft): + + hBw=0.0 + for i in range(self.nCrv): + if(self.Z[self.sCrv[i]]<=draft and self.Z[self.eCrv[i]]>=draft): + + if(self.fittype==0): + tmp=self.y_fz[i](draft) + elif(self.fittype==1): + tmp=interpolate.splev(draft,self.y_fz[i],0) + + hBw=tmp + break + + return hBw + # -------------------------------------------------------------------------------------- + + # -------------------------------------------------------------------------------------- + # Function to Compute half breadths at a given draft + # PyNAT 2021 a : PraveenKumarCH + # -------------------------------------------------------------------------------------- + def getGirthLength(self,draft): + CL=0 + for i in range(self.nCrv): + if(self.Z[self.sCrv[i]]<=draft): + zs=self.Z[self.sCrv[i]] + ze=self.Z[self.eCrv[i]] + if draft=draft): + if(self.fittype==0): + tmp=self.y_fz[i](draft) + elif(self.fittype==1): + tmp=interpolate.splev(draft,self.y_fz[i],0) + hBw=tmp + break + + # Longitudinal Area moment of the Section till draft + LM=Ar*self.X + + # Check against occurance of NaN + # PyNAT 2021 a : Praveen KumarCH + if(Ar!=0): + YC=TM/Ar + ZC=VM/Ar + + # Create and return a list of computed section parameters till given draft + secData=list([hBw,Ar,LM,TM,VM,CL,YC,ZC]) + + return secData + # -------------------------------------------------------------------------------------- + + # -------------------------------------------------------------------------------------- + # Function to interpolate offset at given spacing of drafts + # PyNAT 2021 a : PraveenKumarCH + # -------------------------------------------------------------------------------------- + def getPoints(self,sp): + pts=list() + + zts=np.arange(self.Z[0],self.Z[-1],sp,dtype=np.float64) + yts=np.arange(zts.size,dtype=np.float64) + + for i in range(zts.size): + yts[i]=self.getHalfBreadth(zts[i]) + + pts.append(zts) + pts.append(yts) + + return pts + + # -------------------------------------------------------------------------------------- +