Fix wave current velocity (WEC-Sim#1540)

akeeste · web-flow · commit eddcfd23a7dd · 2025-08-23T16:25:14.000-06:00
* correct current velocity in noWave case

* correct current velocity in regular case

* correct current velocity in irregular case
diff --git a/source/functions/simulink/model/irregWaveMorison.m b/source/functions/simulink/model/irregWaveMorison.m
@@ -46,7 +46,7 @@
         curramp        = currentSpeedDepth;
     end
     %Vel should be a column vector
-    Vel2    = [Vel(1),Vel(2),Vel(3)] + wxr + [curramp*cosd(currentDirection),curramp*sind(currentDirection),0];
+    Vel2    = [Vel(1),Vel(2),Vel(3)] + wxr;
     % Update translational and rotational acceleration
     % dotw refers to \dot{\omega} = rotational acceleration
     Accelt  = [Accel(4),Accel(5),Accel(6)];
@@ -58,7 +58,7 @@
     %% Calculate Directionality
     for aa = 1:length(direction)
 
-        %% Calculate Orbital Velocity
+        %% Calculate Fluid Velocity and Acceleration
         for jj = 1:ff
             waveDirRad      = direction(aa)*pi/180;
             phaseArg        = w(jj,1)*Time - k(jj,1)*(ShiftCg(1)*cos(waveDirRad) + ShiftCg(2)*sin(waveDirRad)) + randPhase(jj,1);
@@ -92,8 +92,9 @@
 
         % Sum the wave components to obtain the x, y, z orbital velocities
         uV = sum(uVt); uA   = sum(uAt); vV = sum(vVt); vA = sum(vAt); wV = sum(wVt); wA = sum(wAt);
-        fluidV              = [uV, vV, wV];
-        fluidA              = [uA, vA, wA];
+        % Total fluid velocity and acceleration
+        fluidV          = [uV, vV, wV] + [curramp*cosd(currentDirection) curramp*sind(currentDirection) 0];
+        fluidA          = [uA, vA, wA];
         if bodyMorison == 2
             %% Decompose Fluid Velocity
             % Tangential Velocity
@@ -133,22 +134,30 @@
             areaRot         = abs(mtimes(Area(ii,:),RotMax));
             CdRot           = mtimes(abs(Cd(ii,:)),RotMax);
             CaRot           = abs(mtimes(Ca(ii,:),RotMax));
+            % Fluid velocity relative to the body
+            uVdiff          = fluidV(1) - Vel2(1);
+            vVdiff          = fluidV(2) - Vel2(2);
+            wVdiff          = fluidV(3) - Vel2(3);
+            % Fluid acceleration relative to the body
+            uAdiff          = fluidA(1) - Accel2(1);
+            vAdiff          = fluidA(2) - Accel2(2);
+            wAdiff          = fluidA(3) - Accel2(3);
             % Forces from velocity drag
-            uVdiff          = uV - Vel2(1); FxuV = (1/2)*abs(uVdiff)*uVdiff*rho*CdRot(1)*areaRot(1);
-            vVdiff          = vV - Vel2(2); FxvV = (1/2)*abs(vVdiff)*vVdiff*rho*CdRot(2)*areaRot(2);
-            wVdiff          = wV - Vel2(3); FxwV = (1/2)*abs(wVdiff)*wVdiff*rho*CdRot(3)*areaRot(3);
+            FxuV            = (1/2)*abs(uVdiff)*uVdiff*rho*CdRot(1)*areaRot(1);
+            FxvV            = (1/2)*abs(vVdiff)*vVdiff*rho*CdRot(2)*areaRot(2);
+            FxwV            = (1/2)*abs(wVdiff)*wVdiff*rho*CdRot(3)*areaRot(3);
             % Forces from body acceleration inertia
-            uAdiff          = uA - Accel2(1); FxuA = uAdiff*rho*Vol(ii,:)*CaRot(1);
-            vAdiff          = vA - Accel2(2); FxvA = vAdiff*rho*Vol(ii,:)*CaRot(2);
-            wAdiff          = wA - Accel2(3); FxwA = wAdiff*rho*Vol(ii,:)*CaRot(3);
+            FxuA            = uAdiff*rho*Vol(ii,:)*CaRot(1);
+            FxvA            = vAdiff*rho*Vol(ii,:)*CaRot(2);
+            FxwA            = wAdiff*rho*Vol(ii,:)*CaRot(3);
             % Forces from fluid acceleration inertia
             FxuAf           = uA*Vol(ii,:)*rho;
             FxvAf           = vA*Vol(ii,:)*rho;
             FxwAf           = wA*Vol(ii,:)*rho;
             % Total Force from the Morison Element
             F               = [FxuV + FxuA + FxuAf,...
-                FxvV + FxvA + FxvAf,...
-                FxwV + FxwA + FxwAf];
+                               FxvV + FxvA + FxvAf,...
+                               FxwV + FxwA + FxwAf];
         end
 
         FMd(aa,:) = F;
