validation section, errors on vertical shear

internalLoads.py

@@ -109,6 +109,7 @@ pvd_name = 'patchFields/'+\
     
     # +'vtu/patchFields_'+ \
     # '0.0994949115796907729'+'.vtu'
+    
 pvd_name = os.path.join(folder, pvd_name)
 
 time_foamstar=[]
@@ -150,9 +151,9 @@ with open(dissMassFile, 'r') as f:
                 count = 0
         k += 1
 
-sections['0']['REFPOINT'][0] = 100
+# sections['0']['REFPOINT'][0] = 100
 
-print('Initialisation : reading cell centers at first time step \n')
+print('Initialisation : reading cell centers within section cut at first time step \n')
 
 file_name = 'patchFields/'+\
  'vtu/body_'+ \
@@ -182,19 +183,7 @@ angles = motionData[ini,4:7].reshape(3)
 # motion in R0 and angles in Rb
 CoR = motion + CoRInitial # CoR in R0
 
-# faceIDs, faceCfs=cellPoints(output, CoR, angles) # TODO check reference frame transformation, start from fixed and transform total force
-
-# alpha=np.array(output.GetCellData().GetArray("alpha.water")) ## les valeurs au centre des cellules
-
-print('Initialisation : selection of cells within section cut')
-
-    
-##Conversion au point des cellules
-# converter = vtk.vtkCellDataToPointData()    
-# converter.SetInputConnection(reader.GetOutputPort() )
-# converter.Update()
 
-# PressurePoints = np.array(converter.GetOutput().GetPointData().GetArray("p")) ##Valeurs points cellules
 
 # faceIDs = np.loadtxt('tmp/faceIDs.txt')
 # faceWeights = np.loadtxt('tmp/faceWeights.txt')
@@ -233,7 +222,12 @@ time = []
 shearSection = np.zeros(( len(time_foamstar) , 3 ))
 momentSection = np.zeros(( len(time_foamstar) , 3 ))
 
-IACS = False
+fluidForces = np.zeros(( len(time_foamstar) , 3 ))
+fluidMoments = np.zeros(( len(time_foamstar) , 3 ))
+
+IACS = True
+ShearStress = False
+reCalcBodyAcc = False
 
 print('\n Cut point : ' + str(cutPoint) + '\n')
 
@@ -256,7 +250,9 @@ for dt in range(0,len(time_foamstar)):
     output = reader.GetOutput()
     
     Pressure=np.array(output.GetCellData().GetArray("p")) ## les valeurs au centre des cellules
-    # Shear=np.array(output.GetCellData().GetArray("wallShearStress")) ## les valeurs au centre des cellules
+    
+    if ShearStress:
+        Shear=np.array(output.GetCellData().GetArray("wallShearStress")) ## les valeurs au centre des cellules
     
     ii = np.where(motionData[:,0]==float(time_foamstar[dt]))
 
@@ -273,26 +269,33 @@ for dt in range(0,len(time_foamstar)):
     
     CoR = motion+np.dot(P0b,CoRInitial)
     
-    fluidForces = np.zeros(3)
-    fluidMoments = np.zeros(3)
+    # fluidForces = np.zeros(3)
+    # fluidMoments = np.zeros(3)
     for iface in range(0,len(faceIDs)):
         faceID = faceIDs[iface]
         
-        # Shear[faceID,:] = np.dot( np.transpose(P0b),\
-        #                   Shear[faceID,:] ) # TODO : check if wall shear stress is given in R0 or Rb
+        if ShearStress:
+            Shear[faceID,:] = np.dot( np.transpose(P0b),\
+                              Shear[faceID,:] ) # TODO : check if wall shear stress is given in R0 or Rb
+            
+        
+
         f1 = Pressure[faceID]*faceSf[iface,:]*faceWeights[iface] # TODO : check connectivity of Pressure field and faceIDs
