Calculate sailing ship drag, or anyway any drag of any shape

[WiredCat]

Hello there, I'm trying to write a sailing model simulator (for ships with sails) i got buoyancy force and gravity force calculated so far.
Before i go any further i need to tell you how i am calculating anything, and what is my ship exactly.
First of all i am making that on a phone, it looks like this:

Now let's consider that this ship on the picture is made of triangles, (set of points, that every 3 next points create closed shape (solid), its convex polyhedra), the most simple 3-dimensional shape is tetrahedron, which is made out of 4 triangles that 'close' the shape - this is important for you to know.

Now to calculate buoyancy force and drag acting on a ship (but only for submerged part) i slice that ship shape (set of triangles) like this: (green thing is my set of triangles on which i calculate buoyancy and drag) (note that i return the shape that is below waterline, i do not close this shape (this means if my ship is a simple tetrahedron and its covered in half with water, i only cut triangles (sides of that tetrahedron) and do not close (create new triangle) at the top of the new shape, like in this picture:

but that's irrevelant.Final shape can look like this:

I'll tell you how i calculate the buoyuancy but you can skip this part since its calculated properly.
For any number of triangles(in that green shape), create a tetrahedron that base is each triangle and apex is the center of mass of the ship, then calculate the volume of each such tetrahedron, then go again through all triangles and do the same but now we perform a cut on tetrahedrons, (we cut by water plane) and we substract each triangle volume by this cut tetrahedron to get total submerged area), after this i calculate geometric center of all submerged tetrahedrons and this gives me the point where i apply buoyancy force.

Forgot to mention that i also calculate area of each submerged triangle.Before i go any further i would like you to see what is the actual shape of my ship (its made out of 24 triangles, i later cut the ship model into 4 sections from front to back, and return with 120 triangles)
any way hull at bottom ain't flat it looks like a triangle, and at front there isn't any curve but flat 2 triangles that are angled:

Back:

Side: (front curve - the one at bottom)

Front:
You will see a quad (2 tris as a flat angled plate)

and from other perspective:

I need to use really small amount of triangles since its realtime simulation on mobile phones, and i need to run it as fast as possible (not mentioning i have there some time consuming things and somewhat simplified physically based rendering)

And i forgot to straight forward what is divided into sections and one thing: i only calculate ships hull not masts etc or this stick in front of it)
sections are made like that:

from 24 triangles to 120 :P (this is important since when the model is rotating drag applied on center of triangle can somehow correspond to real one since in low polygon model(24 tris) 2 triangles defined one side (where cannons are) so actually when such tris are cut horizontally drag will be applied almost at center of the ships hull) with divided tris i can apply drag force to non center of mass point 'a bit more preceisely'
I may not described that properly, sorry for that/Now... given velocity vector 'vel' [m/s] and frame time dt (consider 0.033 )
i calculate drag like this:

for all submerged triangles
{
get geometric center of each triangle, surface and volume of shape that is made out of tetrahedron which apex is on ships center of mass and its cut by waterplane...

get triangle local velocity = ship_linear_vel + ship_angular_vel x vector_form_tri_geom_center_to_ship_center_of_mass;

now i calculate form drag percentage coefficient (to know how much my triangle faces the velocity vector)
form_percentage_coeff = absolute_value( arcus_cosinus( dot_product(triangle_normal, vel_normalized)) / pi ); - reason i return absolute value is that even when triangle is not facing the same direction like vel, the suction force occurs (i managed to calculate that its 0.6666666 of element form drag force).999.1026 - water density
element_form_drag = (-local_velocity_normalized) * (form_percentage_coeff * (999.1026 * 0.5 * local_vel_squared) * triangle_surface_area);now for skin friction drag:
skin_percentage_coeff = absolute_value(1.0 - form_percentage_coeff);
then i do some magic that calculates the choord length for that triangle (projected local velocity vector, onto triangle surface, that comes through triangle geometric center)

0.894 - water dynamic viscocsity..
form this i calculate Reynolds number, and then a coefficient
CFRn = 1.328 / square_root(Rn);

element_skin_drag = (-normalized_local_vel) * (CFRn * ((999.1026 * local_vel_squared) / 2.0 ) * triangle_surface_area * skin_percentage_coeff);
drag_of_element = element_skin_drag + element_form_drag;

total_drag = total_drag + drag_of_element; //for linear motion
}i also sum all torques that are made up by all of tris + torque made by buoyancy force.
(this is another part of the post - at the end i will ask my questions)
I am not sure if my angular acceleration which is rad/s is correct. and if i am using my moment of inertia tensor properly.

