Cavitation¶
Cavitation occurs when the pressure drops below vapour pressure of water, which then turns into gas. The cavitation solver looks for states in the flow where such gas bubbles can exist and change the flow because the pressure on the foil is at vapour pressure and there is no pressure below vapour pressure anywhere on the foil.
There is no actual time-dynamic resolution of the bubble-forming process, nor is the transport or decay of these bubbles modelled. As such, these results can give estimates of the effect on the mean force coefficients, but cannot tell us anything about the time-dynamic behaviour.
The watertunnel data is from [1], and I consider it the best available experimental cavitation data. Notes on the experimental results: As discussed in the paper, the lift slope is below the high quality wind tunnel data, and the moment coefficient is quite off as well. No corrections for possible aspect ratio effects or wall interference have been applied to the experiments, and I will not try to do that here. XFOIL and viiflow are closer to the wind tunnel measurements, and to make this comparison work all simulations are set using an angle of attack matching the non-cavitating lift.
In addition looking at the images of the cavity on the foil section in the reference clearly shows how the cavity does not look the same over the span.
Below is the computation and the results in the form of a comparison of all force coefficients at a range of angles of attack, and we see the effect of reducing sigma is reproduced very well. In fact, it seems such a broad and well reproduced comparison is unmatched in the literature.
However, there are things to note:
The computational results can be "noisy" because there can be multiple solutions. These change depending on how finely one resolves the foil near the cavitation bubble end, as well as the weight factor one uses. Especially the moment coefficient (which depends the most on the exact length of the cavity) is affected. This divergence between branches is - seemingly - confined for the same aoa to a region around maximum lift. Do these appear because these are particularly unsteady conditions? Do the different results compare to different points in an unsteady simulation? The author does not know. If one were to use these computations for any subsequent optimization or simulation process - better sample with different resolutions and decide on a sensible strategy for this noise.
These computations are more demanding than a regular viiflow solve. There are more internal variables (about 2.5x) and each internal step grows with N^4, which tests my patience for above 100 panels. For the sake of this example, they are set low (and I'd say sufficient to estimate the effect of cavitation on foils), but for a publication one might increase this of course.
Low sigma or extreme cavitation bubbles: There is currently no cavity in the wake, and trying to solve a problem that would need to extend there will not converge.
[1] Kermeen, Robert W. "Water tunnel tests of NACA 66_1-012) hydrofoil in noncavitating and cavitating flows." (1956).
%matplotlib inline
%config InlineBackend.figure_format = 'svg'
import os
os.environ["MKL_NUM_THREADS"] = "1"
import os.path
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
import pandas as pd
# Read data into a pandas DataFrame
df = pd.read_csv('tunneldata.csv')
df=df.sort_values(by="K")
online_plot = False
if online_plot:
fig,ax=plt.subplots(2,2)
lines = [None]
# This is a callback function to show relevant properties of the cavitation solver.
# You can define other callbacks if you wish to do so
def plotfun(fig,ax,lines,problem):
p = problem.p
bl = problem.bl
u_max = np.sqrt(1+problem.sigma)
lines[0] = vfp.plot_geometry(ax[0,0],p,bl,lines[0])
NF = len(p.foils)
off = 0
for k in range(NF):
N = p.foils[k].N
index = 0
index += 1
xval = p.foils[k].X[0,:]
yval = 1-np.power(p.gamma_viscid[off:off+N],2)
lines = add_plot(lines,index,xval,yval,ax[0,1],k,NF)
index += 1
yval = 1-np.power(u_max-problem.gam_var[off:off+N],2)
lines = add_plot(lines,index,xval,yval,ax[0,1],k,NF)
index += 1
gambl = np.r_[bl[0].bl_fl.nodes.ue,bl[0].bl_wk.nodes.ue]
for kk in range(2,len(p.foils)):
gambl = np.r_[gambl,bl[kk].bl_fl.nodes.ue,bl[kk].bl_wk.nodes.ue]
yval = 1-np.power(gambl[off:off+N],2)
lines = add_plot(lines,index,xval,yval,ax[0,1],k,NF)
ax[0,1].yaxis.set_inverted(True)
ax[0,1].set_title("CP")
sti = problem.bl[k].sti
legendstr=[]
index += 1
xval = p.foils[k].X[0,0:sti]
yval = problem.vd[off:off+sti]
lines = add_plot(lines,index,xval,yval,ax[1,0],k,NF)
legendstr.append("Cavity Up")
index += 1
yval = problem.bl[k].bl_fl.nodes.delta[off:off+sti]
lines = add_plot(lines,index,xval,yval,ax[1,0],k,NF,style='--')
legendstr.append("BL Up")
index += 1
xval = p.foils[k].X[0,sti:N]
yval = problem.vd[off+sti:off+N]
lines = add_plot(lines,index,xval,yval,ax[1,0],k,NF)
legendstr.append("Cavity Lo")
index += 1
yval = problem.bl[k].bl_fl.nodes.delta[off+sti:off+N]
lines = add_plot(lines,index,xval,yval,ax[1,0],k,NF)
legendstr.append("BL Lo")
if k==0:
ax[1,0].legend(legendstr)
off += N
for key, vec in problem.itr.items():
index += 1
xval = np.arange(0,len(vec))
yval = vec
lines = add_plot(lines,index,xval,yval,ax[1,1],0,NF,True)
ax[0,0].set_title("AOA=%g, u_max=%g, CD=%g)"%(p.alpha,u_max,bl[0].CD+problem.cavdrag))
ax[1,1].legend(problem.itr.keys())
fig.canvas.draw()
return
def add_plot(lines,index,xval,yval,ax,k,NF,semilogy=False, relim=True,style="-"):
if lines is None:
lines = [None]
for j in range(len(lines)-1,index+1):
lines.append(None)
if lines[index] is None:
lines[index] = [None]*NF
if lines[index][k] is None:
if semilogy:
lines[index][k] = ax.semilogy(xval,yval,'.-')
else:
lines[index][k] = ax.plot(xval,yval,style)
else:
lines[index][k][0].set_xdata(xval)
lines[index][k][0].set_ydata(yval)
if relim:
ax.autoscale()
ax.relim()
return lines
N=100
AF = vft.repanel(vft.read_selig("naca_661_012.dat"),N)
# callback = lambda problem: plotfun(fig,ax,lines,problem) # Uncomment to update plot each iteration
callback = None
pol = []
weightList = [.5,1,2,4,8,16] # This list first tries .5, and if unsuccessful the other ones
for angle, group in df.groupby('a (deg)'):
if not (angle in [2,3,4,5,6,7]): continue
RE = 1.1e6*np.max(group['V (fps)'])/40
s = vf.setup(Silent = True,Re = RE,Ma=0,Ncrit=0,Alpha=angle)
s.CmRefX = .25
# Find suitable AOA (because this is not 1-1 match, but we should have same noncavitating lift)
# The lift is given as the lift at largest K/Sigma
s.CLTarget = np.asarray(group['CL'])[np.asarray(np.argmax(group['K']))]
(p,bl,x) = vf.init(AF,s)
vf.iter(x,bl,p,s)
s.Alpha = p.alpha
print("Starting cavitation calculation with non-cavitating CL=%g at AOA=%g for Experimental angle %g"%(s.CLTarget,s.Alpha,angle))
s.CLTarget = np.inf # Disable Lift finding for the following cavitation solver
cd = []
cm = []
cd_cav = []
cd_cf = []
cl = []
sigv = []
flags = []
cd_ref = []
cm_ref = []
cd_cav_ref = []
cl_ref = []
sigv_ref = []
fgfv = []
ilen = None
# Make more points as there are experimantal values in between
SigExp = np.sort(np.asarray(group['K']))
SigNum = np.r_[np.linspace(.3,1,20), np.linspace(1.05,3,5)] # The interesting slope is in the beginning
flag = False
flag_cav = True
for sig in SigNum:
# Ignore if we are already in the non-cavitating zone
if flag_cav:
#print("Solving:%g° - sig=%g"%(s.Alpha,sig))
problem = vf.cav_init(AF,s)
for cavweight in weightList:
flag = vf.cav_solve(problem,problem.s,sig,nsubiter=300,callback=callback,tol=1e-4,restart=1,cavWeight=cavweight)
if flag:
#print("cavweight: %g Drag: %g Lift: %g"%(cavweight,problem.bl[0].CD+problem.cavdrag,problem.p.CL))
print(".",end="")
if online_plot:
plotfun(fig,ax,lines,problem)
# Check if no cavitating anymore
if np.max(problem.vd)==0.0:
flag_cav = False
break
if not flag:
print("x",end="")
cd.append(problem.bl[0].CD)
cd_cav.append(problem.cavdrag)
cl.append(problem.p.CL)
cm.append(problem.p.CM)
fgfv.append(cavweight)
sigv.append(sig)
flags.append(flag)
if online_plot:
plotfun(fig,ax,lines,problem)
print("")
pol.append({'angle': angle, 'cavweight':fgfv, 'alpha': p.alpha, 'sig': sigv, 'Cl': np.asarray(cl),'Cm': np.asarray(cm),
'Cd':np.asarray(cd),'CDCAV':np.asarray(cd_cav),'cp':np.asarray(problem.p.cp.copy()),'flag':np.asarray(flags)})
Starting cavitation calculation with non-cavitating CL=0.193 at AOA=1.76407 for Experimental angle 2 ............. Starting cavitation calculation with non-cavitating CL=0.267 at AOA=2.44376 for Experimental angle 3 ...................... Starting cavitation calculation with non-cavitating CL=0.353 at AOA=3.23832 for Experimental angle 4 ........x.............. Starting cavitation calculation with non-cavitating CL=0.437 at AOA=4.02069 for Experimental angle 5 ........................ Starting cavitation calculation with non-cavitating CL=0.522 at AOA=4.82128 for Experimental angle 6 x........................ Starting cavitation calculation with non-cavitating CL=0.602 at AOA=5.58465 for Experimental angle 7 x........................
matplotlib.rcParams['figure.figsize'] = [12, 7]
if not online_plot:
fig,ax=plt.subplots(2,2)
lines = [None]
# Rerun one (somewhat random) problem
sig = .75
s.Alpha = 2.44662 # that's angle=3 for same non-cavitating CL
for cavweight in weightList:
flag = vf.cav_solve(problem,problem.s,sig,nsubiter=300,callback=callback,tol=1e-4,restart=1,cavWeight=cavweight)
if flag:
plotfun(fig,ax,lines,problem)
break
matplotlib.rcParams['figure.figsize'] = [12, 12]
fig, ax = plt.subplots(3,1)
plot_ref = False
fac_angle_plot = .5 # How far to shift each plot by aoa - Just to make it easier to see how different lines actually follow. Otherwise they collapse.
if not plot_ref:
for i in range(len(pol)):
I = np.nonzero(pol[i]['flag'])[0]
ax[0].plot(np.asarray(pol[i]['sig'])[I]+np.asarray(pol[i]['angle'])*fac_angle_plot,pol[i]['Cd'][I]+pol[i]['CDCAV'][I],'-')
for i in range(len(pol)):
I = np.nonzero(pol[i]['flag'])[0]
ax[1].plot(np.asarray(pol[i]['sig'])[I]+np.asarray(pol[i]['angle'])*fac_angle_plot,pol[i]['Cl'][I],'-')
for i in range(len(pol)):
I = np.nonzero(pol[i]['flag'])[0]
ax[2].plot(np.asarray(pol[i]['sig'])[I]+np.asarray(pol[i]['angle'])*fac_angle_plot,pol[i]['Cm'][I],'-')
ax[0].set_prop_cycle(None)
ax[1].set_prop_cycle(None)
ax[2].set_prop_cycle(None)
else:
for i in range(len(pol)):
ax[0].plot(np.asarray(pol[i]['sigRef'])+np.asarray(pol[i]['angle'])*fac_angle_plot,pol[i]['CdRef']+pol[i]['CDCAVRef'],'-')
for i in range(len(pol)):
ax[1].plot(np.asarray(pol[i]['sigRef'])+np.asarray(pol[i]['angle'])*fac_angle_plot,pol[i]['ClRef'],'-')
for i in range(len(pol)):
ax[2].plot(np.asarray(pol[i]['sigRef'])+np.asarray(pol[i]['angle'])*fac_angle_plot,pol[i]['CmRef'],'-')
ax[0].set_prop_cycle(None)
ax[1].set_prop_cycle(None)
ax[2].set_prop_cycle(None)
ax[0].set_ylabel('CD')
ax[1].set_ylabel('CL')
ax[2].set_ylabel('Cm')
for k in range(3):
ax[k].grid(1)
ax[-1].set_xlabel(r'$\sigma$ + AOA(°)*%1.1g'%fac_angle_plot)
ax[0].set_prop_cycle(None)
ax[1].set_prop_cycle(None)
# Group by 'a (deg)' and plot 'CD' vs 'K' for each group
for angle, group in df.groupby('a (deg)'):
if angle in range(2,8):
ax[0].plot(group['K']+angle*fac_angle_plot, group['CD'], label=f'a = {angle} deg', marker='d', linestyle='',fillstyle='none')
ax[1].plot(group['K']+angle*fac_angle_plot, group['CL'], label=f'a = {angle} deg', marker='d', linestyle='',fillstyle='none')
ax[2].plot(group['K']+angle*fac_angle_plot, group['CM'], label=f'a = {angle} deg', marker='d', linestyle='',fillstyle='none')