-        # f2 = Shear[faceID,:]*faceWeights[iface]
         f2 = np.zeros(3)
-        fluidForces += f1 + f2
-        fluidMoments += np.cross( faceCfs[iface], f1 ) + np.cross( faceCfs[iface], f2 )
         
-    fluidForces[0]  *= 2
-    fluidForces[1]   = 0
-    fluidForces[2]  *= 2
+        if ShearStress:
+            f2 = Shear[faceID,:]*faceWeights[iface]
+
+        fluidForces[dt,:] += f1 + f2
+        fluidMoments[dt,:] += np.cross( faceCfs[iface], f1 ) + np.cross( faceCfs[iface], f2 )
+        
+    fluidForces[dt,0]  *= 2
+    fluidForces[dt,1]   = 0
+    fluidForces[dt,2]  *= 2
     
-    fluidMoments[0]  = 0
-    fluidMoments[1] *= 2
-    fluidMoments[2]  = 0
+    fluidMoments[dt,0]  = 0
+    fluidMoments[dt,1] *= 2
+    fluidMoments[dt,2]  = 0
     
     ####
     
@@ -309,41 +312,47 @@ for dt in range(0,len(time_foamstar)):
     
     mg = np.dot( inertiaSection[:3,:3] , np.dot( np.transpose(P0b),\
                       g ).reshape(3) )
-    
-    ####
-        
-    # mg = massSection*np.dot( np.transpose(P0b),\
-    #                       g ).reshape(3)
-        
-    # TODO : check cogSection, i.e. Rb or R0
-    # mgMoment = massSection * np.dot( lib.skew_matrix(cogSection), np.dot( np.transpose(P0b),\
-    #                       g ) ).reshape(3)
         
-    # inertia = np.dot( P0b, massSection*linearAcc)
-    # inertia -= np.dot( m_uc, angAcc )
-    
-    # ftheta1 = np.dot( np.transpose(-np.dot( P0b, m_uc )) , linearAcc ) # TODO check
-    # ftheta2 = np.dot( inertiaSection[3:,3:] , angAcc )
-    
-    # inertiaAng = ftheta1+ftheta2
-    
-    # Qvtheta = -np.dot( lib.skew_matrix( omg ), np.dot( inertiaSection[3:,3:] , omg ) )
-    
-    # ftheta = inertiaAng - Qvtheta
-    
-    # Qv = np.dot( lib.skew_matrix( np.dot( P0b , omg) ), np.dot( P0b, m_uc) )
-    # Qv = np.dot( Qv, omg )
-    # Qv = -np.dot( np.transpose(P0b), Qv )
     
     ####
     
     time = np.append(time, float(time_foamstar[dt]))
     ##
-    shearSection[dt,:] = fluidForces #+ mg - np.dot( np.transpose(P0b) , inertia[:3] )
+    shearSection[dt,:] = fluidForces[dt,:] + mg - np.dot( np.transpose(P0b) , inertia[:3] )
     ##
-    momentSection[dt,:] = fluidMoments + np.dot( lib.skew_matrix(cogSection) , mg ) -\
+    momentSection[dt,:] = fluidMoments[dt,:] + np.dot( lib.skew_matrix(cogSection) , mg ) -\
         inertia[3:] - np.dot( lib.skew_matrix(uc) , shearSection[dt] )
         