after summing all torques i do something like this:
//angular acceleration.
AngAcc.x = summed_torque.x / Ixx;
AngAcc.y = summed_torque.y / Iyy;
AngAcc.z = summed_torque.z / Izz;

now this is a bit weird for me, because x component describes rotation around x axis, and Ixx defines rotation around which axis uh? i think it defines rotation around Z one.

Anyway next question
after summing doing each iteration step: ship_angular_vel = ship_angular_vel + AngAcc*dt;

is that correct to get each triangle local vel by:
triangle local velocity = ship_linear_vel + ship_angular_vel x vector_form_tri_geom_center_to_ship_center_of_mass;

or maybe my units are bad? in example i add meters per second and radians per second cross some_vector.My ship is unstable when i apply thrust at center of mass (i mean after some time its front starts bobbing like mad (going up and going down in time the amplitude of bobbing increases to finally reach so ridicilous number that processor can't handle it) another thing is when ship is trying to make a turn. it turns but after some time it flips over and acts really weird (i could try to make a video of that if requested)So my questions are:
What i am doing wrong? - i know i am trying to simulate simplified shape (and try to use 'corners' instead of soft curves'
How is form drag and skin friction drag calculated for any surface (or for any shape that is made up out of tiny surfaces (just don't tell me i need to convert all these cut triangles into something else :X)
Is rotation, and angular velcoity correctly calculated.

How to improove the simulation so i can eaisly make turns with ship without causing it to flip over.here is code: (theres also a part with rudder deflection (for turning)

Spoiler: code

void ProcessOceanFrame(TWaterWaves * ocean, float dt) //rotate model due to water
{
    if (dt < 0.0010) return;

//    ALOG("FRAME DT: "+FloatToStr(dt));
    result_force             = vec3(0.0, 0.0, 0.0);
    vec3 thrust_force         = vec3(0.0, 0.0, 1.0) * 1000.0;vec3 elements_torque;
/*
for each element
calculate local velocity related as EAngVel
Calculate all forces acting on body at this frame (add linear force) and add torque
----

after that rotate body
move body
*/
/*
* Buoyancy effect
* ******************************************************************************************
*/

vec3 cog2cob_vec     = vectorAB(pos + ROTATION_MAT * CENTER_OF_GRAVITY, CENTER_OF_BUOYANCY);
vec3 buoyancy_force = vec3(0.0, 1.0, 0.0) * (GForce * SUBMERGED_VOLUME * 999.1026);vec3 ETorque = cog2cob_vec    *  buoyancy_force;

elements_torque = ETorque;
/*
* Eof buoyancy effect
* ******************************************************************************************
*/

vec3 gravity_force = vec3(0.0, -1.0, 0.0) * (mass*9.81);

vec3 drag_force;vec3 form_drag;
vec3 skin_drag;

for (int i=0; i < submerged_tri_count; i++)
if (ACTUAL_FRAME[i].surface > 0.1) //do not coun't anything that is less than 10 cm^2
{
//ALOG("SKIN FRICTION PRESENT");
    vec3 element_form_drag = vec3(0.0, 0.0, 0.0); //somehow i don't trust that compiler due to optimization will zero it
    vec3 element_skin_drag = vec3(0.0, 0.0, 0.0);
    vec3 cog2pC     = vectorAB(pos + ROTATION_MAT * CENTER_OF_GRAVITY, ACTUAL_FRAME[i].pC);
    vec3 local_vel = vel + AngVel * cog2pC;
    vec3 vn = Normalize(local_vel);
    float vv = VectorLength(local_vel);
    float v2 = vv;
    vv = vv * vv;
float form_percentage_coeff = absnf( acos(dot(ACTUAL_FRAME[i].normal, vn)) / float(pi) );
if (dot(ACTUAL_FRAME[i].normal, vn) <= 0.0) form_percentage_coeff = form_percentage_coeff*0.66; //suction force

element_form_drag = (-vn) * (form_percentage_coeff * (999.1026 * 0.5 * vv) * ACTUAL_FRAME[i].surface); //form drag

float skin_percentage_coeff = absnf(1.0 - form_percentage_coeff); //due to float inaccuracy
vec3 surf_n = Normal(ACTUAL_FRAME[i].V[0], ACTUAL_FRAME[i].V[1], ACTUAL_FRAME[i].V[2], true);
float surf_choordlen = getTriangleChoordLength(
-local_vel, surf_n, ACTUAL_FRAME[i].pC,
ACTUAL_FRAME[i].V[0], ACTUAL_FRAME[i].V[1], ACTUAL_FRAME[i].V[2], true, ACTUAL_FRAME[i].pC);

if (surf_choordlen > 20.0) surf_choordlen = 3.0;

//calc reynolds number
float Rn = ReynoldsNum(VectorLength(local_vel), surf_choordlen, 0.894);
bool ah = false;
if (absnf(Rn) <= 0.001 ) {Rn = 1.0; ah= true;}

float CFRn = 1.328 / (sqrt(Rn));
if (ah) CFRn = 0.0;
//float CFRn = 0.075 / (lg * lg); //1.328 / (sqrt(Rn));//

vec3 Vfi = Normalize(ProjectVectorOnPlane(ACTUAL_FRAME[i].normal, -local_vel));
element_skin_drag = (-vn) * (CFRn * ((999.1026 * vv) / 2.0 ) * ACTUAL_FRAME[i].surface * skin_percentage_coeff);
skin_drag = skin_drag + element_skin_drag;
form_drag = form_drag + element_form_drag;

cog2cob_vec = vectorAB(pos + ROTATION_MAT * CENTER_OF_GRAVITY, ACTUAL_FRAME[i].pC);

ETorque = cog2cob_vec    *  ( element_skin_drag + element_form_drag );

elements_torque = elements_torque + ETorque;

}

drag_force = form_drag + skin_drag;
float yv = 30.0;
vec3 vf = YPRangle.GetFrontVector(); vf.y = 0.0;
vec3 thrust_vec = vf * (1.2 * 1.250 * 0.5 * yv*yv * 800.0);
if (pos.z > 520.0)
    if (!rudder_deflection)
{
TurnLeft();
}
// Rudder deflection
if (rudder_deflection)
{

    float deflection_angle;
    if (rudder_left) deflection_angle = 60.0;
    else deflection_angle = -60.0;

    Matrix44<float> rudmat;
    rudmat.RotateY(deflection_angle*imopi);    vec3 cog2cor_vec = vectorAB(pos + ROTATION_MAT * CENTER_OF_GRAVITY, pos + ROTATION_MAT * rudder_pos);    vec3 local_vel = vel + AngVel * cog2cor_vec;
    vec3 vn = Normalize(local_vel);
    float vv = VectorLength(local_vel);
    vv = vv * vv;    rudder_force = -vn * 1.2 *  ( ( 998.0 * vv) / 2.0 ) * 2.33 * 100.0;
    vec3 rETorque = cog2cor_vec    *  rudder_force;

    elements_torque = elements_torque + rETorque;

}

result_force = gravity_force + buoyancy_force + drag_force;// + thrust_vec + rudder_force;

mat4 inv_rot = ROTATION_MAT;
inv_rot.Inverse();

vec3 local_torque = inv_rot * elements_torque;

vec3 EAngAcc;
EAngAcc.x = local_torque.x / Ixx;//Ixx;
EAngAcc.y = local_torque.y / Iyy;
EAngAcc.z = local_torque.z / Izz;

EAngVel = EAngAcc*dt;
//ALOG("EAngAcc: " +POINT_TO_TEXT(EAngAcc));AngVel = AngVel + EAngVel;

vec3 vLocalAngularMoment = (AngVel*dt)*RAD_TO_DEG;
YPRangle.pitch(cos(vLocalAngularMoment.x*imopi), -sin(vLocalAngularMoment.x*imopi));
YPRangle.DoRotation();

YPRangle.yaw(cos(vLocalAngularMoment.y*imopi), sin(vLocalAngularMoment.y*imopi));
YPRangle.DoRotation();

YPRangle.roll(cos(vLocalAngularMoment.z*imopi), -sin(vLocalAngularMoment.z*imopi));
YPRangle.DoRotation();ROTATION_MAT.LoadGLMatrix(YPRangle.AIR_MATRIX);

vec3 accel = result_force / mass;

vel = vel + accel*dt;
pos = pos + vel*dt;

}void StepSim(TWaterWaves * ocean, float dt)
{

SliceHullByWater(ocean);
SliceAndCalculateSubmergedHullArea();
CalculateBuoyancy();
ShootOutSubmergedShape();
ProcessOceanFrame(ocean, dt);
}

Please help the kitten ;]

 

[A.T.]

My ship is unstable when i apply thrust at center of mass (i mean after some time its front starts bobbing like mad

You need more damping, that dissipates the energy.

 

[anorlunda]

Drag is irrelevant to stability. Did you put ballast in the ship?

To be stable, the center of mass must lie below the center of buoyancy at all times.

 

[SteamKing]

Drag is irrelevant to stability. Did you put ballast in the ship?

To be stable, the center of mass must lie below the center of buoyancy at all times.

That may be true mathematically (at least it must be true for fully submerged vessels to remain stable), but it is rarely true for actual vessels which are floating.

For stability at small angles, the center of gravity of the vessel must be below a point called the metacenter, whose location is determined in part by the location of the center of buoyancy of the displaced volume and by the shape of the waterplane.

https://en.wikipedia.org/wiki/Metacentric_height

 

[anorlunda]

OK, metacener is a more precise term. However, it's dead wrong to say that the COB is rarely above the COG.

Did you put ballast in the vessel? Is the metacenter above the COG?

Drag is irrelevant for an object that is not moving, or moving very slowly.

 

[SteamKing]

OK, metacener is a more precise term. However, it's dead wrong to say that the COB is rarely above the COG.

In my experience as a naval architect with over 35 years of experience doing floating vessel stability, I can say that the COB of all the vessels I have worked on has been located safely below the waterline, where it should be. I cannot say, with similar confidence, that the COG has been located below the COB, because it has not.

For certain vessels, namely deck barges loaded with cargo, the COG of the loaded vessel is often located above the main deck, well above the location of the COB, yet the barge is perfectly stable. Even in the unloaded condition, the COG of the barge hull is typically located between 50% and 60% of the depth of the hull above the bottom, while the COB is located below this point due to the small draft typical of these vessels in the light condition.

For many other vessels, the operating COG is located somewhere around the location of the main deck, in other words, above the waterline.

If you wanted to build a vessel whose COG remained below the COB in the operating condition, it would be a strange-looking beast, shallow and broad, having little or no superstructure, with lots of fixed ballast strapped to the keel.

You certainly can't say that a modern cruise ship has a COG below the COB:

​

 

[anorlunda]

Let's return to the OP. You said that the vessel you designed is unstable and you don't know why. We are trying to help via questions like ballast.

Porpoising tendencies are sometimes solved by redistributing ballast.

By the way, I double checked on COB relative to COG. I had it backwards. Apologies.

 

[SteamKing]

Let's return to the OP. You said that the vessel you designed is unstable and you don't know why. We are trying to help via questions like ballast.

Porpoising tendencies are sometimes solved by redistributing ballast.

By the way, I double checked on COB relative to COG. I had it backwards. Apologies.

It might be more advantageous at this point to find out where the OP has located the center of gravity for this vessel, before worrying about ballast and whatnot.
It's still not clear if the OP knows how to determine if a vessel is stable transversely or longitudinally by calculation.

The OP seems to be trying to design a 3-D motion vessel simulation by the seat of his pants, and I'm not sure he hasn't made some major errors in doing so. The motions of vessels in seaways is a difficult problem to solve, but there are codes which can predict the motion of vessels if given a sufficiently accurate description of the loading of the vessel.

Maneuverability of vessels still must be assessed using either turning trials of the actual vessel or a model test in a special turning basin; I'm aware of no software which allows for the prediction of the maneuverability characteristics of a vessel based solely on a geometric description of the hull.

 

[WiredCat]

The hull i described is uniform density and its mass is 120 000 kg so i do not use anything that is ballast. Ship is floating on water when its not moved, when in motion it acts like it shouldn't. and most of the times center of mass is located above center of buoyancy (in meaning that center of buoayncy is gemoetric center of submerged volume)

You need more damping, that dissipates the energy.

yeah it may have something with wave resistance, and suction force.Any additional ideas?

 

[SteamKing]

The hull i described is uniform density and its mass is 120 000 kg so i do not use anything that is ballast. Ship is floating on water when its not moved, when in motion it acts like it shouldn't. and most of the times center of mass is located above center of buoyancy (in meaning that center of buoayncy is gemoetric center of submerged volume)yeah it may have something with wave resistance, and suction force.Any additional ideas?

Yes, what is the height of the center of mass above the keel of the boat?
Have you calculated the location of the metacenter when the vessel is sitting still in the water?
What is the longitudinal location of the center of mass relative to the center of buoyancy?

In other words, have you checked the stability characteristics of this vessel when it is sitting still in still water?

 

[WiredCat]

Yes, what is the height of the center of mass above the keel of the boat?

- 1.66 of meter

Have you calculated the location of the metacenter when the vessel is sitting still in the water?

- i am not quite sure how to do that

What is the longitudinal location of the center of mass relative to the center of buoyancy?

-sorry my english is bad i do not understand that question

In other words, have you checked the stability characteristics of this vessel when it is sitting still in still water?

no, i can't make ship still without calculating proper drag.


I changed center of mass position to be lower, now on still water is situated somewhat 30 centimeters (0.3 m) above COB (i say somewhat because wrong drag calculations destabilize the model) i will post videos, as far as i get hands on another cellphone coz my one is too slow to run such simulation and additionally record it.

i still think i should change the shape of the bottom of the hull from flipped triangle to something more flat.
Anyways i need to calculate drag forces, then i will be able to fix center of mass position.

 

[SteamKing]

Why do you need to calculate drag forces to make the ship sit still? Why is the ship moving in the first place?

Calculating the resistance of an arbitrary hull is one of those problems which doesn't still have a complete analytic solution, which is why model testing of vessel designs is not completely obsolete.

 

[WiredCat]

Ah maybe you don't know that but applying only gravity and buoyancy forces, make the ship (lets consider we drop ship from 1 meter above waterline) bob forever. well maybe not forever i understand that gravity opposes buoyancy but still its not enough to just put it on waterplane and hope that after 2 hours it will stablizie, additionally buo force ain't acting at center of mass so it gives additionall angular momentum which destabilizes everything, in the evening (GMT+1) i will post videos of that and other.

it penetrates some water then its moved up then down then up then down up down forever. that's why i need drag. People are able somehow to calculate drag for some shapes. To ease calculations i would only consider laminar flow.

Calculating the resistance of an arbitrary hull is one of those problems which doesn't still have a complete analytic solution, which is why model testing of vessel designs is not completely obsolete.

Yeah anyway i would like to know how do they do that. the best would be paper that doesn't include calculus and any other mathematical symbols that i do not understand. Like for shape A divide it into x small areas S and do for each: calculate skin drag like this **, and form drag like this ** - i think i am asking too much now

 

[SteamKing]

If you iterate the location of the waterline of the vessel properly, it should reach an equilibrium between the gravity force and the buoyant force after a while, assuming the vessel is stable in the first place.

If you haven't studied the basic stability of a floating vessel, I don't hold out much optimism that you are going to realistically simulate either the resistance or motion of a vessel under arbitrary conditions. The fluid mechanics are just too complex, and the theoretical development of turbulent flow is too rudimentary to provide a complete solution without use of some empirical data, like that obtained from model testing.

 

[WiredCat]

If you haven't studied the basic stability of a floating vessel, I don't hold out much optimism..

Ah geez.
What do you mean by basic stability of a floating vessel, i have read some papers about making sailing model, one of it was talking about metacenter.

I am more likely treating that ship like a body that is moved by external forces.

like i said earlier ill post videos:
i would like to add, i draw there a dynamic water, but simulation is run that there's is only one water plane that means n 0,1,0 (that is pointing up) and d = 0 so i cut all tris by that plane, the reason i am showing dynamic water is that drawing simple quad is more complicated than it should be (when using phones) you can't just define 4 verts and call draw them you need to create buffers and else useless things, it would take too much time so i decided to draw dynamic water.and sorry for camera shake its hard too move camera (in game) and hold still phone camera in other hand
ok so first Byouancy + Gravity
Buo + Grav + Drag (main reason why i posted about how to calc drag)
(camera stabilized:P + one fancy effect do not be mad if you will see black screen for few seconds)
And wierdo thing at 1:10 it starts to bob (still you may see that ship penetrates water and its below it but water drawn is dynamic, i test by one static plane through all the sim time)
Buo + Grav + Drag + Thrust (applied at center of mass)


If one manages to see text (on the half of phone screen height) "y: number" its the y (up/down) distance from center of buo to center of mass

 

[SteamKing]

Ah geez.
What do you mean by basic stability of a floating vessel, i have read some papers about making sailing model, one of it was talking about metacenter.

I am more likely treating that ship like a body that is moved by external forces.

Basic stability means that the vessel floats upright on a fairly even keel with enough positive stability to return to the upright position if it is disturbed by these external forces.

You were complaining that the vessel pitched into the water and then disappeared below the surface, or fell over in a turn. Neither effect has anything to do with drag.

It suggests, however, that the vessel is unstable, which means that GM is too small or is in fact negative. I am also assuming that no loads have been applied to this model from any forces developed by the sails, which suggests that the vessel is not stable enough to remain upright in any kind of wind.

 

[WiredCat]

ok that's a start, what is GM?

 

[SteamKing]

ok that's a start, what is GM?

That's the metacentric height of the vessel.