diff --git a/source/functions/simulink/model/noWaveMorison.m b/source/functions/simulink/model/noWaveMorison.m
@@ -42,13 +42,25 @@
         curramp        = currentSpeedDepth;
     end
     %Vel should be a column vector
-    Vel2    = [Vel(1),Vel(2),Vel(3)] + wxr + [curramp*cosd(currentDirection),curramp*sind(currentDirection),0];
+    Vel2    = [Vel(1),Vel(2),Vel(3)] + wxr;
     % Update translational and rotational acceleration
     % dotw refers to \dot{\omega} = rotational acceleration
     Accelt  = [Accel(4),Accel(5),Accel(6)];
     dotwxr  = cross(Accelt,r(ii,:));
     wxwxr   = cross(Velt,wxr);
     Accel2  = [Accel(1),Accel(2),Accel(3)] + dotwxr + wxwxr;
+    %% Calculate Fluid Velocity and Acceleration
+    % Orbital Velocity
+    uV              = 0;
+    vV              = 0;
+    wV              = 0;
+    % Orbital Acceleration
+    uA              = 0;
+    vA              = 0;
+    wA              = 0;
+    % Total fluid velocity and acceleration
+    fluidV          = [uV, vV, wV] + [curramp*cosd(currentDirection) curramp*sind(currentDirection) 0];
+    fluidA          = [uA, vA, wA];
     if bodyMorison == 2
         %% Decompose Body Velocity
         % Tangential Velocity
@@ -74,14 +86,22 @@
         areaRot     = abs(mtimes(Area(ii,:),RotMax));
         CdRot       = mtimes(abs(Cd(ii,:)),RotMax);
         CaRot       = abs(mtimes(Ca(ii,:),RotMax));
-        % Forces from body velocity drag
-        FxuV        = (1/2)*abs(-1*Vel2(1))*-1*Vel2(1)*rho*CdRot(1)*areaRot(1);
-        FxvV        = (1/2)*abs(-1*Vel2(2))*-1*Vel2(2)*rho*CdRot(2)*areaRot(2);
-        FxwV        = (1/2)*abs(-1*Vel2(3))*-1*Vel2(3)*rho*CdRot(3)*areaRot(3);
+        % Fluid velocity relative to the body
+        uVdiff      = fluidV(1) - Vel2(1);
+        vVdiff      = fluidV(2) - Vel2(2);
+        wVdiff      = fluidV(3) - Vel2(3);
+        % Fluid acceleration relative to the body
+        uAdiff      = fluidA(1) - Accel2(1);
+        vAdiff      = fluidA(2) - Accel2(2);
+        wAdiff      = fluidA(3) - Accel2(3);
+        % Forces from velocity drag
+        FxuV        = (1/2)*abs(uVdiff)*uVdiff*rho*CdRot(1)*areaRot(1);
+        FxvV        = (1/2)*abs(vVdiff)*vVdiff*rho*CdRot(2)*areaRot(2);
+        FxwV        = (1/2)*abs(wVdiff)*wVdiff*rho*CdRot(3)*areaRot(3);
         % Forces from body acceleration inertia
-        FxuA        = rho*Vol(ii,:)*CaRot(1)*-1*Accel2(1);
-        FxvA        = rho*Vol(ii,:)*CaRot(2)*-1*Accel2(2);
-        FxwA        = rho*Vol(ii,:)*CaRot(3)*-1*Accel2(3);
+        FxuA        = uAdiff*rho*Vol(ii,:)*CaRot(1);
+        FxvA        = vAdiff*rho*Vol(ii,:)*CaRot(2);
+        FxwA        = wAdiff*rho*Vol(ii,:)*CaRot(3);
         % Sum the force contributions
         F           = [FxuV + FxuA,FxvV + FxvA,FxwV + FxwA];
     end
diff --git a/source/functions/simulink/model/regWaveMorison.m b/source/functions/simulink/model/regWaveMorison.m
@@ -43,14 +43,14 @@
         curramp     = currentSpeedDepth;
     end
     %Vel should be a column vector
