Viscous Inverse Design¶

This notebook demonstrates the use of gradients from viiflow for fully viscous inverse design. It defines a target pressure distribution from one airfoil and, coming from another airfoil, tries to find the shape necessary to arrive at this target pressure. It uses virtual displacements, which do not necessitate the recalculation of the panel operator. Instead, it uses the same model used for the effect of boundary layer thickness onto the flow for modification of the airfoil shape.

The heart of this notebook is a Gauss-Newton iteration which solves for these virtual displacements. Instead of trying to solve the pressure distribution exactly, the iteration sovles a least-squares problem that joins the pressure difference with regularizing terms. Fully viscous inverse design is not a straightforward problem. There are several ways an optimizer may cheat, for example

  • The velocity is defined by the inviscid solution of the airfoil shape plus boundary layer thickness. An optimizer can therefore choose to reduce the thickness of the airfoil if for some reason a thick boundary layer leads to the target velocity distribution.
  • Kinks in the desired velocity are, in the case below, due to laminar-turbulent transition. However, an optimizer can choose to model this kink by an actual kink in the airfoil.

To alleviate this, the pressure error is appended by a regularizing term that penalizes non-smooth displacements - simply by adding $ \frac{\mathrm{d}^2}{\mathrm{d}^2 s} \delta_{virtual}(s) $ at every point along the airfoil surface coordinate $s$ to the Least-Squares problem.

The parameters chosen to increrase/decrease the penalties were chosen ad-hoc by trial and error. In addition, the nodes very close to the stagnation point are not modified.

In addition, the residual $r$ of the viiflow solver itself is added to the Least-Squares problem and scaled such that at convergence its error is sufficiently low. Every iteration then performs for displacements $y$ and the viiflow variables $x$ $$ \begin{aligned} y^{k+1} = y^k - \lambda {\Delta y}^k \\ x^{k+1} = x^k - \lambda {\Delta x}^k \\ {\Delta y}^k, {\Delta x}^k = \min_{\Delta y,\Delta x} \| F(y^k,x^k) - \frac{\partial F}{\partial y}(y^k,x^k) \Delta y - \frac{\partial F}{\partial x}(y^k,x^k) \Delta x\|^2\\ \|F(y,x)\|^2 = \gamma_{cp}^2\|ue(y)-ue_{target}\|^2 + \gamma_y^2\| \frac{\mathrm{d}^2}{\mathrm{d}^2 s} y \|^2 + \gamma_r^2 \|r(y,x)\|^2 \end{aligned} $$ This may seem like a large problem, but the effort for solving the overdetermined least-squares problem grows largely with the degrees of freedom, not the amount of equations.

Below, this procedure is used to morph the S805 airfoil into the S825 airfoil. Even with the regularizing terms, little dips that enforce the laminar-turbulent transition can still be seen when zooming in.

While this solves for an airfoil shape of a specified pressure distribution, it is probably not a very smart idea to use this for actual design. A better idea is to use first an inviscid inverse design method, e.g. conformal mapping [1, 2], and remove the discrepancies using a fully viscid iteration. The benefit of this Gauss-Newton approach is how straightforward additional constraints can be included, e.g. only fit the suction side from .1c onwards or fit multiple target distributions at multiple angles of attack.

In [1]:
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
In [29]:
# Analysis Settings
RE = 1e6
ncrit = 6
Mach = 0.0
alpha = 2.0

N = 200

# Read Airfoils
BASE = vft.repanel(vft.read_selig("S805.dat"),N,KAPFAC=2)
TARGET = vft.repanel(vft.read_selig("S825.dat"),N,KAPFAC=2)
# Solve target for our target cp (or more precisely edge velocity)
s = vf.setup(Re=RE,Ma=Mach,Ncrit=ncrit,Alpha=alpha,Silent=False)

# Internal iterations
s.Itermax = 100

# Set-up and initialize based on inviscid panel solution
[p,bl,x] = vf.init([TARGET],s)
res = None
grad = None



# Solve aerodynamic problem of target airfoil
vf.iter(x,bl,p,s)
XT0 = p.foils[0].X[0,:].copy()
UT = p.gamma_viscid[0:p.foils[0].N].copy()

# Set-up and initialize based on inviscid panel solution
[p,bl,x0] = vf.init([BASE],s)
res = None
grad = None

