# Stagnation flow

This notebook will use Dedalus to create a minimal working example of the solution to viscous flow past a stagnation point, and the volume penalized approximation

\begin{align}
(\partial_{x} u)^2 - \partial_x^2 u \, u - \frac{1}{\text{Re}} \partial_x^3 u -  \frac{1}{\text{Re}\,\varepsilon^2}\Gamma u &= 1
\end{align}

where the penalty mask $\Gamma$ can be chosen by shifting and rescaling optimal masks $\Gamma^*$

\begin{align}
\Gamma(x) &= \Gamma^*\left(\frac{x - \varepsilon \ell}{\varepsilon \delta}\right)
\end{align}

and the physical solution $u$ satisfies the boundary conditions

\begin{align}
x &= -1 & u &= 0\\
x &= -1 & \partial_ x u &= 0\\
x &\to +\infty & \partial_x u &\to 1
\end{align}

# Imports


```python
import numpy as np
import matplotlib.pyplot as plt
import dedalus.public as de

from matplotlib import rc
rc('font',**{'family':'serif','serif':['Computer Modern Roman']})
rc('text', usetex=True)
```

# Ideal no-slip solution 

We will first solve the problem numerically using true no-slip boundaries


```python
N = 256
Re = 1

xbasis = de.Chebyshev('x', N, interval=(0,10), dealias=3/2)
domain = de.Domain([xbasis], grid_dtype=np.float64)
x0 = xbasis.grid(*domain.dealias)

stag = de.NLBVP(domain,variables=['u','ux','uxx'])
stag.meta[:]['x']['dirichlet'] = True
stag.parameters['Re'] = Re
stag.add_equation('dx(u) - ux = 0')
stag.add_equation('dx(ux) - uxx = 0')
stag.add_equation('dx(uxx)/Re = -1 - uxx*u + ux*ux')
stag.add_bc('right(ux) = 1')
stag.add_bc('left(u) = 0')
stag.add_bc('left(ux) = 0')

solver = stag.build_solver()
u0, ux0, uxx0 = (solver.state[name] for name in stag.variables)
for field in [u0, ux0, uxx0]: field.set_scales(domain.dealias)

u0['g'], ux0['g'], uxx0['g'] = x0, 1, 0

tolerance = 1e-10
pert = solver.perturbations.data
pert.fill(1+tolerance)

while np.sum(np.abs(pert)) > tolerance:
    solver.newton_iteration()
```

# Volume penalized solution 

We now solve the volume penalized problem.


```python
N = [32,64]
Re = 1
η = 0.01
ϵ = np.sqrt(η/Re)
l = ϵ

xb0 = de.Chebyshev('x0', N[0], interval=(-1,l), dealias=3/2)
xb1 = de.Chebyshev('x1', N[1], interval=(l,5), dealias=3/2)
xbasis = de.Compound('x',[xb0,xb1])
domain = de.Domain([xbasis], grid_dtype=np.float64)
x = xbasis.grid(*domain.dealias)

Γ = domain.new_field(scales=domain.dealias)
Γ['g'] = 1.*(x<l)

stag = de.NLBVP(domain,variables=['u','ux','uxx'])
stag.meta[:]['x']['dirichlet'] = True
stag.parameters['Re'] = Re
stag.parameters['Γ'] = Γ
stag.parameters['ε'] = ϵ
stag.add_equation('dx(u) - ux = 0')
stag.add_equation('dx(ux) - uxx = 0')
stag.add_equation('dx(uxx)/Re - (Γ/(Re*ε**2))*ux = - 1 - uxx*u + ux*ux')
stag.add_bc('right(ux) = 1')
stag.add_bc('left(u) = 0')
stag.add_bc('left(ux) = 0')

solver = stag.build_solver()
u, ux, uxx = (solver.state[name] for name in stag.variables)
for field in [u, ux, uxx]: field.set_scales(domain.dealias)

u['g'], ux['g'], uxx['g'] = x, 1, 0

tolerance = 1e-10
pert = solver.perturbations.data
pert.fill(1+tolerance)

while np.sum(np.abs(pert)) > tolerance:
    solver.newton_iteration()
```


```python
# Comparison of penalized and reference solutions
fig, ax = plt.subplots()
ax.plot(x,u['g'],'C2',label='Penalized')
ax.plot(x0,u0['g'],'k--',label='Reference')
ax.fill_between(x,0,1.5*Γ['g'],color='lightgray')
ax.set(aspect=1,xlim=[-1,2],ylim=[0,1.5],yticks=np.arange(5),
       xlabel='$x$',ylabel='$u$')
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
ax.legend()
```

e372e5e31ed2dfee5d54289ad1cf4bc61d0ae0a9
ipynb
Jupyter Notebook
stagnation-flow.ipynb
ericwhester/volume-penalty-code
66a1745daeae2ad71bda0bc9299c8b8271a9871f
MIT
2020-03-14T19:40:40.000Z
2022-03-18T03:02:33.000Z
stagnation-flow.ipynb
ericwhester/volume-penalty-code
66a1745daeae2ad71bda0bc9299c8b8271a9871f
MIT
stagnation-flow.ipynb
ericwhester/volume-penalty-code
66a1745daeae2ad71bda0bc9299c8b8271a9871f
MIT
Qwen/Qwen-72B
1. YES
2. YES