+    ####
+    
+    if reCalcBodyAcc: # TODO correction on this part
+        
+        tol = 1e-6
+        maxiter = 1000
+        error = 1
+        iteration = 1
+        while error > tol and iteration < maxiter:
+            
+            aux1 = np.dot( P0b ,  - shearSection[dt,:] + fluidForces[dt,:] + mg)
+            aux2 = fluidMoments[dt,:] + np.dot( lib.skew_matrix(cogSection) , mg ) -\
+                np.dot( lib.skew_matrix(uc) , shearSection[dt] ) - momentSection[dt,:]
+                
+            aux = np.hstack(( aux1 , aux2 ))
+            aux = aux + Qv
+            AccRecalc = np.linalg.solve( inertiaSection , aux )
+            
+            inertia = np.dot( inertiaSection , AccRecalc ) - Qv
+            shearSection[dt,:] = fluidForces[dt,:] + mg - np.dot( np.transpose(P0b) , inertia[:3] )
+            ##
+            momentSection[dt,:] = fluidMoments[dt,:] + np.dot( lib.skew_matrix(cogSection) , mg ) -\
+            inertia[3:] - np.dot( lib.skew_matrix(uc) , shearSection[dt] )
+            
+            error = np.linalg.norm(Acc-AccRecalc)
+            Acc = AccRecalc
+            print('Recalculate Body Acceleration error :'+str(error))
+            iteration += 1
+            
+        
         
     ## IACS convention
     
@@ -367,66 +376,134 @@ ax = fig.add_subplot(111, projection='3d')
 ax.plot3D(faceCfs[:,0], faceCfs[:,1], faceCfs[:,2], 'bo')
 ax.set_box_aspect([ub - lb for lb, ub in (getattr(ax, f'get_{a}lim')() for a in 'xyz')])
 
-# res = os.path.join(folder, 'fz.dat')
-# res = np.loadtxt(res, skiprows=4)
+res = os.path.join(folder, 'fz.dat')
+res = np.loadtxt(res, skiprows=4)
 
-# t = res[:, 0]
-# fz = res[:,5]
+t = res[:, 0]
+fz = res[:,5]
 
-# plt.figure()
-# plt.plot(time, -shearSection[:,2], 'r--')
-# plt.plot(t, fz, 'k--')
-# plt.xlim(min(time), max(time))
-# # plt.ylim(-300, 300)
-# plt.grid()
-# plt.xlabel('Time (s)')
-# plt.ylabel('Shear (N)')
+plt.figure()
+plt.plot(time, shearSection[:,2], 'r--')
+plt.plot(t, fz, 'k--')
+plt.xlim(min(time), max(time))
+plt.ylim(-200, 200)
+plt.grid()
+plt.xlabel('Time (s)')
+plt.ylabel('Shear z (N)')
+plt.legend(['internalLoads-python','foamStar'])
 
+####
 
-# res = os.path.join(folder, 'my.dat')
-# res = np.loadtxt(res, skiprows=4)
+res = os.path.join(folder, 'fx.dat')
+res = np.loadtxt(res, skiprows=4)
+
+t = res[:, 0]
+fx = res[:,5]
+
+plt.figure()
+plt.plot(time, shearSection[:,0], 'r--')
+plt.plot(t, fx, 'k--')
+plt.xlim(min(time), max(time))
+# plt.ylim(-200, 200)
+plt.grid()
+plt.xlabel('Time (s)')
+plt.ylabel('Shear x (N)')
+plt.legend(['internalLoads-python','foamStar'])
 
-# t = res[:, 0]
-# my = res[:,5]
 
-# plt.figure()
-# plt.plot(time, momentSection[:,1], 'r--')
-# plt.plot(t, -my, 'k--')
-# plt.xlim(min(time), max(time))
-# plt.ylim(-600, 400)
-# plt.grid()
-# plt.xlabel('Time (s)')
-# plt.ylabel('Moment (N.m)')
 
-res = os.path.join(folder, 'forces.dat')
-res = np.loadtxt(res, skiprows=3)
+####
 
-t = res [:, 0]
-forces = res[:,1:4]
 