# Solve aerodynamic problem of current airfoil and save for later plotting
[x0,_,res,grad,_] = vf.iter(x0,bl,p,s,None,None)
XC0 = p.foils[0].X[0,:].copy()
UC = p.gamma_viscid[0:p.foils[0].N].copy()

# To interpolate from one grid to the next, suction and pressure side must have unique grid points
# That is why below a grid is created where the pressure side is appended with *-1 at the nose
XT = XT0.copy()
XC = XC0.copy()
XT[np.argmin(XT0)+1::] = 2*XT0[np.argmin(XT0)]-XT0[np.argmin(XT0)+1::]
XC[np.argmin(XC0)+1::] = 2*XC0[np.argmin(XC0)]-XC0[np.argmin(XC0)+1::]

# Interpolate target pressure onto current airfoil grid
UT = np.interp(-XC.flatten(),-XT.flatten(),np.asarray(UT[:,0]).flatten())

Iteration 17, |res| 0.000034, lam 1.000000
Iteration 14, |res| 0.000030, lam 1.000000
In [57]:
# Weighting factors for Gauss-Newton
facx = 10 # Penalty for smooth dioscplacement
fac_err = .5 #Weighting of cp error w.r.t. above penalties
fac_res = 1e4
s.Gradients = True

NAERO = x.shape[0]
NVD = len(XC)

# Set-up and initialize based on inviscid panel solution
[p,bl,x0] = vf.init([BASE],s)
res = None
grad = None

# Solve aerodynamic problem to convergence
[x,_,_,_,_] = vf.iter(x0,bl,p,s,None,None)
fprev = np.inf
# Find ST and do not change near there
II = np.logical_and(np.fabs(XT-XT[bl[0].sti])>0.001,p.foils[0].X[0,:].ravel()>np.amin(p.foils[0].X[0,:].ravel()))
II[0]=False
II[NVD-1]=False
iter = 0
lam = 1.0
y = np.zeros(NVD)
while True:
    iter+=1

    # Solve Aerodynamic problem
    s.Itermax = 0
    s.Silent = True
    [_,_,res,grad,gradients] = vf.iter(x,bl,p,s,None,None,y)

    # Residual 
    RESy = fac_err*(p.gamma_viscid[0:p.foils[0].N,0]-UT)
    dRESydy = fac_err*gradients.partial.gam_vd[0:NVD,:]
    dRESydx = fac_err*gradients.partial.gam_x[0:NVD,:]
    
    # Penalty for smooth displacement
    difforder = 2
    REGdelta = np.diff(y,difforder)*facx
    dREGdeltady = np.diff(np.eye(NVD),difforder,0)*facx
    dREGdeltadx = np.zeros((len(REGdelta),len(x)))
    
    # Gauss-Newton step from all terms
    F = np.r_[RESy,REGdelta,res*fac_res]
    fcurr = np.sqrt(F.T@F)
    
    y0 = y
    fprev = fcurr
    # Find ST and do not change near there
    II = np.logical_and(np.fabs(XT-XT[bl[0].sti])>0.001,p.foils[0].X[0,:].ravel()>np.amin(p.foils[0].X[0,:].ravel()))
    II[0]=False
    II[NVD-1]=False
    
    dFdy = np.r_[dRESydy,dREGdeltady,gradients.partial.res_vd*fac_res]
    dFdx = np.r_[dRESydx,dREGdeltadx,grad*fac_res]
    dF = np.c_[dFdy[:,II],dFdx]
    dX = -np.linalg.lstsq(dF,F,rcond=None)[0]
    dy = dX[0:np.sum(II)]
    dx = dX[np.sum(II)::]
    lam = 1
    
    # Print
    resaero = np.sqrt(np.matmul(res,res.T))

    
    # Ad-hoc Damping
    for k in range(len(dy)):
        lam = np.fmin(lam,0.003/abs(dy[k])) # Do not move virtual displacement more than 1mm
    for k in range(len(x)):
        lam = np.fmin(lam,.2/(abs(dx[k]/x[k])))
        
    print("iter %u res p:%f resaero: %f dvd:%f lam:%f"%(iter, np.sqrt(np.matmul(F,F.T)), \
                resaero,np.sqrt(np.matmul(dy,dy.T)),lam))
        
    if np.sqrt(np.matmul(dy,dy.T))<1e-4 and resaero<1e-4:
        print("Converged")
        break
        
    if iter>200:
        print("Not Converged (iteration)")
        break
        
    j =0
    for k in np.argwhere(II):
        y[k] += lam*dy[j]
        j+=1
    x += lam*dx

iter 1 res p:11305.033001 resaero: 1.130503 dvd:0.228020 lam:0.080688
iter 2 res p:10328.224175 resaero: 1.032822 dvd:0.209047 lam:0.084548
iter 3 res p:9438.900520 resaero: 0.943890 dvd:0.194219 lam:0.078743
iter 4 res p:8704.320998 resaero: 0.870432 dvd:0.175230 lam:0.069111
iter 5 res p:7908.335966 resaero: 0.790834 dvd:0.135116 lam:0.039980
iter 6 res p:7632.007492 resaero: 0.763201 dvd:0.125562 lam:0.035724
iter 7 res p:7386.510312 resaero: 0.738651 dvd:0.123468 lam:0.049718
iter 8 res p:7050.623860 resaero: 0.705062 dvd:0.121002 lam:0.046532
iter 9 res p:6759.003911 resaero: 0.675900 dvd:0.119752 lam:0.041989
iter 10 res p:6505.956603 resaero: 0.650596 dvd:0.116854 lam:0.035197
iter 11 res p:6290.880531 resaero: 0.629088 dvd:0.116294 lam:0.048219
iter 12 res p:6014.220121 resaero: 0.601422 dvd:0.114600 lam:0.050947
iter 13 res p:5740.229215 resaero: 0.574023 dvd:0.112556 lam:0.051888
iter 14 res p:5475.621722 resaero: 0.547562 dvd:0.109549 lam:0.053041
iter 15 res p:5222.029939 resaero: 0.522203 dvd:0.107223 lam:0.054684
iter 16 res p:4986.220836 resaero: 0.498622 dvd:0.104677 lam:0.056858
iter 17 res p:4766.314868 resaero: 0.476631 dvd:0.102026 lam:0.059676
iter 18 res p:4559.556939 resaero: 0.455956 dvd:0.107796 lam:0.031394
iter 19 res p:4733.323402 resaero: 0.473332 dvd:0.098253 lam:0.062638
iter 20 res p:4508.371640 resaero: 0.450837 dvd:0.095269 lam:0.078171
iter 21 res p:4248.131321 resaero: 0.424813 dvd:0.095056 lam:0.041807
iter 22 res p:4322.146994 resaero: 0.432215 dvd:0.088141 lam:0.081787
iter 23 res p:4067.778364 resaero: 0.406778 dvd:0.085139 lam:0.066255
iter 24 res p:3988.273711 resaero: 0.398827 dvd:0.080102 lam:0.135227
iter 25 res p:3566.353301 resaero: 0.356635 dvd:0.071289 lam:0.185355
iter 26 res p:3136.901848 resaero: 0.313690 dvd:0.060076 lam:0.169985
iter 27 res p:6465.571023 resaero: 0.646557 dvd:0.051849 lam:0.223166
iter 28 res p:5200.625916 resaero: 0.520063 dvd:0.041839 lam:0.152437
iter 29 res p:4526.469867 resaero: 0.452647 dvd:0.036726 lam:0.299740
iter 30 res p:3513.462762 resaero: 0.351346 dvd:0.026821 lam:0.142870
iter 31 res p:3146.002769 resaero: 0.314600 dvd:0.023360 lam:0.397235
iter 32 res p:2351.200052 resaero: 0.235120 dvd:0.016212 lam:0.419092
iter 33 res p:1578.624575 resaero: 0.157862 dvd:0.010331 lam:0.460962
iter 34 res p:1085.743172 resaero: 0.108574 dvd:0.007741 lam:0.221209
iter 35 res p:948.818562 resaero: 0.094882 dvd:0.004912 lam:0.582935
iter 36 res p:488.285292 resaero: 0.048829 dvd:0.006175 lam:0.115428
iter 37 res p:531.779479 resaero: 0.053178 dvd:0.002123 lam:0.722645
iter 38 res p:739.068449 resaero: 0.073907 dvd:0.006015 lam:0.127418
iter 39 res p:696.198730 resaero: 0.069620 dvd:0.002016 lam:0.774806
iter 40 res p:456.333183 resaero: 0.045633 dvd:0.000894 lam:0.972758
iter 41 res p:194.843761 resaero: 0.019484 dvd:0.003830 lam:0.177509
iter 42 res p:266.831234 resaero: 0.026683 dvd:0.001045 lam:0.621735
iter 43 res p:915.206836 resaero: 0.091521 dvd:0.001251 lam:0.711178
iter 44 res p:308.416426 resaero: 0.030842 dvd:0.001273 lam:0.559955
iter 45 res p:214.921865 resaero: 0.021492 dvd:0.001425 lam:0.412011
iter 46 res p:214.957751 resaero: 0.021496 dvd:0.001019 lam:0.659308
iter 47 res p:193.617923 resaero: 0.019362 dvd:0.002261 lam:0.253245
iter 48 res p:243.532997 resaero: 0.024353 dvd:0.000874 lam:0.855105
iter 49 res p:263.166584 resaero: 0.026317 dvd:0.003293 lam:0.185278
iter 50 res p:287.202603 resaero: 0.028720 dvd:0.000847 lam:1.000000
iter 51 res p:117.143286 resaero: 0.011714 dvd:0.002973 lam:0.214944
iter 52 res p:153.667825 resaero: 0.015367 dvd:0.000794 lam:0.815472
iter 53 res p:154.316409 resaero: 0.015432 dvd:0.003069 lam:0.195672
iter 54 res p:217.396835 resaero: 0.021740 dvd:0.000588 lam:1.000000
iter 55 res p:84.145800 resaero: 0.008415 dvd:0.003242 lam:0.196991
iter 56 res p:142.383891 resaero: 0.014238 dvd:0.000692 lam:0.945732
iter 57 res p:143.936712 resaero: 0.014394 dvd:0.003239 lam:0.186751
iter 58 res p:210.651348 resaero: 0.021065 dvd:0.000497 lam:1.000000
iter 59 res p:49.909226 resaero: 0.004991 dvd:0.002352 lam:0.365554
iter 60 res p:204.289179 resaero: 0.020429 dvd:0.000919 lam:1.000000
iter 61 res p:319.960958 resaero: 0.031996 dvd:0.000371 lam:1.000000
iter 62 res p:81.425119 resaero: 0.008143 dvd:0.000402 lam:1.000000
iter 63 res p:60.072560 resaero: 0.006007 dvd:0.000236 lam:1.000000
iter 64 res p:29.166717 resaero: 0.002917 dvd:0.000188 lam:1.000000
iter 65 res p:21.104512 resaero: 0.002110 dvd:0.000117 lam:1.000000
iter 66 res p:12.624065 resaero: 0.001262 dvd:0.000088 lam:1.000000
iter 67 res p:8.819571 resaero: 0.000882 dvd:0.000051 lam:1.000000
iter 68 res p:6.003555 resaero: 0.000600 dvd:0.000023 lam:1.000000
iter 69 res p:1.771110 resaero: 0.000177 dvd:0.000035 lam:1.000000
iter 70 res p:3.436607 resaero: 0.000344 dvd:0.000010 lam:1.000000
iter 71 res p:0.822305 resaero: 0.000082 dvd:0.000024 lam:1.000000
Converged
In [58]:
%matplotlib inline
%config InlineBackend.figure_format = 'svg'
matplotlib.rcParams['figure.figsize'] = [11, 5.5]
fig,ax = plt.subplots(1,1)
ax.plot(p.foils[0].X[0,:],np.power(UC,2)-1,'-k')
ax.plot(p.foils[0].X[0,:],np.power(p.gamma_viscid[0:p.foils[0].N],2)-1,'-',color=(0.6,0.6,0.6))
ax.plot(p.foils[0].X[0,:],np.power(UT,2)-1,'2k')
ax.legend(['Initial Pressure','Found Pressure','Target Pressure'])
xlim = ax.get_xlim()

fig,ax = plt.subplots(1,1)
lines = None
ax.plot(TARGET[0,:],TARGET[1,:],'2k')
lines = vfp.plot_geometry(ax,p,bl,lines)
ax.legend(['Target Airfoil','Initial Geometry','Found Geometry'])
ax.set_xlim(xlim)
Out[58]:
(-0.049969536899920834, 1.0499985493761868)
No description has been provided for this image
No description has been provided for this image

[1] Selig, Michael S., and Mark D. Maughmer. Generalized multipoint inverse airfoil design. AIAA journal 30.11 (1992): 2618-2625.

[2] Drela, Mark. XFOIL: An analysis and design system for low Reynolds number airfoils. Low Reynolds number aerodynamics. Springer, Berlin, Heidelberg, 1989. 1-12.