Animating streamlines of the flow past a rotating Joukowski airfoil

The streamlines of the flow past a Joukowski airfoil are generated following this chapter from the Internet Book of Fluid Dynamics.

Visualization of streamlines is based on the property of the complex flow with respect to a conformal transformation:

If w is the complex plane of the airfoil, z is the complex plane of the circle, as the section in a circular cylinder, and $w=w(z)$ is a conformal tranformatiom mapping the circle to the airfoil, then the complex flow, $F$, past the airfoil is related to the complex flow, $f$, past the circle(cylinder) by: $F(w)=f(z(w))$ or equivalently $F(w(z))=f(z)$.

The streamlines of each flow are defined as contour plots of the imaginary part of the complex flow. In our case, due to the latter relation, we plot the contours of $Imag{(f)}$ over $w(z)$, where $w(z)$ is the Joukowski transformation, that maps a suitable circle onto the airfoil.

Algorithm to generate a frame:

m- Choose the center of the circle as in the Internet Book, but with respect to a new coordinate $\tilde{z}$, $\tilde{z}=e^{i \alpha}z$, i.e $\tilde{z}_c=\lambda-R e^{i\beta}$ .

  • The Joukowski transformation $\tilde{J}(\tilde{z})=\tilde{z}+\lambda^2/\tilde{z}=\tilde{w}$ maps this circle to an airfoil referenced to the coordinate $\tilde{w}=e^{i \alpha}w$.

  • The circle center has the coordinate $z_c=e^{-i\alpha}\tilde{z}_c$, in the z-complex plane.

  • The Joukowski transformation from the complex plane z, to w, is conjugated via a rotation with the Joukowski tranformation, $\tilde{J}$, i.e. $J(z)= e^{-i \alpha} \underbrace{\tilde{J}(e^{i \alpha}z)}_{\tilde{w}}$

Hence, $J(z)=z+\lambda^2e^{-2 i\alpha}/z$

Each animation frame displays the flow streamlines of angle of attack 0 with respect to the z-coordinate. Due to the inclination of the airfoil the corresponding streamlines are in fact the streamlines of the flow expressed in the $\tilde{z}$-plane, with angle of attack $\alpha$.


In [3]:
import plotly.graph_objects as go
import numpy as np
import numpy.ma as ma
import matplotlib.pyplot as plt

In [4]:
def J(z, lam, alpha):
    return z+(np.exp(-1j*2*alpha)*lam**2)/z

In [5]:
def circle(C, R):
    t = np.linspace(0,2*np.pi, 200)
    return C + R*np.exp(1j*t)

In [6]:
def deg2radians(deg):
    return deg*np.pi/180

In [7]:
def flowlines(alpha=10, beta=5, V_inf=1, R=1, ratio=1.2):
    #alpha, beta are given in degrees
    #ratio =R/lam
    alpha = deg2radians(alpha)# angle of attack
    beta = deg2radians(beta)# -beta is the argument of the complex no (Joukowski parameter - circle center)
    if ratio <= 1: #R/lam must be >1
        raise ValueError('R/lambda must be > 1')
    lam = R/ratio#lam is the parameter of the Joukowski transformation
   
    center_c = np.exp(-1j*alpha)*(lam-R*np.exp(-1j*beta))

    Circle = circle(center_c,R)
    Airfoil = J(Circle, lam, alpha)
    X = np.arange(-3, 3, 0.1)
    Y = np.arange(-3, 3, 0.1)

    x,y = np.meshgrid(X, Y)
    z = x+1j*y
    z = ma.masked_where(np.absolute(z-center_c)<=R, z)
    w = J(z, lam, alpha)
    beta = beta+alpha
    Z = z-center_c

    Gamma = -4*np.pi*V_inf*R*np.sin(beta)#circulation
    U = np.zeros(Z.shape, dtype=np.complex)
    with np.errstate(divide='ignore'):#
        for m in range(Z.shape[0]):
            for n in range(Z.shape[1]):# due to this numpy bug https://github.com/numpy/numpy/issues/8516
                                       #we evaluate  this term of the flow elementwise
                U[m,n] = Gamma*np.log((Z[m,n])/R)/(2*np.pi)
    c_flow = V_inf*Z + (V_inf*R**2)/Z - 1j*U #the complex flow 
   
    
    return w, c_flow.imag, Airfoil

