From 223e49af2f5ab501d69a4becee76bbe22a6bf9e4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jose=20Luis=20Cerc=C3=B3s=20pita?= Date: Fri, 9 Mar 2012 13:38:55 +0100 Subject: [PATCH] Added stability hydrostatics --- src/Mod/Ship/shipHydrostatics/Plot.py | 36 ++- src/Mod/Ship/shipHydrostatics/Tools.py | 355 ++++++++++++++++++++++++- 2 files changed, 387 insertions(+), 4 deletions(-) diff --git a/src/Mod/Ship/shipHydrostatics/Plot.py b/src/Mod/Ship/shipHydrostatics/Plot.py index 72e81a0a8..5e4f4d327 100644 --- a/src/Mod/Ship/shipHydrostatics/Plot.py +++ b/src/Mod/Ship/shipHydrostatics/Plot.py @@ -59,6 +59,7 @@ class Plot(object): if self.execute(): return ImageGui.open(self.path + 'volume.png') + ImageGui.open(self.path + 'stability.png') def createDirectory(self): """ Create needed folder to write data and scripts. @@ -97,13 +98,16 @@ class Plot(object): Output.write(" # 3: Wetted surface [m2]\n") Output.write(" # 4: 1cm triming ship moment [ton m]\n") Output.write(" # 5: Bouyance center x coordinate\n") + Output.write(" # 6: Floating area\n") + Output.write(" # 7: KBt\n") + Output.write(" # 8: BMt\n") Output.write(" #\n") Output.write(" #################################################################\n") # Print data for i in range(0,len(drafts)): draft = drafts[i] point = Tools.Point(ship,draft,trim) - string = "%f %f %f %f %f\n" % (point.disp, point.draft, point.wet, point.mom, point.xcb) + string = "%f %f %f %f %f %f %f %f\n" % (point.disp, point.draft, point.wet, point.mom, point.xcb, point.farea, point.KBt, point.BMt) Output.write(string) # Close file Output.close() @@ -160,13 +164,28 @@ class Plot(object): Output.write("set style 1 line linetype 1 linewidth 1 colour rgb (0):(0):(0)\n") Output.write("set style 2 line linetype 1 linewidth 1 colour rgb (1):(0):(0)\n") Output.write("set style 3 line linetype 1 linewidth 1 colour rgb (0):(0):(1)\n") - Output.write("set style 4 line linetype 1 linewidth 1 colour rgb (1):(0):(1)\n") + Output.write("set style 4 line linetype 1 linewidth 1 colour rgb (0.1):(0.5):(0.1)\n") # Write plot call Output.write("# Plot\n") - Output.write("plot '%s' using 2:1 title '$\\bigtriangleup$' axes x1y1 with lines style 1, \\\n" % (self.dataFile)) + Output.write("plot '%s' using 2:1 title 'Draft' axes x1y1 with lines style 1, \\\n" % (self.dataFile)) Output.write(" '' using 3:1 title 'Wetted area' axes x2y1 with lines style 2, \\\n") Output.write(" '' using 4:1 title '1cm trim moment' axes x3y1 with lines style 3, \\\n") Output.write(" '' using 5:1 title 'XCB' axes x4y1 with lines style 4\n") + # Prepare second plot + Output.write("set output '%s'\n" % (self.path + 'stability.eps')) + Output.write("# X axis\n") + Output.write("set x2label '\\textit{Floating area} / $\\mathrm{m}^2$'\n") + Output.write("set x2tic\n") + Output.write("set x3label '$KB_{T}$ / $\\mathrm{m}$'\n") + Output.write("set x3tic\n") + Output.write("set x4label '$BM_{T}$ / $\\mathrm{m}$'\n") + Output.write("set x4tic\n") + # Write plot call + Output.write("# Plot\n") + Output.write("plot '%s' using 2:1 title 'Draft' axes x1y1 with lines style 1, \\\n" % (self.dataFile)) + Output.write(" '' using 6:1 title 'Floating area' axes x2y1 with lines style 2, \\\n") + Output.write(" '' using 7:1 title '$KB_{T}$' axes x3y1 with lines style 3, \\\n") + Output.write(" '' using 8:1 title '$BM_{T}$' axes x4y1 with lines style 4\n") # Close file self.layoutFile = filename Output.close() @@ -176,6 +195,7 @@ class Plot(object): """ Calls pyxplot in order to plot an save an image. @return True if error happens. """ + # Plot filename = self.path + 'volume' comm = "pyxplot %s" % (self.layoutFile) if os.system(comm): @@ -184,6 +204,16 @@ class Plot(object): msg = Translator.translate("Plot will not generated\n") FreeCAD.Console.PrintError(msg) return True + # Convert volume + comm = "gs -r300 -dEPSCrop -dTextAlphaBits=4 -sDEVICE=png16m -sOutputFile=%s.png -dBATCH -dNOPAUSE %s.eps" % (filename,filename) + if os.system(comm): + msg = Translator.translate("Can't execute ghostscript. Maybe is not installed?\n") + FreeCAD.Console.PrintError(msg) + msg = Translator.translate("Generated image will not converted to png\n") + FreeCAD.Console.PrintError(msg) + return True + # Convert stability + filename = self.path + 'stability' comm = "gs -r300 -dEPSCrop -dTextAlphaBits=4 -sDEVICE=png16m -sOutputFile=%s.png -dBATCH -dNOPAUSE %s.eps" % (filename,filename) if os.system(comm): msg = Translator.translate("Can't execute ghostscript. Maybe is not installed?\n") diff --git a/src/Mod/Ship/shipHydrostatics/Tools.py b/src/Mod/Ship/shipHydrostatics/Tools.py index 25737e7de..b36ed8d2e 100644 --- a/src/Mod/Ship/shipHydrostatics/Tools.py +++ b/src/Mod/Ship/shipHydrostatics/Tools.py @@ -355,6 +355,350 @@ def Moment(ship, draft, trim, disp, xcb): mom1 = -data[1]*data[2] return mom1 - mom0 +def FloatingArea(ship, draft, trim): + """ Calculate ship floating area. + @param ship Selected ship instance + @param draft Draft. + @param trim Trim in degrees. + @return Ship floating area. + """ + angle = math.radians(trim) + sections = Instance.sections(ship) + xCoord = ship.xSection[:] + lines = [] + area = 0.0 + if not sections: + return 0.0 + for i in range(0, len(sections)): + # Get the section + section = sections[i] + if len(section) < 2: # Empty section + lines.append(0.0) + continue + # Get the position of the section + x = xCoord[i] + # Get the maximum Z value + Z = draft - x*math.tan(angle) + # Format section + section = convertSection(section, x, Z) + if not section: + lines.append(0.0) + continue + # Get floating line length + line = 0.0 + flag = True # Even lines compute for floating areas, odd no + j = len(section)-1 + k = len(section[j])-1 + while k>0: + k = k-1 + if flag: + y0 = abs(section[j][k-1].y) + y1 = abs(section[j][k].y) + line = line + (y1 - y0) + flag = not flag + if flag: # Central body computation lefts + y = abs(section[j][0].y) + line = line + y + lines.append(2.0*line) # 2x because only half ship is represented + # Add area if proceed + if i > 0: + dx = xCoord[i] - xCoord[i-1] + x = 0.5*(xCoord[i] + xCoord[i-1]) + line = 0.5*(lines[i] + lines[i-1]) + area = area + line*dx + return area + +def KBT(ship, draft, trim, roll=0.0): + """ Calculate ship Keel to Bouyance center transversal distance. + @param ship Selected ship instance + @param draft Draft. + @param trim Trim in degrees. + @param roll Roll angle in degrees. + @return [KBTx, KBTy]: \n + KBTy : TRansversal KB y coordinate \n + KBTz : TRansversal KB z coordinate + """ + angle = math.radians(trim) + rAngle = math.radians(roll) + sections = Instance.sections(ship) + xCoord = ship.xSection[:] + areas = [] + vol = 0.0 + vMy = 0.0 + vMz = 0.0 + My = [] + Mz = [] + kb = [0.0,0.0] + if not sections: + return [[],0.0,0.0] + for i in range(0, len(sections)): + # ------------------------------------ + # Board + # ------------------------------------ + # Get the section + section = sections[i] + if len(section) < 2: # Empty section + areas.append(0.0) + continue + # Get the position of the section + x = xCoord[i] + # Get the maximum Z value + Z = draft - x*math.tan(angle) + # Format section + aux = convertSection(section, x, Z) + if not aux: + areas.append(0.0) + My.append(0.0) + Mz.append(0.0) + continue + # Correct by roll angle + pts = aux[len(aux)-1] + beam = pts[0].y + for i in range(1,len(pts)): + beam = max(beam, pts[i].y) + Z = Z - beam*math.tan(rAngle) + # Format section + section = convertSection(section, x, Z) + if not section: + areas.append(0.0) + My.append(0.0) + Mz.append(0.0) + continue + # Integrate area + area = 0.0 + momy = 0.0 + momz = 0.0 + for j in range(0, len(section)-1): + for k in range(0, min(len(section[j])-1, len(section[j+1])-1)): + # y11,z11 ------- y01,z01 + # | | + # | | + # | | + # y10,z10 ------- y00,z00 + y00 = abs(section[j][k].y) + z00 = section[j][k].z + y10 = abs(section[j][k+1].y) + z10 = section[j][k+1].z + y01 = abs(section[j+1][k].y) + z01 = section[j+1][k].z + y11 = abs(section[j+1][k+1].y) + z11 = section[j+1][k+1].z + dy = 0.5*((y00 - y10) + (y01 - y11)) + dz = 0.5*((z01 - z00) + (z11 - z10)) + y = 0.25*(y00 + y10 + y01 + y11) + z = 0.25*(z01 + z00 + z11 + z10) + area = area + dy*dz + momy = momy + y*dy*dz + momz = momz + z*dy*dz + if(len(section[j]) < len(section[j+1])): + # y01,z01 ------- y11,z11 + # | __/ + # | __/ + # | / + # y00,z00 + k = len(section[j])-1 + y00 = abs(section[j][k].y) + z00 = section[j][k].z + y01 = abs(section[j+1][k].y) + z01 = section[j+1][k].z + y11 = abs(section[j+1][k+1].y) + z11 = section[j+1][k+1].z + dy = y01 - y11 + dz = z01 - z00 + y = 0.33*(y00 + y01 + y11) + z = 0.33*(z01 + z00 + z11) + area = area + 0.5*dy*dz + momy = momy + y*0.5*dy*dz + momz = momz + z*0.5*dy*dz + elif(len(section[j]) > len(section[j+1])): + # y01,z01 + # | \__ + # | \__ + # | \ + # y00,z00 ------- y10,z10 + k = len(section[j+1])-1 + y00 = abs(section[j][k].y) + z00 = section[j][k].z + y10 = abs(section[j][k+1].y) + z10 = section[j][k+1].z + y01 = abs(section[j+1][k].y) + z01 = section[j+1][k].z + dy = y00 - y10 + dz = z01 - z00 + y = 0.33*(y00 + y01 + y10) + z = 0.33*(z01 + z00 + z10) + area = area + 0.5*dy*dz + momy = momy + y*0.5*dy*dz + momz = momz + z*0.5*dy*dz + elif(len(section[j]) == 1): + # y1,z1 ------- + # | + # | + # | + # y0,z0 ------- + k = 0 + y0 = abs(section[j][k].y) + z0 = section[j][k].z + y1 = abs(section[j+1][k].y) + z1 = section[j+1][k].z + dy = 0.5 * (y0 + y1) + dz = z1 - z0 + y = 0.25*(y1 + y0) + z = 0.25*(z1 + z0) + area = area + dy*dz + momy = momy + y*dy*dz + momz = momz + z*dy*dz + # ------------------------------------ + # StarBoard + # ------------------------------------ + # Get the section + section = sections[i] + # Get the maximum Z value + Z = draft - x*math.tan(angle) + # Format section + aux = convertSection(section, x, Z) + if not aux: + areas.append(0.0) + My.append(0.0) + Mz.append(0.0) + continue + # Correct by roll angle + pts = aux[len(aux)-1] + beam = pts[0].y + for i in range(1,len(pts)): + beam = max(beam, pts[i].y) + Z = Z + beam*math.tan(rAngle) + # Format section + section = convertSection(section, x, Z) + if not section: + areas.append(0.0) + My.append(0.0) + Mz.append(0.0) + continue + # Integrate area + for j in range(0, len(section)-1): + for k in range(0, min(len(section[j])-1, len(section[j+1])-1)): + # y11,z11 ------- y01,z01 + # | | + # | | + # | | + # y10,z10 ------- y00,z00 + y00 = abs(section[j][k].y) + z00 = section[j][k].z + y10 = abs(section[j][k+1].y) + z10 = section[j][k+1].z + y01 = abs(section[j+1][k].y) + z01 = section[j+1][k].z + y11 = abs(section[j+1][k+1].y) + z11 = section[j+1][k+1].z + dy = 0.5*((y00 - y10) + (y01 - y11)) + dz = 0.5*((z01 - z00) + (z11 - z10)) + y = 0.25*(y00 + y10 + y01 + y11) + z = 0.25*(z01 + z00 + z11 + z10) + area = area + dy*dz + momy = momy - y*dy*dz + momz = momz + z*dy*dz + if(len(section[j]) < len(section[j+1])): + # y01,z01 ------- y11,z11 + # | __/ + # | __/ + # | / + # y00,z00 + k = len(section[j])-1 + y00 = abs(section[j][k].y) + z00 = section[j][k].z + y01 = abs(section[j+1][k].y) + z01 = section[j+1][k].z + y11 = abs(section[j+1][k+1].y) + z11 = section[j+1][k+1].z + dy = y01 - y11 + dz = z01 - z00 + y = 0.33*(y00 + y01 + y11) + z = 0.33*(z01 + z00 + z11) + area = area + 0.5*dy*dz + momy = momy - y*0.5*dy*dz + momz = momz + z*0.5*dy*dz + elif(len(section[j]) > len(section[j+1])): + # y01,z01 + # | \__ + # | \__ + # | \ + # y00,z00 ------- y10,z10 + k = len(section[j+1])-1 + y00 = abs(section[j][k].y) + z00 = section[j][k].z + y10 = abs(section[j][k+1].y) + z10 = section[j][k+1].z + y01 = abs(section[j+1][k].y) + z01 = section[j+1][k].z + dy = y00 - y10 + dz = z01 - z00 + y = 0.33*(y00 + y01 + y10) + z = 0.33*(z01 + z00 + z10) + area = area + 0.5*dy*dz + momy = momy - y*0.5*dy*dz + momz = momz + z*0.5*dy*dz + elif(len(section[j]) == 1): + # y1,z1 ------- + # | + # | + # | + # y0,z0 ------- + k = 0 + y0 = abs(section[j][k].y) + z0 = section[j][k].z + y1 = abs(section[j+1][k].y) + z1 = section[j+1][k].z + dy = 0.5 * (y0 + y1) + dz = z1 - z0 + y = 0.25*(y1 + y0) + z = 0.25*(z1 + z0) + area = area + dy*dz + momy = momy - y*dy*dz + momz = momz + z*dy*dz + areas.append(area) + My.append(momy) + Mz.append(momz) + # Add volume & moment if proceed + if i > 0: + dx = xCoord[i] - xCoord[i-1] + x = 0.5*(xCoord[i] + xCoord[i-1]) + area = 0.5*(areas[i] + areas[i-1]) + momy = 0.5*(My[i] + My[i-1]) + momz = 0.5*(Mz[i] + Mz[i-1]) + vol = vol + area*dx + vMy = vMy + momy*dx + vMz = vMz + momz*dx + # Compute KBT + kb[0] = vMy / vol + kb[1] = vMz / vol + return kb + +def BMT(ship, draft, trim): + """ Calculate ship Bouyance center transversal distance. + @param ship Selected ship instance + @param draft Draft. + @param trim Trim in degrees. + @return BM Bouyance to metacenter height. + """ + nRoll = 2 + maxRoll = 7.0 + kb0 = KBT(ship,draft,trim) + BM = 0.0 + for i in range(0,nRoll): + roll = (maxRoll/nRoll)*(i+1) + kb = KBT(ship,draft,trim,roll) + # * M + # / \ + # / \ BM ==|> BM = (BB/2) / tan(alpha/2) + # / \ + # *-------* + # BB + BB = [kb[0] - kb0[0], kb[1] - kb0[1]] + BB = math.sqrt(BB[0]*BB[0] + BB[1]*BB[1]) + BM = BM + 0.5*BB/math.tan(math.radians(0.5*roll)) / nRoll # nRoll is the weight function + return BM + class Point: """ Hydrostatics point, that conatins: \n draft Ship draft [m]. \n @@ -362,7 +706,10 @@ class Point: disp Ship displacement [ton]. \n xcb Bouyance center X coordinate [m]. wet Wetted ship area [m2]. - mom triming 1cm ship moment [ton m]. + mom Triming 1cm ship moment [ton m]. + farea Floating area [m2]. + KBt Transversal KB height [m]. + BMt Transversal BM height [m]. @note Moment is positive when produce positive trim. """ def __init__(self, ship, draft, trim): @@ -376,10 +723,16 @@ class Point: areasData = Displacement(ship,draft,trim) wettedArea = WettedArea(ship,draft,trim) moment = Moment(ship,draft,trim,areasData[1],areasData[2]) + farea = FloatingArea(ship,draft,trim) + kb = KBT(ship,draft,trim) + bm = BMT(ship,draft,trim) # Store final data self.draft = draft self.trim = trim self.disp = areasData[1] self.xcb = areasData[2] self.wet = wettedArea + self.farea = farea self.mom = moment + self.KBt = kb[1] + self.BMt = bm