-for i in range(0, len(t)):
+res = os.path.join(folder, 'my.dat')
+res = np.loadtxt(res, skiprows=4)
+
+t = res[:, 0]
+my = res[:,5]
+
+plt.figure()
+plt.plot(time, momentSection[:,1], 'r--')
+plt.plot(t, my, 'k--')
+plt.xlim(min(time), max(time))
+plt.ylim(400, -200)
+plt.grid()
+plt.xlabel('Time (s)')
+plt.ylabel('Vertical Bending Moment VBM (N.m)')
+plt.legend(['internalLoads-python','foamStar'])
+
+
+########################################
+
+res = os.path.join(folder, 'fFluid.dat')
+res = np.loadtxt(res, skiprows=4)
+
+t = res[:, 0]
+forces = res[:,13:16]
+
+fluidForcesR0 = np.zeros(( len(time_foamstar) , 3 ))
+
+for i in range(0, len(time_foamstar)):
     
-    angles    = motionData[i,4:7].reshape(3)
+    ii = np.where(motionData[:,0]==float(time_foamstar[i]))
+    
+    angles    = motionData[ii,4:7].reshape(3)
 
     P0b = lib.rotation_matrix(angles[0],angles[1],angles[2])
     
-    forces[i,:] = np.dot( np.transpose(P0b) , forces[i,:] )
+    fluidForces[i,:] = np.dot( P0b , fluidForces[i,:] )
+    
 
 plt.figure()
-plt.plot(time, shearSection[:,2]/2, 'r--')
+plt.plot(time, fluidForces[:,2], 'r--')
 plt.plot(t, forces[:,2], 'k--')
 plt.xlim(min(time), max(time))
+plt.ylim(1000, 1500)
 plt.grid()
 plt.xlabel('Time (s)')
-plt.ylabel('Force Z (N)')
+plt.ylabel('Fluid force Z (N)')
+
 
 plt.figure()
-plt.plot(time, shearSection[:,0]/2, 'r--')
+plt.plot(time, fluidForces[:,0], 'r--')
 plt.plot(t, forces[:,0], 'k--')
 plt.xlim(min(time), max(time))
+plt.ylim(0, 200)
 plt.grid()
 plt.xlabel('Time (s)')
-plt.ylabel('Force X (N)')
+plt.ylabel('Fluid force X (N)')
+
+
+########################################
+
+# res = os.path.join(folder, 'forces.dat')
+# res = np.loadtxt(res, skiprows=3)
+
+# t = res [:, 0]
+# forces = res[:,1:4]
+
+# for i in range(0, len(t)):
+    
+#     angles    = motionData[i,4:7].reshape(3)
+
+#     P0b = lib.rotation_matrix(angles[0],angles[1],angles[2])
+    
+#     forces[i,:] = np.dot( np.transpose(P0b) , forces[i,:] )
+
+# plt.figure()
+# plt.plot(time, shearSection[:,2]/2, 'r--')
+# plt.plot(t, forces[:,2], 'k--')
+# plt.xlim(min(time), max(time))
+# plt.grid()
+# plt.xlabel('Time (s)')
+# plt.ylabel('Force Z (N)')
+
+# plt.figure()
+# plt.plot(time, shearSection[:,0]/2, 'r--')
+# plt.plot(t, forces[:,0], 'k--')
+# plt.xlim(min(time), max(time))
+# plt.grid()
+# plt.xlabel('Time (s)')
+# plt.ylabel('Force X (N)')
 
 
 

libs.py

@@ -228,7 +228,6 @@ def find_and_keep_faces(output, cutPoint, CoR, angles):
     # zCut = cutPoint[2]  # Assuming cutPoint is a 3D vector
     xCut = cutPoint[0]  # Assuming cutPoint is a 3D vector
 
-
     # Lists to store face IDs and weights for each patch
     faceWeights = []
     faceIDs = []
@@ -262,7 +261,7 @@ def find_and_keep_faces(output, cutPoint, CoR, angles):
             faceIDs=np.append(faceIDs, faceI)
             
             if count == cellIds.GetNumberOfIds():
-    
+                
                 faceWeights=np.append(faceWeights, 1.0)
                 
                 center = center/cellIds.GetNumberOfIds()