-    Vel2            = [Vel(1),Vel(2),Vel(3)] + wxr + [curramp*cosd(currentDirection),curramp*sind(currentDirection),0];
+    Vel2            = [Vel(1),Vel(2),Vel(3)] + wxr;
     % Update translational and rotational acceleration
     % dotw refers to \dot{\omega} = rotational acceleration
     Accelt          = [Accel(4),Accel(5),Accel(6)];
     dotwxr          = cross(Accelt,r(ii,:));
     wxwxr           = cross(Velt,wxr);
     Accel2          = [Accel(1),Accel(2),Accel(3)] + dotwxr + wxwxr;
-    %% Calculate Orbital Velocity
+    %% Calculate Fluid Velocity and Acceleration
     waveDirRad      = direction*pi/180;
     phaseArg        = w*Time - k*(ShiftCg(1)*cos(waveDirRad) + ShiftCg(2)*sin(waveDirRad));
     % Vertical Variation
@@ -73,11 +73,12 @@
     uV              =  ramp*coeffHorz*cos(phaseArg)*g*k*(1/w)*cos(waveDirRad);
     vV              =  ramp*coeffHorz*cos(phaseArg)*g*k*(1/w)*sin(waveDirRad);
     wV              = -ramp*coeffVert*sin(phaseArg)*g*k*(1/w);
-    fluidV          = [uV, vV, wV];
     % Orbital Acceleration
     uA              = -ramp*coeffHorz*sin(phaseArg)*g*k*cos(waveDirRad);
     vA              = -ramp*coeffHorz*sin(phaseArg)*g*k*sin(waveDirRad);
     wA              = -ramp*coeffVert*cos(phaseArg)*g*k;
+    % Total fluid velocity and acceleration
+    fluidV          = [uV, vV, wV] + [curramp*cosd(currentDirection) curramp*sind(currentDirection) 0];
     fluidA          = [uA, vA, wA];
     if bodyMorison == 2
         %% Decompose Fluid Velocity
@@ -118,14 +119,22 @@
         areaRot     = abs(mtimes(Area(ii,:),RotMax));
         CdRot       = mtimes(abs(Cd(ii,:)),RotMax);
         CaRot       = abs(mtimes(Ca(ii,:),RotMax));
+        % Fluid velocity relative to the body
+        uVdiff      = fluidV(1) - Vel2(1);
+        vVdiff      = fluidV(2) - Vel2(2);
+        wVdiff      = fluidV(3) - Vel2(3);
+        % Fluid acceleration relative to the body
+        uAdiff      = fluidA(1) - Accel2(1);
+        vAdiff      = fluidA(2) - Accel2(2);
+        wAdiff      = fluidA(3) - Accel2(3);
         % Forces from velocity drag
-        uVdiff      = uV - Vel2(1); FxuV = (1/2)*abs(uVdiff)*uVdiff*rho*CdRot(1)*areaRot(1);
-        vVdiff      = vV - Vel2(2); FxvV = (1/2)*abs(vVdiff)*vVdiff*rho*CdRot(2)*areaRot(2);
-        wVdiff      = wV - Vel2(3); FxwV = (1/2)*abs(wVdiff)*wVdiff*rho*CdRot(3)*areaRot(3);
+        FxuV        = (1/2)*abs(uVdiff)*uVdiff*rho*CdRot(1)*areaRot(1);
+        FxvV        = (1/2)*abs(vVdiff)*vVdiff*rho*CdRot(2)*areaRot(2);
+        FxwV        = (1/2)*abs(wVdiff)*wVdiff*rho*CdRot(3)*areaRot(3);
         % Forces from body acceleration inertia
-        uAdiff      = uA - Accel2(1); FxuA = uAdiff*rho*Vol(ii,:)*CaRot(1);
-        vAdiff      = vA - Accel2(2); FxvA = vAdiff*rho*Vol(ii,:)*CaRot(2);
-        wAdiff      = wA - Accel2(3); FxwA = wAdiff*rho*Vol(ii,:)*CaRot(3);
+        FxuA        = uAdiff*rho*Vol(ii,:)*CaRot(1);
+        FxvA        = vAdiff*rho*Vol(ii,:)*CaRot(2);
+        FxwA        = wAdiff*rho*Vol(ii,:)*CaRot(3);
         % Forces from fluid acceleration inertia
         FxuAf       = uA*Vol(ii,:)*rho;
         FxvAf       = vA*Vol(ii,:)*rho;
