First of all define a boat hull, lets say you define it by triangles, now you cut those triangles with waterplanes to get the area that is under water, you list those triangles (newly created), then you calculate skin drag and form drag and apply that to rotation velocity, after that you need a damping function either way to cut the close to zero related problems and smoothly deaccelerate things for the simulation to run efficently.
Now pay attention
vec3 form_drag;
vec3 skin_drag;
float whole_sub_surf = 0.0;
for (int i=0; i < submerged_tri_count; i++)
//now iterate through form_percentage_coefff
if (ACTUAL_FRAME[i].surface > 0.1) //do not coun't anything that is less than 10 cm^2
{
whole_sub_surf = whole_sub_surf + ACTUAL_FRAME[i].surface;
vec3 element_form_drag = vec3(0.0, 0.0, 0.0); //somehow i dont 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 = 0; //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 = 1.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;
result_force = gravity_force + buoyancy_force + drag_force + thrust_vec;// + rudder_force;// + drag_force + thrust_vec + rudder_force;
mat4 inv_rot = ROTATION_MAT;
inv_rot.Inverse();
vec3 local_torque = inv_rot * elements_torque;
EAngVel = EAngAcc*dt;
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();
float yawrate = 5.0*dt;
YPRangle.yaw(cos(yawrate*imopi), sin(yawrate*imopi));
YPRangle.DoRotation();
ROTATION_MAT.LoadGLMatrix(YPRangle.AIR_MATRIX);
BUO_COG = vectorAB(CENTER_OF_BUOYANCY, pos + ROTATION_MAT * CENTER_OF_GRAVITY);
vec3 accel = result_force / mass;
vel = vel + accel*dt;
pos = pos + vel*dt;
PC stands for pressure center point
epsilon = 0.001;
inline void damp( vec3 * ptr, float amount, float epsilon)
{
(*ptr) = (*ptr) - (*ptr)/amount;
vec3 p = vec3(absnf<float>((*ptr).x),absnf<float>((*ptr).y),absnf<float>((*ptr).z));
if (IsNan(p.x)) p.x = 0.0;
if (IsNan(p.y)) p.y = 0.0;
if (IsNan(p.z)) p.z = 0.0;
if ( (p.x > 0) && (p.x < epsilon) ) (*ptr).x = 0.0;
if ( (p.y > 0) && (p.y < epsilon) ) (*ptr).y = 0.0;
if ( (p.z > 0) && (p.z < epsilon) ) (*ptr).z = 0.0;
}
inline void damp( float * ptr, float amount, float epsilon)
{
(*ptr) = (*ptr) - (*ptr)/amount;
if (IsNan((*ptr))) (*ptr) = 0.0;
if ( (absnf<float>(*ptr) > 0) && (absnf<float>(*ptr) < epsilon) ) (*ptr) = 0.0;
}
//---------------------------------------------------------------------------
#ifndef aeroH
#define aeroH
//---------------------------------------------------------------------------
//Mix of crappyness
//don't even think it is correct in 0.0001%
#include <math.h>
#include "DxcMath.h"
const float AIR_STD_KINEMATIC_VISCOSITY = 1.460 * pow(10.0, -5.0); // m^2/s
const float WATER_STD_KINEMATIC_VISCOSITY = 1.3083 * pow(10.0, -6.0); // m^2/s at temp = 10*C
inline float ReynoldsNum(float V, float choord_len, float kinematic_viscosity_of_fluid)
{
return (V*choord_len) / kinematic_viscosity_of_fluid;
}
//0.894
inline vec3 FindTriangleAerodynamicCenter(vec3 A, vec3 B, vec3 C)
{
}
inline float getChoordLength(vec3 velocity, vec3 surf_normal, vec3 pC, vec3 * verts, int vert_count, bool project, vec3 point_on_surf)
{
//Project velocity vector onto surface
vec3 projected = Normalize(ProjectVectorOnPlane(surf_normal, velocity));
//find the biggest side and square it - this will ensure us that ray will always hit two sides
float adst = -999.0;
for (int i=0; i < vert_count; i++)
{
int next = i + 1;
if (next == vert_count) next = 0;
adst = Max(adst, n3ddistance(verts[i],verts[next]));
}
adst = adst*adst;
//compute ray
vec3 rA = pC - projected*adst;
vec3 rB = pC + projected*adst;
vec3 cp;
for (int i=0; i < vert_count; i++)
cp = cp + verts[i];
cp = cp / float(vert_count);
vec3 n = surf_normal*1000.0;
vec3 colp;
//we compute center point to determine whenever a side face is at front or back to the center point (because i don't need to waste time and run math on paper i can just check that and use it for further processing)
int cpside;
for (int i=0; i < vert_count; i++)
{
int next = i + 1;
if (next == vert_count) next = 0;
vec3 normal = Normal(verts[i], verts[next], verts[next] + n, true);
float dst = getplaneD(normal, verts[i]);
//cpside will now determine the
cpside = classifyapointagainstaplane(cp, normal, dst); // which side faces 'the inside' of polygon (always either front or back never on plane) - solid convex polygon
if (cpside != isBack) //when center point is not behind side plane
{
normal = -normal;
dst = getplaneD(normal, verts[i]);
}
colp = rA;
if (SegmentPlaneIntersection(normal, dst, rA, rB, colp))
{
//if ray and face normal 'look' (face) at each other
if (dot(Normalize(vectorAB(rA,rB)), normal) < 0) rA = colp;
else rB = colp;
}
}
return n3ddistance(rA, rB);
}
inline float getTriangleChoordLength(vec3 vel, vec3 surf_n, vec3 pC, vec3 A, vec3 B, vec3 C, bool project, vec3 pop)
{
vec3 ahue[3];
ahue[0] = A;
ahue[1] = B;
ahue[2] = C;
return getChoordLength(vel, surf_n, pC, ahue, 3, project, pop);
}
inline vec3 ShearStress(vec3 Force, float Area)
{
return Force / Area;
}
//1/ (pressure*v^2*Area) * ShearStress * dot(surf_normal, tangent_surf_dir)
inline float SkinFrictionDrag( float dens, float V, float Area, float ShearStress, t3dpoint<double> surf_normal, t3dpoint<double> surf_tangent_relative_to_airflow)
{
return (dens*V*V*Area) * ShearStress * absnf(dot(surf_normal, surf_tangent_relative_to_airflow));
}
inline float BodyDrag(float V, float Area, float fluidDens, t3dpoint<double> SurfNormal, t3dpoint<double> FluidFlowDir, bool front_only)
{
//now im not sure if its only the percentage of the rotation or if body is rotated the area increases
//anyways im not increasing the area.
float ratio = -dot(FluidFlowDir, SurfNormal);
if(!front_only) ratio = absnf(ratio);
return (0.5*(V*V)*fluidDens*Area*ratio);
}
#endif
And heres a vid but keep in mind i applied rotation force and thrust force each frame so they are moving almost constantly
https://youtu.be/7YXtG73KY_k
Also the magic number to be applied might need to vary with the size of the boat, I'll have to test it out.. and I could turn off the righting if the boat has capsized.
