This example shows how to use viflow for fluid-structure interaction. The airfoil as the famous NACA0012, and we assume that it somehow behaves like a beam that is fixed in the front part up to .2c. To to this, we will not rely on viiflow to solve the aerodynamic problem, but just use its computed gradient and residual to couple those with our structure problem for a joint Newton step.
This necessitates:
This is why this notebook is rather long.
%matplotlib inline
%config InlineBackend.figure_format = 'svg'
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
# Settings
RE = 1e6
ncrit = 9.0
Mach = 0.0
alpha = 8.0
# Read Airfoil Data
NACA0012 = vft.repanel(vft.read_selig("NACA0012STE.dat"),200)
Let $s \in [0,L]$ be the variable along the length of the beam and $q$ be the load on this beam at these elements. The simple beam fixed at the s(0) is defined by $$ \frac{d^2}{ds^2} EI \frac{d^2}{ds^2} y = q $$ $$ y(0) = 0, \frac{d}{ds} y(0) = 0 $$ $$ \frac{d^2}{ds^2} y(L) = 0, \frac{d^3}{ds^3} y(L) = 0 $$
To keep things sane, we will use $M = - EI \frac{d^2}{ds^2} y$ as an additional variable.
$$ - \frac{d^2}{ds^2} M = q, M(L) = 0, \frac{d}{ds} M(L) = 0 \ (I) $$$$ - EI \frac{d^2}{ds^2} y = M (I) , y(0) = 0, \frac{d}{ds} y(0) = 0 \ (II) $$Since we use the foil coordinates from the top, the actual free end is at $s=0$, while the fixed and is at the end of $s$.
The load $q$ is due to the aerodynamic pressure difference between top and bottom side. From the external velocities $u_{top}$ and $u_{bot}$ on the airfoil sides it is calculated by $q=u_{top}^2-u_{bot}^2$.
# Define Structure Model
def beam(s,q):
# s is the variable along the length of the beam
# q is the load on this beam at these elements
EI = 1 # Random scale. Too low: So much movement that you see transient wiggling. Too high: Very little change in geometry.
N = len(s) # Number of elements
# Second order differential operator
# d²/ds²
D = np.zeros((N,N))
for k in range(1,N-1):
dm = s[k]-s[k-1]
dp = s[k+1]-s[k]
D[k,k-1] = 2/dm/(dm+dp)
D[k,k] = -2/dp/(dm+dp)-2/dm/(dm+dp)
D[k,k+1] = 2/dp/(dm+dp)
# Operator (I)
O1 = -D
O1[0,0] = 1
O1[-1,[0, 1]] = [-1,1]/(s[1]-s[0])
# Operator (II)
O2 = -np.multiply(EI,D)
O2[0,N-1] = 1
O2[-1,[N-2, N-1]] = [-1,1]/(s[N-1]-s[N-2])
# Solve for displacement
E0 = np.eye(N)
E0[0,0] = 0
E0[-1,-1] = 0
M = np.linalg.solve(O1,np.matmul(E0,q))
dMdq = np.linalg.solve(O1,E0)
y = np.linalg.solve(O2,np.matmul(E0,M))
dydq = np.linalg.solve(O2,np.matmul(E0,dMdq))
return [y,dydq]
def af_root(p):
# Find root of airfoil
iroot = 0
for k in range(1,p.foils[0].N):
if p.foils[0].X[0,k]>p.foils[0].X[0,k-1]:
iroot = k
break
return iroot
# How to calculate pressure difference on this beam based on the aerodynamics
def q_beam(SB,p,dgam):
iroot = af_root(p)
N = p.foils[0].N
# interpolate velocities to beam SB
utop = np.interp(-SB,-p.foils[0].X[0,0:iroot],np.asarray(p.gamma_viscid[0:iroot]).flatten())
ubot = np.interp(SB,p.foils[0].X[0,iroot::],np.asarray(p.gamma_viscid[iroot:N]).flatten())
# Pressure difference
q = utop**2-ubot**2
# Gradient. Only true for linear interpolation!
dq = None
if not dgam is None:
dq = np.zeros((q.shape[0],dgam.shape[1]))
# dq = dq_utop*dutop + dq_ubot*dubot
# = 2*utop*dutop + ...
for i in range(dgam.shape[1]):
dq[:,i]+= 2*utop*np.interp(-SB,-p.foils[0].X[0,0:iroot-1],np.asarray(dgam[0:iroot-1,i]).flatten())
dq[:,i]+= 2*ubot*np.interp(SB,p.foils[0].X[0,iroot::],np.asarray(dgam[iroot:N,i]).flatten())
return [q,dq]
First, we solve the aerodynamic problem without the structural part. The solution is used as a starting point for the subsequent calculations.
# Setup for viiflow
s = vf.setup(Re=RE,Ma=Mach,Ncrit=ncrit,Alpha=alpha)
# Internal iterations
#s.VirtualGradients = True # You need to set this for FSI problems.
s.Gradients = True # You need to set this for FSI problems.
# If iterate_wakes is False, the wake does not change when moving which may lead to better convergence.
# If True, the wake changes. Because the gradient does not take movement of nodes into account,
# this may leasd to poorer convergence.
s.IterateWakes = True
# Set-up and initialize based on inviscid panel solution
[p,bl,x] = vf.init([NACA0012],s)
res = None
grad = None
# Solve aerodynamic problem
[x,_,_,_,_] = vf.iter(x,bl,p,s,None,None)
cp0 = p.cp.copy();
# FSI Setup
# This is our beam, the x values of the top side from TE until 0.2
ii = 0 # Index where the beam starts
for k in range(p.foils[0].N):
if p.foils[0].X[0,k]<0.2:
ii = k
break
SB = p.foils[0].X[0,0:ii]
ND = len(SB)
# How do we get the complete foil displacement from xtop?
def vdfull(p,xtop):
NXT = xtop.shape[0]
NP = p.foils[0].N
iroot = af_root(p)
vd = np.zeros(p.foils[0].N)
vd[0:NXT] = xtop
vd[iroot::] = -np.interp(-p.foils[0].X[0,iroot::],-p.foils[0].X[0,0:iroot],vd[0:iroot])
# Gradient
dvd_xtop = np.zeros((NP,NXT))
E = np.eye(NXT)
dvd_xtop[0:NXT,0:NXT] = E
for k in range(NXT):
dvd_xtop[iroot::,k] = -np.interp(-p.foils[0].X[0,iroot::],-p.foils[0].X[0,0:NXT],E[:,k])
return [vd, dvd_xtop]
The next section combines viiflow with the above structural model to solve for stationary conditions. That is, due to the pressure difference $q$ the airfoil is deformed, which changes its pressure variation and hence $q$. In the end, we need a structural displacement $y_{var}$ that is equivalent to the deformation due to the loads due to this deformation $y_q(q(y_{var}))=y_{var}$.
# Newton Iteration
tol = 1e-5
# Additional unknowns: The displacement on the top side
# On the bottom, we have the negative displacement at the same x coordinates
xtop = np.zeros(ii) # Initialize with 0
xj = np.r_[x,xtop] # Joined aerodynamic and structural variables
NAERO = x.shape[0]
NTOP = ii
s.Itermax = 0 # Disable internal iterations. Important!
s.Silent = True
vres_s = [] # Vector of stuctural residuals
vres_a = [] # Vector of aerodynamic residuals
for k in range(50):
res = None
grad = None
# Current virtual displacements
[vd, dvd_xj] = vdfull(p,xj[NAERO::])
# Solve Aerodynamics
[xj[0:NAERO],flag,res,grad,grad_vii] = vf.iter(xj[0:NAERO],bl,p,s,res,grad,vd)
# Gradient of aerodynamic problem w.r.t. xj
dres_xj = np.c_[grad,np.matmul(grad_vii.partial.res_vd,dvd_xj)]
# Get structural loads
[q,dq_dxaero] = q_beam(SB,p,grad_vii.partial.gam_x)
[_,dq_dvd] = q_beam(SB,p,grad_vii.partial.gam_vd)
# Solve structural part
[y,dydq] = beam(SB,q)
# Gradient of structural problem w.r.t. xj
dy_xaero = np.matmul(dydq,dq_dxaero)
dy_xtop = np.matmul(dydq,np.matmul(dq_dvd,dvd_xj))
# Interaction goal: y - xtop = 0
# Aerodynamic goal: res = 0
res_s = y-xj[NAERO::]
F = np.r_[res_s,res]
DF = np.r_[np.c_[dy_xaero,dy_xtop-np.eye(NTOP)],dres_xj]
# Newton Step
dxj = -np.linalg.solve(DF,F)
# Some Newton step-size restriction
lam = 1.0
for k in range(NAERO-3):
if abs(dxj[k])>0:
lam = min(lam,.1/(abs(dxj[k])/abs(xj[k])))
for k in range(NAERO,NAERO+NTOP):
if abs(dxj[k])>0:
lam = min(lam,.1/(abs(dxj[k])/max(abs(xj[k]),0.02)))
xj += lam*dxj
vres_s.append(np.sqrt(np.sum(res_s*res_s)))
vres_a.append(np.sqrt(np.sum(res*res)))
if np.sqrt(np.sum(F*F))<tol or np.sqrt(np.sum(dxj*dxj))<tol:
print("Converged")
break
matplotlib.rcParams['figure.figsize'] = [12, 6]
# Plot current geometry
fig,ax=plt.subplots(1,1)
ax.plot(p.foils[0].X[0,:],-.1*cp0+.2,'-',color='tab:orange')
ax.plot(p.foils[0].X[0,:],-.1*p.cp+.2,'-',color='tab:red')
vfp.plot_geometry(ax,p,bl)
# Plot pressure in same plot
ax.legend(('Initial pressure','Final pressure','Initial Geometry','Deformed Geometry'))
ax.set_yticks([])
# Plot Residuals with iterations
fig, ax = plt.subplots()
ax.semilogy(vres_s)
ax.semilogy(vres_a)
ax.legend(('Structural Residual','Aerodynamic Residual'))
ax.set_xlabel('Iteration');