In [8]:
def get_contours(mplcont):
    conts = mplcont.allsegs  #  get the segments of line computed via plt.contour
    xline = []
    yline = []

    for  cont in conts:
        if len(cont) != 0:
            for arr in cont: 
                
                xline += arr[:,0].tolist()
                yline += arr[:,1].tolist()
                xline.append(None) 
                yline.append(None)

    return xline, yline

In [16]:
levels = np.arange(-3, 3.7, 0.25).tolist()
plt.figure(figsize=(0.05,0.05))
plt.axis('off')
Alpha=   list(range(0, 19)) + list(range(17, -19, -1)) + list(range(-17, 1))

frames = []

for k, alpha in enumerate(Alpha):
    Jz, stream_func, Airfoil = flowlines(alpha=alpha)
    #Define an instance of the mpl class contour
    cp = plt.contour(Jz.real, Jz.imag, stream_func, levels=levels, colors='blue')
   
    xline, yline = get_contours(cp)

    frames.append( go.Frame(data=[go.Scatter(x=xline, 
                                             y=yline),
                                  go.Scatter(x=Airfoil.real, 
                                             y=Airfoil.imag),
                                 ],
                         traces=[0, 1],
                         name=f'frame{k}'  
                        ) )



In [10]:
data = [go.Scatter( 
           x=frames[0].data[0].x,
           y=frames[0].data[0].y,
           mode='lines',
           line=dict(color='blue', width=1),
           name=''
          ),
        go.Scatter(
           x=frames[0].data[1].x,
           y=frames[0].data[1].y,
           mode='lines',
           line=dict(color='blue', width=2),
           name=''
          )
      ]

In [11]:
def get_sliders(Alpha, n_frames, fr_duration=100, x_pos=0.0, y_pos=0, slider_len=1.0):
    # n_frames= number of frames
    #fr_duration=duration in milliseconds of each frame
    #x_pos x-coordinate where the slider starts
    #slider_len is a number in (0,1] giving the slider length as a fraction of x-axis length 
    return [dict(steps= [dict(method= 'animate',#Sets the Plotly method to be called when the
                                                #slider value is changed.
                              args= [ [ f'frame{k}'],#Sets the arguments values to be passed to 
                                                              #the Plotly method set in method on slide
                                      dict(mode= 'immediate',
                                           frame= dict( duration=fr_duration, redraw= False ),
                                           transition=dict( duration= 0)
                                          )
                                    ],
                              label=str(alpha)
                             ) for k, alpha in enumerate(Alpha)], 
                transition= dict(duration= 0 ),
                x=x_pos,
                y=y_pos, 
                currentvalue=dict(font=dict(size=12), 
                                  prefix="Angle of attack: ", 
                                  visible=True, 
                                  xanchor= "center"
                                 ),  
                len=slider_len)
           ]

In [12]:
axis = dict(showline=True, zeroline=False, ticklen=4, mirror=True, showgrid=False)


layout = go.Layout(title_text="Streamlines of a flow past a rotating Joukowski airfoil",
                   title_x=0.5,
            font=dict(family='Balto'),
            showlegend=False, 
            autosize=False, 
            width=600, 
            height=600, 
            xaxis=dict(axis, **{'range': [ma.min(Jz.real), ma.max(Jz.real)]}),
            yaxis=dict(axis, **{'range':[ma.min(Jz.imag), ma.max(Jz.imag)]}),
            
            plot_bgcolor='#c1e3ff',
            hovermode='closest',
            
            updatemenus=[dict(type='buttons', showactive=False,
                                y=1,
                                x=1.15,
                                xanchor='right',
                                yanchor='top',
                                pad=dict(t=0, l=10),
                                buttons=[dict(
                                            label='Play',
                                            method='animate',
                                            args=[None, dict(frame=dict(duration=50, redraw=False), 
                                                             transition=dict(duration=0),
                                                             fromcurrent=True,
                                                             mode='immediate')])]
                               )
                          ])

layout.update(sliders=get_sliders(Alpha, len(frames), fr_duration=50))
fig = go.Figure(data=data,layout=layout, frames=frames)

In [13]:
import chart_studio.plotly as py
py.iplot(fig, filename='streamlJouk1485353936.44' )


Out[13]:

In [14]:
from IPython.core.display import HTML
def  css_styling():
    styles = open("./custom.css", "r").read()
    return HTML(styles)
css_styling()


Out[14]: