2. Potential flow with an impenetrable surface

This part introduces the treatment of rigid impenetrable surfaces in GridPotentialFlow. We can impose the no-penetration constraint by setting the fluid streamfunction equal to that of the surface $\mathfrak{s}_b$, up to a uniform value $s_0$. The discrete no-penetration constraint is thus

\[\mathsf{Es}=\mathfrak{s}_{b}-\mathsf{E} \mathsf{s}_{\infty}-\mathfrak{s}_{0},\]

where $\mathsf{E}$ is the interpolation operator that interpolates grid data to the surface points, $\mathsf{s}_{\infty}$ is the streamfunction of the uniform flow. For a body translating at velocity $(U,V)$ and rotating at angular velocity $\Omega$, this streamfunction would be

\[\mathfrak{s}_{b,k} = U \mathfrak{r}_y - V \mathfrak{r}_x - \frac{1}{2} \Omega (\mathfrak{r}_x^2+\mathfrak{r}_y^2).\]

For later shorthand, we will denote the difference between the body motion streamfunction and interpolated uniform flow streamfunction by $\mathfrak{s}'_b \equiv \mathfrak{s}_b - \mathsf{E} \mathsf{s}_{\infty}$. The no-penetration constraint is enforced in the basic potential flow problem with the help of a vector of Lagrange multipliers, $\mathfrak{f}$, on the surface points. The modified potential flow problem is thus

\[\mathsf{Ls} + \mathsf{R}\mathfrak{f} = -\mathsf{w},\]

where $\mathsf{R}$ is the regularization operator that transfers data from surface data to nodes. $\mathsf{E}$ can be constructed (and we will assume it has) so that it is the transpose of the interpolation operator, $\mathsf{E} = \mathsf{R}^{T}$.

From the previous equation, it is clear (by simple comparison with the vorticity on the right-hand side) that the vector $\mathfrak{f}$ represents the strength of the discrete bound vortex sheet on the surface. Suppose we consider the bound vortex sheet $\gamma(s)$ that emerges from the analogous continuous problem on the undiscretized surface, where $s$ is the arc-length parameter along the surface. At each point $p$, the discrete solution $\mathfrak{f}$ is approximately equal to this continuous solution, multiplied by the length $\delta S_p$ of the small segment surrounding the point:

\[\mathfrak{e}_{p}^{T} \mathfrak{f} \approx \gamma(s) \delta S_p.\]

Thus, the potential flow problem in the presence of the impenetrable surface is

\[\begin{bmatrix} \mathsf{L} & \mathsf{R} \\ \mathsf{E} & 0 \end{bmatrix} \begin{pmatrix} \mathsf{s} \\ \mathfrak{f} \end{pmatrix} = \begin{pmatrix} -\mathsf{w} \\ \mathfrak{s}'_b - \mathfrak{s}_0 \end{pmatrix}.\]

This problem has the structure of a generic saddle-point problem. GridPotentialFlow.jl automatically formulates the no-penetration constraint on the streamfunction and solves the saddle system when it needs to. See the next page for the role of $\mathfrak{s}_0$ and how it will be treated internally.

We will now consider two examples. In these examples, we use the following grid.

Δx = 0.02
Lx = 4.0
xlim = (-Lx/2,Lx/2)
ylim = (-Lx/2,Lx/2)
g = PhysicalGrid(xlim,ylim,Δx);

Vortex near a cylinder

As a basic example, consider now a point vortex near a circular cylinder.

A body in GridPotentialFlow.jl is represented by PotentialFlowBody. This is a wrapper around a Body from RigidBodyTools.jl and can be endowed with a linear velocity, rotational velocity, circulation, and other properties using the appropriate keywords during construction or by mutating its fields using its setter methods. For our circular cylinder, we use the Circle shape.

Rc = Lx/4
Δs = 2Δx
body = PotentialFlowBody(Circle(Rc,Δs))

Potential flow body with 0 regularized edges
Shape and position: Circular body with 156 points and radius 1.0
   Current position: (0.0,0.0)
   Current angle (rad): 0.0
Linear velocity: (0.0, 0.0)
Angular velocity: 0.0
Current circulation: 0.0

We choose an initial position $R_v$ of the vortex at $3/2 R_c$ from the origin.

Rv = 3/2*Rc
Γv = 1.0
v = Vortex(Rv,0.0,Γv);

And we set an equal and opposite circulation about the body

setΓ(body,-Γv)

-1.0

We construct the VortexModel by using the optional bodies keyword. In this case, we are interested in the bound vortex sheet strength $\mathfrak{f}$, instead of just the streamfunction field. For this reason, we use solve, which in this case returns a ConstrainedIBPoissonSolution.

model = VortexModel(g,vortices=[v],bodies=[body]);
sol = solve(model);

We can then easily retrieve the streamfunction and the bound vortex sheet strength from the fields of the solution variable.

plot(sol.ψ,g)
plot!(body,fillcolor=:black,fillrange=0,fillalpha=0.25,linecolor=:black,linewidth=2)
scatter!([v.x],[v.y],color=:black,markersize=2,xlabel="x",ylabel="y")

For this example, we can easily calculate the analytical bound vortex sheet strength $\gamma$ and compare it with our numerical solution $f$. See the notebook in the examples folder for the analytical solution.

plot(sol.f./Δs,label="f/ds",xlabel="body point index")
plot!(γ,label="gamma")

Finally, we can also check if the vortex will move around the cylinder on a circular path if we advance the system in time, as predicted by analytical potential flow theory. The period of the circular motion can be easily determined by calculating the velocity induced by the image vortex at $R_c^2/R_v$:

Vv = Γv/(2π*(Rv-Rc^2/Rv))
T = 2π*Rv/Vv

49.34802200544679

For the time stepping, we again use the fourth-order Runge Kutta scheme from OrdinaryDiffEq.jl.

import OrdinaryDiffEq
function rhs(X,model,t)
    setvortexpositions!(model,X)
    Ẋ = vortexvelocities!(model)
    return Ẋ
end
X = getvortexpositions(model)
prob = OrdinaryDiffEq.ODEProblem(rhs,X,(0.0,T),model);
sol = OrdinaryDiffEq.solve(prob,dt=0.1,OrdinaryDiffEq.RK4(),dense=false,adaptive=false);
plot(body,fillcolor=:black,fillrange=0,fillalpha=0.25,linecolor=:black,linewidth=2)
plot!(map(s->s.u[1],sol.u),map(s->s.v[1],sol.u))
scatter!([v.x],[v.y],color=:black,markersize=2,xlabel="x",ylabel="y")

Moving cylinder

In this example we will create a model with a moving body.

Here we use the keyword U to add create a body with a linear velocity.

body = PotentialFlowBody(Circle(Rc,Δs), U=(1.0,0.0))
model = VortexModel(g,bodies=[body]);
sol = solve(model);

plot(sol.ψ,g)
plot!(body,fillcolor=:black,fillrange=0,fillalpha=0.25,linecolor=:black,linewidth=2,xlabel="x",ylabel="y")

Again, we can compare the discrete bound vortex sheet strength with the analytical vortex sheet strength.

plot(sol.f./Δs,label="f/ds",xlabel="body point index")
plot!(γ,label="gamma")

This page was generated using Literate.jl.