Files
create/src/Mod/Ship/shipHydrostatics/Tools.py

601 lines
21 KiB
Python

#***************************************************************************
#* *
#* Copyright (c) 2011, 2016 *
#* Jose Luis Cercos Pita <jlcercos@gmail.com> *
#* *
#* This program is free software; you can redistribute it and/or modify *
#* it under the terms of the GNU Lesser General Public License (LGPL) *
#* as published by the Free Software Foundation; either version 2 of *
#* the License, or (at your option) any later version. *
#* for detail see the LICENCE text file. *
#* *
#* This program is distributed in the hope that it will be useful, *
#* but WITHOUT ANY WARRANTY; without even the implied warranty of *
#* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
#* GNU Library General Public License for more details. *
#* *
#* You should have received a copy of the GNU Library General Public *
#* License along with this program; if not, write to the Free Software *
#* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 *
#* USA *
#* *
#***************************************************************************
import math
import random
from FreeCAD import Vector, Rotation, Matrix, Placement
import Part
import Units
import FreeCAD as App
import FreeCADGui as Gui
from PySide import QtGui, QtCore
import Instance
from shipUtils import Math
import shipUtils.Units as USys
DENS = 1.025 # [tons/m3], salt water
COMMON_BOOLEAN_ITERATIONS = 10
def placeShipShape(shape, draft, roll, trim):
"""Move the ship shape such that the free surface matches with the plane
z=0. The transformation will be applied on the input shape, so copy it
before calling this method if it should be preserved.
Position arguments:
shape -- Ship shape
draft -- Ship draft
roll -- Roll angle
trim -- Trim angle
Returned value:
The same transformed input shape. Just for debugging purposes, you can
discard it.
"""
# Roll the ship. In order to can deal with large roll angles, we are
# proceeding as follows:
# 1.- Applying the roll with respect the base line
# 2.- Recentering the ship in the y direction
# 3.- Readjusting the base line
shape.rotate(Vector(0.0, 0.0, 0.0), Vector(1.0, 0.0, 0.0), roll)
base_z = shape.BoundBox.ZMin
shape.translate(Vector(0.0, draft * math.sin(math.radians(roll)), -base_z))
# Trim the ship. In this case we only need to correct the x direction
shape.rotate(Vector(0.0, 0.0, 0.0), Vector(0.0, -1.0, 0.0), trim)
shape.translate(Vector(draft * math.sin(math.radians(trim)), 0.0, 0.0))
shape.translate(Vector(0.0, 0.0, -draft))
return shape
def getUnderwaterSide(shape):
"""Get the underwater shape, simply cropping the provided shape by the z=0
free surface plane.
Position arguments:
shape -- Solid shape to be cropped
Returned value:
Cropped shape. It is not modifying the input shape
"""
# Convert the shape into an active object
Part.show(shape)
orig = App.ActiveDocument.Objects[-1]
bbox = shape.BoundBox
xmin = bbox.XMin
xmax = bbox.XMax
ymin = bbox.YMin
ymax = bbox.YMax
zmin = bbox.ZMin
zmax = bbox.ZMax
# Create the "sea" box to intersect the ship
L = xmax - xmin
B = ymax - ymin
H = zmax - zmin
box = App.ActiveDocument.addObject("Part::Box","Box")
length_format = USys.getLengthFormat()
box.Placement = Placement(Vector(xmin - L, ymin - B, zmin - H),
Rotation(App.Vector(0,0,1),0))
box.Length = length_format.format(3.0 * L)
box.Width = length_format.format(3.0 * B)
box.Height = length_format.format(- zmin + H)
App.ActiveDocument.recompute()
common = App.activeDocument().addObject("Part::MultiCommon",
"UnderwaterSideHelper")
common.Shapes = [orig, box]
App.ActiveDocument.recompute()
if len(common.Shape.Solids) == 0:
# The common operation is failing, let's try moving a bit the free
# surface
msg = QtGui.QApplication.translate(
"ship_console",
"Boolean operation failed when trying to get the underwater side."
" The tool is retrying such operation slightly moving the free"
" surface position",
None,
QtGui.QApplication.UnicodeUTF8)
App.Console.PrintWarning(msg + '\n')
random_bounds = 0.01 * H
i = 0
while len(common.Shape.Solids) == 0 and i < COMMON_BOOLEAN_ITERATIONS:
i += 1
box.Height = length_format.format(
- zmin + H + random.uniform(-random_bounds, random_bounds))
App.ActiveDocument.recompute()
out = common.Shape
App.ActiveDocument.removeObject(common.Name)
App.ActiveDocument.removeObject(orig.Name)
App.ActiveDocument.removeObject(box.Name)
App.ActiveDocument.recompute()
return out
def areas(ship, n, draft=None,
roll=Units.parseQuantity("0 deg"),
trim=Units.parseQuantity("0 deg")):
""" Compute the ship transversal areas
Position arguments:
ship -- Ship object (see createShip)
n -- Number of points to compute
Keyword arguments:
draft -- Ship draft (Design ship draft by default)
roll -- Roll angle (0 degrees by default)
trim -- Trim angle (0 degrees by default)
Returned value:
List of sections, each section contains 2 values, the x longitudinal
coordinate, and the transversal area. If n < 2, an empty list will be
returned.
"""
if n < 2:
return []
if draft is None:
draft = ship.Draft
# Manipulate a copy of the ship shape
shape = getUnderwaterSide(placeShipShape(ship.Shape.copy(),
draft,
roll,
trim))
# Sections distance computation
bbox = shape.BoundBox
xmin = bbox.XMin
xmax = bbox.XMax
dx = (xmax - xmin) / (n - 1.0)
# Since we are computing the sections in the total length (not in the
# length between perpendiculars), we can grant that the starting and
# ending sections have null area
areas = [(Units.Quantity(xmin, Units.Length),
Units.Quantity(0.0, Units.Area))]
# And since we just need to compute areas we will create boxes with its
# front face at the desired transversal area position, computing the
# common solid part, dividing it by faces, and getting only the desired
# ones.
App.Console.PrintMessage("Computing transversal areas...\n")
App.Console.PrintMessage("Some Inventor representation errors can be"
" shown, please ignore them.\n")
for i in range(1, n - 1):
App.Console.PrintMessage("{0} / {1}\n".format(i, n - 2))
x = xmin + i * dx
try:
f = Part.Face(shape.slice(Vector(1,0,0), x))
except Part.OCCError:
areas.append((Units.Quantity(x, Units.Length),
Units.Quantity(0.0, Units.Area)))
continue
# It is a valid face, so we can add this area
areas.append((Units.Quantity(x, Units.Length),
Units.Quantity(f.Area, Units.Area)))
# Last area is equal to zero (due to the total length usage)
areas.append((Units.Quantity(xmax, Units.Length),
Units.Quantity(0.0, Units.Area)))
App.Console.PrintMessage("Done!\n")
return areas
def displacement(ship, draft, roll=0.0, trim=0.0, yaw=0.0):
""" Compute the ship displacement.
@param ship Ship instance.
@param draft Ship draft.
@param roll Ship roll angle.
@param trim Ship trim angle.
@param yaw Ship yaw angle. Ussually you don't want to use this
value.
@return [disp, B, Cb], \n
- disp = Ship displacement [ton].
- B = Bouyance center [m].
- Cb = Block coefficient.
@note Bouyance center will returned as a FreeCAD.Vector instance.
@note Returned Bouyance center is in the non modified ship coordinates
"""
# We will take a duplicate of ship shape in order to conviniently
# manipulate it
shape = ship.Shape.copy()
_draft = draft * Units.Metre.Value
# Roll the ship. In order to can deal with large roll angles, we are
# proceeding as follows:
# 1.- Applying the roll with respect the base line
# 2.- Recentering the ship in the y direction
# 3.- Readjusting the base line
shape.rotate(Vector(0.0, 0.0, 0.0), Vector(1.0, 0.0, 0.0), roll)
base_z = shape.BoundBox.ZMin
shape.translate(Vector(0.0,
_draft * math.sin(math.radians(roll)),
-base_z))
# Trim and yaw the ship. In this case we only need to correct the x
# direction
shape.rotate(Vector(0.0, 0.0, 0.0), Vector(0.0, -1.0, 0.0), trim)
shape.translate(Vector(_draft * math.sin(math.radians(trim)), 0.0, 0.0))
shape.rotate(Vector(0.0, 0.0, 0.0), Vector(0.0, 0.0, 1.0), yaw)
shape.translate(Vector(0.0, 0.0, -_draft))
Part.show(shape)
ship_shape = App.ActiveDocument.Objects[-1]
bbox = shape.BoundBox
xmin = bbox.XMin
xmax = bbox.XMax
ymin = bbox.YMin
ymax = bbox.YMax
zmin = bbox.ZMin
zmax = bbox.ZMax
# Create the "sea" box to intersect the ship
L = xmax - xmin
B = ymax - ymin
H = zmax - zmin
box = App.ActiveDocument.addObject("Part::Box","Box")
length_format = USys.getLengthFormat()
box.Placement = Placement(Vector(xmin - L, ymin - B, zmin - H),
Rotation(App.Vector(0,0,1),0))
box.Length = length_format.format(3.0 * L)
box.Width = length_format.format(3.0 * B)
box.Height = length_format.format(- zmin + H)
App.ActiveDocument.recompute()
common = App.activeDocument().addObject("Part::MultiCommon",
"DisplacementHelper")
common.Shapes = [ship_shape, box]
App.ActiveDocument.recompute()
if len(common.Shape.Solids) == 0:
# The common operation is failing, let's try moving a bit the free
# surface
msg = QtGui.QApplication.translate(
"ship_console",
"Boolean operation failed. The tool is retrying that slightly"
" moving the free surface position",
None,
QtGui.QApplication.UnicodeUTF8)
App.Console.PrintWarning(msg + '\n')
random_bounds = 0.01 * H
i = 0
while len(common.Shape.Solids) == 0 and i < COMMON_BOOLEAN_ITERATIONS:
i += 1
box.Height = length_format.format(
- zmin + H + random.uniform(-random_bounds, random_bounds))
App.ActiveDocument.recompute()
vol = 0.0
cog = Vector()
if len(common.Shape.Solids) > 0:
for solid in common.Shape.Solids:
vol += solid.Volume / Units.Metre.Value**3
sCoG = solid.CenterOfMass
cog.x = cog.x + sCoG.x * solid.Volume / Units.Metre.Value**4
cog.y = cog.y + sCoG.y * solid.Volume / Units.Metre.Value**4
cog.z = cog.z + sCoG.z * solid.Volume / Units.Metre.Value**4
cog.x = cog.x / vol
cog.y = cog.y / vol
cog.z = cog.z / vol
Vol = L * B * abs(bbox.ZMin) / Units.Metre.Value**3
App.ActiveDocument.removeObject(common.Name)
App.ActiveDocument.removeObject(ship_shape.Name)
App.ActiveDocument.removeObject(box.Name)
App.ActiveDocument.recompute()
# Undo the transformations
B = Part.Point(Vector(cog.x, cog.y, cog.z))
m = Matrix()
m.move(Vector(0.0, 0.0, draft))
m.rotateZ(-math.radians(yaw))
m.move(Vector(-draft * math.sin(math.radians(trim)), 0.0, 0.0))
m.rotateY(math.radians(trim))
m.move(Vector(0.0,
-draft * math.sin(math.radians(roll)),
base_z / Units.Metre.Value))
m.rotateX(-math.radians(roll))
B.transform(m)
# Return the computed data
return [DENS*vol, Vector(B.X, B.Y, B.Z), vol/Vol]
def wettedArea(shape, draft, trim):
""" Calculate wetted ship area.
@param shape Ship external faces instance.
@param draft Draft.
@param trim Trim in degrees.
@return Wetted ship area.
"""
area = 0.0
nObjects = 0
shape = shape.copy()
_draft = draft * Units.Metre.Value
shape.rotate(Vector(0.0, 0.0, 0.0), Vector(0.0, -1.0, 0.0), trim)
shape.translate(Vector(0.0, 0.0, -_draft))
bbox = shape.BoundBox
xmin = bbox.XMin
xmax = bbox.XMax
ymin = bbox.YMin
ymax = bbox.YMax
# Create the "sea" box
L = xmax - xmin
B = ymax - ymin
p = Vector(xmin - L, ymin - B, bbox.ZMin - 1.0)
try:
box = Part.makeBox(3.0 * L, 3.0 * B, - bbox.ZMin + 1.0, p)
except Part.OCCError:
return 0.0
for f in shape.Faces:
try:
common = box.common(f)
except Part.OCCError:
continue
area = area + common.Area
return area / Units.Metre.Value**2
def moment(ship, draft, trim, disp, xcb):
""" Calculate triming 1cm ship moment.
@param ship Selected ship instance
@param draft Draft.
@param trim Trim in degrees.
@param disp Displacement at selected draft and trim.
@param xcb Bouyance center at selected draft and trim.
@return Moment to trim ship 1cm (ton m).
@note Moment is positive when produce positive trim.
"""
factor = 10.0
angle = factor * math.degrees(math.atan2(
0.01,
0.5 * ship.Length.getValueAs('m').Value))
newTrim = trim + angle
data = displacement(ship, draft, 0.0, newTrim, 0.0)
mom0 = -disp * xcb
mom1 = -data[0] * data[1].x
return (mom1 - mom0) / factor
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, and floating coefficient.
"""
area = 0.0
cf = 0.0
maxX = 0.0
minX = 0.0
maxY = 0.0
minY = 0.0
shape = ship.Shape.copy()
_draft = draft * Units.Metre.Value
shape.rotate(Vector(0.0, 0.0, 0.0), Vector(0.0, -1.0, 0.0), trim)
shape.translate(Vector(0.0, 0.0, -_draft))
bbox = shape.BoundBox
xmin = bbox.XMin
xmax = bbox.XMax
ymin = bbox.YMin
ymax = bbox.YMax
# Create the "sea" box
L = xmax - xmin
B = ymax - ymin
p = Vector(xmin - L, ymin - B, bbox.ZMin - 1.0)
try:
box = Part.makeBox(3.0 * L, 3.0 * B, - bbox.ZMin + 1.0, p)
except Part.OCCError:
return [area, cf]
maxX = bbox.XMin / Units.Metre.Value
minX = bbox.XMax / Units.Metre.Value
maxY = bbox.YMin / Units.Metre.Value
minY = bbox.YMax / Units.Metre.Value
for s in shape.Solids:
try:
common = box.common(s)
except Part.OCCError:
continue
if common.Volume == 0.0:
continue
# Recompute the object adding it to the scene. OpenCASCADE must be
# performing an internal tesellation doing that
try:
Part.show(common)
except (TypeError,Part.OCCError):
continue
# Divide the solid by faces and filter the well placed ones
faces = common.Faces
for f in faces:
faceBounds = f.BoundBox
# Orientation filter
if faceBounds.ZMax - faceBounds.ZMin > 0.00001:
continue
# Position filter
if abs(faceBounds.ZMax) > 0.00001:
continue
area = area + f.Area / Units.Metre.Value**2
maxX = max(maxX, faceBounds.XMax / Units.Metre.Value)
minX = min(minX, faceBounds.XMin / Units.Metre.Value)
maxY = max(maxY, faceBounds.YMax / Units.Metre.Value)
minY = min(minY, faceBounds.YMin / Units.Metre.Value)
App.ActiveDocument.removeObject(App.ActiveDocument.Objects[-1].Name)
dx = maxX - minX
dy = maxY - minY
if dx*dy > 0.0:
cf = area / (dx * dy)
return [area, cf]
def BMT(ship, draft, trim=0.0):
""" Calculate ship Bouyance center transversal distance.
@param ship Ship instance.
@param draft Ship draft.
@param trim Ship trim angle.
@return BM Bouyance to metacenter height [m].
"""
nRoll = 2
maxRoll = 7.0
B0 = displacement(ship, draft, 0.0, trim, 0.0)[1]
BM = 0.0
for i in range(nRoll):
roll = (maxRoll / nRoll)*(i + 1)
B1 = displacement(ship, draft, roll, trim, 0.0)[1]
# * M
# / \
# / \ BM ==|> BM = (BB/2) / sin(alpha/2)
# / \
# *-------*
# BB
BB = [B1.y - B0.y, B1.z - B0.z]
BB = math.sqrt(BB[0] * BB[0] + BB[1] * BB[1])
# nRoll is acting as the weight function
BM = BM + 0.5 * BB / math.sin(math.radians(0.5 * roll)) / nRoll
return BM
def mainFrameCoeff(ship, draft):
""" Calculate main frame coefficient.
@param ship Selected ship instance
@param draft Draft.
@return Main frame coefficient
"""
cm = 0.0
maxY = 0.0
minY = 0.0
shape = ship.Shape.copy()
shape.translate(Vector(0.0, 0.0, -draft * Units.Metre.Value))
x = 0.0
area = 0.0
bbox = shape.BoundBox
xmin = bbox.XMin
xmax = bbox.XMax
ymin = bbox.YMin
ymax = bbox.YMax
# Create the "sea" box
L = xmax - xmin
B = ymax - ymin
p = Vector(xmin - L, ymin - B, bbox.ZMin - 1.0)
try:
box = Part.makeBox(1.5 * L, 3.0 * B, - bbox.ZMin + 1.0, p)
except Part.OCCError:
return cm
maxY = bbox.YMin / Units.Metre.Value
minY = bbox.YMax / Units.Metre.Value
for s in shape.Solids:
try:
common = box.common(s)
except Part.OCCError:
continue
if common.Volume == 0.0:
continue
# Recompute the object adding it to the scene. OpenCASCADE must be
# performing an internal tesellation doing that
try:
Part.show(common)
except (TypeError,Part.OCCError):
continue
# Divide the solid by faces and filter the well placed ones
faces = common.Faces
for f in faces:
faceBounds = f.BoundBox
# Orientation filter
if faceBounds.XMax - faceBounds.XMin > 0.00001:
continue
# Position filter
if abs(faceBounds.XMax - x) > 0.00001:
continue
area = area + f.Area / Units.Metre.Value**2
maxY = max(maxY, faceBounds.YMax / Units.Metre.Value)
minY = min(minY, faceBounds.YMin / Units.Metre.Value)
App.ActiveDocument.removeObject(App.ActiveDocument.Objects[-1].Name)
dy = maxY - minY
if dy * draft > 0.0:
cm = area / (dy * draft)
return cm
class Point:
""" Hydrostatics point, that conatins: \n
draft Ship draft [m]. \n
trim Ship trim [deg]. \n
disp Ship displacement [ton]. \n
xcb Bouyance center X coordinate [m].
wet Wetted ship area [m2].
mom Triming 1cm ship moment [ton m].
farea Floating area [m2].
KBt Transversal KB height [m].
BMt Transversal BM height [m].
Cb Block coefficient.
Cf Floating coefficient.
Cm Main frame coefficient.
@note Moment is positive when produce positive trim.
"""
def __init__(self, ship, faces, draft, trim):
""" Use all hydrostatics tools to define a hydrostatics
point.
@param ship Selected ship instance
@param faces Ship external faces
@param draft Draft.
@param trim Trim in degrees.
"""
# Hydrostatics computation
dispData = displacement(ship, draft, 0.0, trim, 0.0)
if not faces:
wet = 0.0
else:
wet = wettedArea(faces, draft, trim)
mom = moment(ship, draft, trim, dispData[0], dispData[1].x)
farea = FloatingArea(ship, draft, trim)
bm = BMT(ship, draft, trim)
cm = mainFrameCoeff(ship, draft)
# Store final data
self.draft = draft
self.trim = trim
self.disp = dispData[0]
self.xcb = dispData[1].x
self.wet = wet
self.farea = farea[0]
self.mom = mom
self.KBt = dispData[1].z
self.BMt = bm
self.Cb = dispData[2]
self.Cf = farea[1]
self.Cm = cm