Mirror and Free Surface Effect¶
This notebook shows the free surface and mirror functionality and compares the results with experiments.
Free Surface Case¶
To model a free surface we set s = setup(FreeSurface=True,FroudeNumber=Fr) or change s.
The location of the airfoil with respect to the origin is important.
The water surface or ground is the ray through the origin parallel to the free flow (which is parallel to the ground as you know it, if aoa=0).
Because that may be a bit weird to follow for validation and because it makes visualizing things less intuitive, below the airfoil is rotated and everything is simulated with aoa=0.
Regarding validation data: As far as I have seen, the data of Ausman [1] is the best experimental data of a submerged hydrofoil close to the surface - such that the pressure distribution is significantly affected - despite being also the first experiment that measured the effect on the pressure distribution. Now, there is also the big set of Duncan [2], but the NACA0012 profile is rather deeply submerged and the effect on the pressure distribution small [3] (since the topic is not really pressure or effect on hydrofoils but waves and when they collapse which we would not see here since we model linear waves that cannot collapse). I do not have that thesis [1] from the 50s and instead digitized the experimental data as printed in [4]. This is not incidental, the method in [4] is also employed for the panel method in viiflow, with some modifications. Because I do not have the original publication I must guess the Reynolds number. In [5] I found the hint of a "6 inch" model in [1] which together with the Froude number lets me guess the speed (1.25m/s) and thus the Reynolds number (roughly 150000). I also adopted the coordinate system of [5], since their goal was to add to the data by Ausman at higher Froude number. That is, the submergence depth is measured at the trailing edge - but I do not know whether this is correct.
import viiflow as vf
import viiflowtools.vf_tools as vft
import viiflowtools.vf_plots as vfp
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
%matplotlib inline
%config InlineBackend.figure_format = 'svg'
def rotate(deg,center,x):
c, s = np.cos(deg*np.pi/180), np.sin(deg*np.pi/180)
A = np.array([[c, s], [-s, c]])
return ((A@((x.T-center.T).T)).T+center.T).T
# Read Airfoil Data
N = 160
rotmode = True
rotcenter = 1.0
AF=vft.read_selig("NACA4412.dat")
Fr = 1.03
aoa = 5
RE = 1.5e5
AF = vft.repanel(AF,N)
results = {}
heights = np.asarray([.94,.6])
# Mirror codes:
# 0: No mirror
# 1: Classic ground effect mirror
# -1: "Anti mirror" for linearized cp=const free surface case approximation
for h in heights:
results[h] = {} # Sub-Dictionary of results
for ind in [0,1,-1,2]:#[2,1,-1,0]:
# Settings
# In water lower values are more common I hear, like 1-3. Though, from the looks of it, something like .3 would be best for this experiment.
s = vf.setup(Re=RE,Ncrit=.3,Alpha=aoa if not rotmode else 0,Silent=False)
s.FreeSurface = ind==2
# This defines where some velocities are calculated for the plots below, meaning how much chord beyond the foil we want to see.
# It is irrelevant for accuracy
s.FreeSurfMultLeft = 2
s.FreeSurfMultRight = 5
s.FroudeNumber = Fr
s.Mirror = 0 if ind==2 else ind
s.IterateWakes = False
s.WakeLength = 3 # Just to plot below, no relevant for accuracy
results[h][ind] = {} # Sub-Dictionary of results
AFV=rotate(aoa*rotmode,np.asarray([rotcenter,0]),AF.copy())
AFV[1,:]-=h
results[h][ind]["CP"]=[]
results[h][ind]["CPI"]=[]
[p,bl,x] = vf.init(AFV,s)
[_,flag,_,_,_] = vf.iter(x,bl,p,s)
results[h][ind]["AOA"] =aoa
results[h][ind]["CL"] = p.CL
results[h][ind]["CD"] = bl[0].CD
results[h][ind]["CP"].append(p.cp.copy())
results[h][ind]["CPI"].append(p.cp_inviscid.copy())
results[h]["GEO"] = AFV.copy()
results[h]["WAKE"] = p.wakes[0].X.copy()
if s.FreeSurface:
results[h][ind]["SurfaceNodes"] = p.SurfaceNodes[0,:].copy()
results[h][ind]["SurfaceHeight"] = p.SurfaceHeight.copy() # Just -p.FroudeNumber**2*(p.SurfaceVel[0,:]-1)
# Look for nonlinear surface condition
# 1/2|(u,v)-(1,0)|^2+1/Fr^2*dz =0
# It's 1D, we use simple finite-difference Newton iteration
# Don't need to do this of course, but maybe you're interested.
v = np.zeros((2,1))
v2 = np.zeros((2,1))
xp = np.zeros((2,len(p.SurfaceNodes[0,:])))
xn = np.zeros((2,1))
xn2 = np.zeros((2,1))
for j in range(len(p.SurfaceNodes[0,:])):
xn[0,0] = p.SurfaceNodes[0,j]
xn2[0,0] = p.SurfaceNodes[0,j]
delz = 1e-9
for it in range(30): # Actually takes 1-2 iterations
# Evaluate
vf.inviscid_velocity(p,xn[:,0:1],v[:,0:1])
F = 1/2*((v[0])**2+(v[1])**2 - 1)+1/Fr**2*xn[1,0:1]
# Finite difference derivative
xn2[1,0] = xn[1,0]+delz
vf.inviscid_velocity(p,xn2[:,0:1],v2[:,0:1])
dvdz = (v2-v)/delz
dFdz = (v[0])*dvdz[0] + v[1]*dvdz[1] + 1/Fr**2
if abs(dFdz)<1e-14:
break
# Newton step
dz = -F/dFdz
xn[1,0] += dz[0]
if abs(dz)<1e-9:
break
xp[0,j] = xn[0,0]
xp[1,j] = xn[1,0]
results[h][ind]["SurfaceHeightNonlinear"] = xp.copy()
Iteration 22, |res| 0.000069, lam 0.500000 Iteration 19, |res| 0.000024, lam 1.000000 Iteration 17, |res| 0.000063, lam 1.000000 Iteration 15, |res| 0.000019, lam 1.000000 Iteration 22, |res| 0.000069, lam 0.500000 Iteration 20, |res| 0.000044, lam 1.000000 Iteration 20, |res| 0.000081, lam 0.500000 Iteration 16, |res| 0.000072, lam 0.500000
matplotlib.rcParams['figure.figsize'] = [12, 6]
EXPRES=np.genfromtxt("NACA4412HData.csv",delimiter=";")
figCP,axCP = plt.subplots(1,2)
figGEO,axGEO = plt.subplots( 2,1)
ii=0
for h in heights:
i = ii%2
axGEO[i].plot(results[h]["GEO"][0,:],results[h]["GEO"][1,:])
axGEO[i].plot(results[h]["WAKE"][0,:],results[h]["WAKE"][1,:])
axGEO[i].plot(results[h][2]["SurfaceNodes"],results[h][2]["SurfaceHeight"],label="Linear")
axGEO[i].plot(results[h][2]["SurfaceHeightNonlinear"][0,:],results[h][2]["SurfaceHeightNonlinear"][1,:]*1,'-',label="Nonlinear")
axGEO[i].axis("equal")
axGEO[i].set_xlim([-2,6])
axGEO[i].set_title('Depth %g'%h)
axGEO[i].grid(1)
axGEO[i].legend()
#axCP[i].plot(p.foils[0].X[0,:],-results[h][0]["CP"][0],label="Free flow")
axCP[i].plot(p.foils[0].X[0,:],-results[h][1]["CP"][0],'-.',label="Mirror")
axCP[i].plot(p.foils[0].X[0,:],-results[h][-1]["CP"][0],'-.',label="Anti-Mirror")
for kk in range(len(results[h][2]["CP"])):
axCP[i].plot(p.foils[0].X[0,:],-results[h][2]["CP"][kk],label="Free Surface")
axCP[i].plot(p.foils[0].X[0,:],-results[h][2]["CPI"][kk],label="Free Surface (inviscid)")
axCP[i].set_title('Depth %g'%h)
if ii<2:
axCP[i].plot(EXPRES[:,2*i],EXPRES[:,2*i+1],'.',label="Experiment")
axCP[i].legend(loc=4)
axCP[i].grid(1)
ii+=1
Discussion¶
The change in pressure can obviously not be modelled by any mirror case (which are the "high Fr" and "low Fr" limit cases). As for the free surface pressures, the effect of the waves is reasonably well predicted I'd say. The inviscid pressure is in line with other published works (e.g. [3]), the viscous calculation corrects that a bit into the "right" direction, including a transition point one might see in the experimental data. The effect of submergence is a bit underpredicted for the larger submergence depth. In contrast to the following comparison, here the lift coefficient (from surface pressure) seems to be very well modelled at .6 submergence depth.
Lift and Drag for NACA 0012¶
The following case looks at the wave drag and the change in lift in the free surface model in viiflow. The data is from [6], where the other two variants from two other publications have been printed as well and used in the graphs below. Please refer to the publication [6] for more details - interestingly all of them are unsteady RANS comutations which track the water-air interface with some means, but with surprisingly variable results.
# Read Airfoil Data
N = 160
rotmode = True
rotcenter = 1.0
AF=vft.read_selig("NACA0012.dat")
Fr = .571
aoa = 5
RE = 1.59e5
heights = np.r_[np.linspace(.3,1.3,25), np.linspace(1.3,4,10)]
results = {}
results["CL"] = []
results["CD"] = []
for h in heights:
for ind in [2]:
s = vf.setup(Re=RE,Ncrit=.3,Alpha=aoa if not rotmode else 0,Silent=True)
s.FreeSurface = True
s.FroudeNumber = Fr
AFV=rotate(aoa*rotmode,np.asarray([rotcenter,0]),AF.copy())
AFV[1,:]-=h
[p,bl,x] = vf.init(AFV,s)
[_,flag,_,_,_] = vf.iter(x,bl,p,s)
results["CL"].append(p.CL)
results["CD"].append(bl[0].CD+p.CDwave)
# Data in the study [6]
import io
# 1. Store the CSV data as a string
csvURANS = """Submergence_Ratio_h/c,Lift_Coefficient_CL,Drag_Coefficient_CD
4.0,0.53,0.015
3.5,0.53,0.015
3.0,0.53,0.015
2.5,0.54,0.015
2.0,0.54,0.015
1.5,0.56,0.017
1.3,0.59,0.019
1.15,0.61,0.023
1.0,0.63,0.028
0.91,0.63,0.033
0.8,0.60,0.040
0.7,0.40,0.058
0.55,0.20,0.060
0.4, -0.06,0.050
0.3,-0.23,0.037
0.2,-0.32,0.032
0.1,-0.40,0.027"""
csvAli = """Submergence_Ratio_h/c,Lift_Coefficient_CL,Drag_Coefficient_CD
4.0,0.53,0.011
3.0,0.53,0.011
2.5,0.53,0.013
2.0,0.54,0.012
1.5,0.56,0.014
0.91,0.61,0.022"""
csvPrasad = """Submergence_Ratio_h/c,Lift_Coefficient_CL,Drag_Coefficient_CD
1.3,0.60,0.028
1.15,0.62,0.031
1.0,0.63,0.036
0.91,0.64,0.040
0.8,0.61,0.054"""
# 3. Read the data directly into a NumPy array
# skiprows=1 skips the header row
dataURANS = np.loadtxt(io.StringIO(csvURANS), delimiter=",", skiprows=1)
dataAli = np.loadtxt(io.StringIO(csvAli), delimiter=",", skiprows=1)
dataPra = np.loadtxt(io.StringIO(csvPrasad), delimiter=",", skiprows=1)
fig, ax = plt.subplots(1,2)
ax[0].plot(heights,results["CL"],label='Viiflow Linear Free Surface')
ax[1].plot(heights,results["CD"],label='Viiflow Linear Free Surface')
for k in range(1,3):
ax[k-1].plot(dataURANS[:,0],dataURANS[:,k],'d--',label='Pernod et al. [6] (URANS)')
ax[k-1].plot(dataAli[:,0],dataAli[:,k],'d--',label='Ali et al. [6] (URANS)')
ax[k-1].plot(dataPra[:,0],dataPra[:,k],'d--',label='Prasad et al. [6] (URANS)')
ax[k-1].grid(1)
ax[k-1].legend()
ax[k-1].set_xlabel("h/c")
ax[0].set_ylim([-.6,1])
ax[0].set_ylabel("CL")
ax[1].set_ylim([0,.1])
ax[1].set_ylabel("CD")
Text(0, 0.5, 'CD')
Discussion¶
The four methods somewhat agree on the increase in drag and lift. However, the linear wave model in viiflow does not capture wave breakdown and overestimates the increase in lift compared to the URANS methods. The URANS methods also agree on an earlier onset of wave drag increase. The maximum wave drag at .5 seems to be visible in viiflow as well, only to be superseeded by a sharp rise in drag that we do not seen in the URANS data. This is likely due to wave breakdown, but as [7] muses: "The subsequent decrease in drag for h/c < 0.5 is not fully understood. As the amplitude of the following wave train is reduced, it may be conjectured that the wave drag decreases, but some energy is certainly dissipated in the wave breaking". At this point, I also do not understand the dramatic rise in drag below .5 in the linear method.
I would conclude: For submergence depths and Froude numbers far from wave breaking, the increase in drag and lift is quite well modelled by viiflow. For lower submergence, wave breaking may occur and the linear model may overpredict a lift increase, and very shallow depths should probably not be trusted. However, there are no experimental values to compare these different results to.
[1] Ausman, J.S. "Pressure limitation on the upper surface of a hydrofoil." Ph. D. Thesis in Mechanics Engineering at the University of California.
[2] Duncan, James H. "The breaking and non-breaking wave resistance of a two-dimensional hydrofoil." Journal of fluid mechanics 126 (1983): 507-520.
[3] Chen, H., and H. Hefazi. "A panel method for two-dimensional nonlinear free surface flows." WIT Transactions on Engineering Sciences 36 (2002).
[4] Chen, Zhi-min, and W. G. Price. "Dissipative free-surface solver for potential flow around hydrofoil distributed with doublets." Applied Mathematics and Mechanics 33 (2012): 1467-1480.
[5] Parkin, Blaine R., Byrne Perry, and T. Yao‐tsu Wu. "Pressure distribution on a hydrofoil running near the water surface." Journal of Applied physics 27.3 (1956): 232-240.
[6] Pernod, Laetitia, et al. "Free-surface effects on two-dimensional hydrofoils by RANS-VOF simulations." Journal of Sailing Technology 8.01 (2023): 24-38.