# 6. Force and the added mass

## Force

This part introduces the calculation of the force and moment on the body through the negative rate of change of impulse in the fluid. The continuous expressions for linear and angular impulse (about the origin) are, in two dimensions,

$\begin{align}
    \boldsymbol{P} &= \int_{V_f} \boldsymbol{x}\times \boldsymbol{\omega}\,\mathrm{d}V + \int_{S_b} \boldsymbol{x}\times \left( \boldsymbol{n}\times\boldsymbol{v}\right)\,\mathrm{d}S \\
    \boldsymbol{\Pi}_0 &= \frac{1}{2}\int_{V_f} \boldsymbol{x}\times \left(\boldsymbol{x}\times\boldsymbol{\omega}\right)\,\mathrm{d}V + \frac{1}{2} \int_{S_b} \boldsymbol{x}\times \left[ \boldsymbol{x}\times\left( \boldsymbol{n} \times \boldsymbol{v} \right) \right] \,\mathrm{d}S
\end{align}$

If there is only a single body, then the force and moment (about the origin) exerted by the fluid on that body are given by

$\boldsymbol{F} = -\rho \frac{\mathrm{d}\boldsymbol{P}}{\mathrm{d}t}, \qquad \boldsymbol{M}_{0} = -\rho \frac{\mathrm{d}\boldsymbol{\Pi}_{0}}{\mathrm{d}t},$

where $\rho$ is the fluid density. In the two-dimensional applications of this package, angular impulse and the moment have only a single component, e.g., $\boldsymbol{\Pi}_0 =\Pi_0\boldsymbol{e}_z$, where $\boldsymbol{e}_z$ is the unit vector out of the plane.

#### Point vortices past a cylinder

To verify our method of calculating the impulse, we create a model of two point vortices of equal and opposite circulation positioned at both sides of the $x$-axis, with the $x$-axis as axis of symmetry, and left of a circular cylinder positioned at the origin. If the top vortex has a positive circulation, the vortices propel each other towards the cylinder and convect past it. This configuration has an analytical solution for the trajecory and impulse.

We use a rectangular grid.


```julia
using GridPotentialFlow
using Plots
Lx = 10.0
Ly = 4.0
xlim = (-Lx/2,Lx/2)
ylim = (-Ly/2,Ly/2)
g = PhysicalGrid(xlim, ylim, 0.04)
```




    PhysicalGrid{2}((252, 104), (126, 52), 0.04, ((-5.0, 5.0), (-2.04, 2.04)))




```julia
Rc = 1
circle = Circle(Rc,2*cellsize(g))
pfb = PotentialFlowBody(circle)

Δs = dlength(circle);
```

The initial spacing between the vortices is `y∞`. After the vortices pass the cylinder, they should return to this spacing.


```julia
y∞ = Rc/2
v1 = Vortex(-Lx/2+10*cellsize(g),y∞,1.0);
v2 = Vortex(-Lx/2+10*cellsize(g),-y∞,-1.0);
```

Analytical trajectory


```julia
x(r) = sqrt(r^2-(y∞^2*(r^2-1)^2)/((r^2-1)^2-4*y∞^2))
y(r) = y∞*(r^2-1)/sqrt((r^2-1)^2-4*y∞^2)
s = 0:0.001:π/2
rmin = 1.4375649
c1 = rmin/(-Lx/2+rmin)
c2 = -Lx/2*c1
r = c2 ./ (cos.(s).-c1)
x_trajectory = -reverse(x.(r)); append!(x_trajectory,x.(r))
y_trajectory_upper = reverse(y.(r)); append!(y_trajectory_upper,y.(r))
y_trajectory_lower = -y_trajectory_upper;
```

Analytical impulse


```julia
Px_func(x,y) = y - y/(x^2+y^2) + y - y/(x^2+y^2);
Py_func(x,y) = -x + x/(x^2+y^2) + x - x/(x^2+y^2);
Px_exact = Px_func.(x_trajectory,y_trajectory_upper);
Py_exact = Py_func.(x_trajectory,y_trajectory_upper);
```

We create the vortex model with these two point vortices and the circle and advance the position of the point vortices over some time. We use `OrdinaryDiffEq.RK4()` for the time stepping and we compute the impulse history with `impulse` in a post-processing step.


```julia
model = VortexModel(g,bodies=[pfb],vortices=[v1,v2]);

function rhs(X,model,t)
    setvortexpositions!(model,X)
    Ẋ = vortexvelocities!(model)
    return Ẋ
end

import OrdinaryDiffEq
X = getvortexpositions(model)
prob = OrdinaryDiffEq.ODEProblem(rhs,X,(0.0,64.0),model);
sol = OrdinaryDiffEq.solve(prob,dt=0.1,OrdinaryDiffEq.RK4(),dense=false,adaptive=false);
```

Using the history of the vortex positions in `sol.u`, we can calculate the impulse history.


```julia
Px_numerical_hist = []
Py_numerical_hist = []
for u in sol.u
    setvortexpositions!(model,u)
    Px, Py = impulse(model)
    push!(Px_numerical_hist,Px)
    push!(Py_numerical_hist,Py)
end
```

When we compare the trajectories and impulse history, the numerical and analytical solution should match closely, which is indeed the case.


```julia
plot(circle,fillcolor=:black,fillrange=0,fillalpha=0.25,linecolor=:black,linewidth=2,xlabel="x",ylabel="y")
scatter!(model.vortices.x,model.vortices.y,color=:red)
plot!(x_trajectory,y_trajectory_upper,linecolor=:red,label="exact")
plot!(x_trajectory,y_trajectory_lower,linecolor=:red,label="")
plot!(map(s->s.u[1],sol.u),map(s->s.v[1],sol.u),color=:blue,linestyle=:dash,label="simulated")
plot!(map(s->s.u[2],sol.u),map(s->s.v[2],sol.u),color=:blue,linestyle=:dash,label="")
```


```julia
plot(x_trajectory,Px_exact,color=:red,label="exact",xlabel="x",ylabel="Px")
plot!(map(s->s.u[1],sol.u),Px_numerical_hist,color=:blue,linestyle=:dash,label="simulated")
```

#### Flat plate

We can apply this method also to our flat plate example and compare it to the Biot-Savart method from the `PotentialFlow.jl` package.

Because we create the `PotentialFlow.jl` model with a moving plate instead of a uniform flow, we specify a translational velocity. The `GridPotentialFlow.jl` model will still use a moving body coordinate system and thus a uniform flow that equals the negative translational velocity.


```julia
c = 1.0;   # chord length

ċ = -1.0;  # translational velocity
α = -π/3;   # angle of attack

σLE = 0.0;
σTE = 0.0;

Δt = 5e-2;
tf = 1.0;
T = 0.0:Δt:tf;
```

Potential Flow


```julia
using PotentialFlow
function compute_ẋ!(ẋ, x, t)
    plate, ambient_sys = x
    motion = ẋ[1]
    motion.ċ, motion.c̈, motion.α̇ = motion.kin(t)
    Plates.enforce_no_flow_through!(plate, motion, ambient_sys, t)
    reset_velocity!(ẋ, x)
    self_induce_velocity!(ẋ, x, t)
end
function shed_new_vorticity!(blobs, plate, motion, t, lesp = 0.0, tesp = 0.0)
    z₊ = (blobs[end-1].z + 2plate.zs[end])/3
    z₋ = (blobs[end].z + 2plate.zs[1])/3
    blob₊ = PotentialFlow.Vortex.Blob(z₊, 1.0, δ)
    blob₋ = PotentialFlow.Vortex.Blob(z₋, 1.0, δ)
    Plates.enforce_no_flow_through!(plate, motion, blobs, t)
    Γ₊, Γ₋, _, _ = Plates.vorticity_flux!(plate, blob₊, blob₋, t, lesp, tesp);
    push!(blobs, PotentialFlow.Vortex.Blob(z₊, Γ₊, blobs[1].δ), PotentialFlow.Vortex.Blob(z₋, Γ₋, blobs[1].δ))
end
N = 128 # number of plate control points (distributed along a extrema Chebyshev grid)
δ = 0.01
plate = PotentialFlow.Plate(N, c, zero(ComplexF64), α)
motion = Plates.RigidBodyMotion(ċ, 0.0);
```

    WARNING: using PotentialFlow.Vortex in module ##1127 conflicts with an existing identifier.


Initial step


```julia
Δz₀ = im*3Δt*exp(im*plate.α) # vectors perpendicular to the plate
z₋, z₊ = plate.zs[[1,N]] # LE and TE
blobs = PotentialFlow.Vortex.Blob.(Δz₀ .+ [z₊, z₋], 1.0, δ) # First two point vortices are placed close to the LE and TE with unit strength
Plates.enforce_no_flow_through!(plate, motion, (), 0)
Γ₊, Γ₋, _, _ = Plates.vorticity_flux!(plate, blobs[1], blobs[2], 0.0, σLE, σTE); # Determine strength of first two vortices
blobs = PotentialFlow.Vortex.Blob.(Δz₀ .+ [z₊, z₋], [Γ₊, Γ₋], δ) # Create first two point vortices now with calculated strengths
sys₀ = (plate, blobs)
sys = deepcopy(sys₀)
sys₊ = deepcopy(sys₀) # Used for storage during time-marching
ẋs = (motion, allocate_velocity(blobs))
imp = ComplexF64[];
t_hist = Float64[]
global t = 0
push!(t_hist,t)
```




    1-element Vector{Float64}:
     0.0



Time stepping


```julia
for tloc in T[2:end]
    global t += Δt
    global sys
    global sys₊
    push!(imp,Elements.impulse(sys))
    push!(t_hist,t)
    local plate, ambient_ω = sys
    local motion, ambient_u = ẋs
    resize!(sys₊[2], length(sys[2]))
    forward_euler!(sys₊, sys, t, Δt, compute_ẋ!, advect!, ẋs)
    sys, sys₊ = sys₊, sys
    shed_new_vorticity!(sys[2], sys[1], ẋs[1], t, σLE, σTE)
end
force = -diff(imp)/Δt;

Δx = 0.01
xlim = (-0.5,2)
ylim = (-1,1)
g = PhysicalGrid(xlim,ylim,Δx);
```


```julia
dsdx = 2
plate = RigidBodyTools.Plate(c,dsdx*Δx)
pfb = PotentialFlowBody(plate,edges=[1,length(plate)],σ=[SuctionParameter(σLE),SuctionParameter(σTE)])
transform = RigidTransform((0.0,0.0),-π/3)
transform(plate);
Δs = dlength(plate)
maximum(Δs/Δx);
```

As before, we create the initial vortices based on the time step and uniform flow and define a function that create the new point vortices for the subsequent time steps.


```julia
firstvLE = GridPotentialFlow.Vortex(plate.x[1]+3Δt*(-ċ)*cos(plate.α+π/2),plate.y[1]+3Δt*(-ċ)*sin(plate.α+π/2),0.0);
firstvTE = GridPotentialFlow.Vortex(plate.x[end]+3Δt*(-ċ)*cos(plate.α+π/2),plate.y[end]+3Δt*(-ċ)*sin(plate.α+π/2),0.0);
```


```julia
function createsheddedvortices(plate,oldvortices)

    vLE = GridPotentialFlow.Vortex(2/3*plate.x[1]+1/3*oldvortices[end-1].x,2/3*plate.y[1]+1/3*oldvortices[end-1].y,0.0)
    vTE = GridPotentialFlow.Vortex(2/3*plate.x[end]+1/3*oldvortices[end].x,2/3*plate.y[end]+1/3*oldvortices[end].y,0.0)

    return vLE,vTE
end
```




    createsheddedvortices (generic function with 1 method)



The model free stream should be the negative of the translational velocity of the plate.


```julia
model = VortexModel(g,bodies=[pfb],vortices=[firstvLE,firstvTE],U∞=(-ċ,0.0));
```

Then we enter a time stepping loop and record the impulse every time we create new vortices. Note that because we release new vortices at every time step, it is not straightforward to use a multi-step time integration method like RK4 to advance the positions of the point vortices. Instead, we use forward Euler and the size of the positions vector increases every time step by the number of new vortices we insert. To advance the system one time step, we first need to solve the saddle-point system, set the strengths of the new vortices, subtract that circulation from the bodies, and calculate the vortex velocities with the saddle-point system solution. For convenience, these four operations are carried out internally if we would call `vortexvelocities!(model)`. However, if we want to have access to the solution of the saddle-point system, we can explicitly perform each operation as is done in the following time stepping loop. Note that we have to set the strengths of the new vortices before calling the `impulse` function, otherwise the new vortices won't be included in the calculation.


```julia
Px_hist = Float64[];
Py_hist = Float64[];
sol = solve(model);
for tloc in T[2:end]
    X = getvortexpositions(model) # gets bigger every time step because we add vortices
    Ẋ = deepcopy(X)

    solve!(sol, model)
    setvortexstrengths!(model, sol.δΓ_vec, length(X.u)-1:length(X.u))
    subtractcirculation!(model.bodies, sol.δΓ_vec)
    Px, Py = impulse(model)

    vortexvelocities!(Ẋ, model, sol.ψ)
    X .= X .+ Ẋ*Δt
    setvortexpositions!(model, X)

    vLE, vTE = createsheddedvortices(plate,model.vortices)
    pushvortices!(model,vLE,vTE)
    push!(Px_hist,Px)
    push!(Py_hist,Py)
end
```

    ┌ Warning: Assignment to `X` in soft scope is ambiguous because a global variable by the same name exists: `X` will be treated as a new local. Disambiguate by using `local X` to suppress this warning or `global X` to assign to the existing global variable.
    └ @ string:5


The force history can be obtained by taking finite differences of the impulse history.


```julia
Fx_hist = -diff(Px_hist)/Δt;
Fy_hist = -diff(Py_hist)/Δt;
```

We can now compare the positions of the point vortices by shifting the `PotentialFlow.jl` solution by `-tf*ċ` such that the origin of the frame of reference coincides with the center plate. Superimposing the `GridPotentialFlow.jl` solution then shows that the positions of the point vortices agree very well.


```julia
plot(plate,fillcolor=:black,fillrange=0,fillalpha=0.25,linecolor=:black,linewidth=2,xlim=xlim,ylim=ylim,xlabel="x",ylabel="y")
scatter!(real.((v->v.z).(sys[2])).-tf*ċ,imag.((v->v.z).(sys[2])),color=:red,markersize=4,label="PotentialFlow.jl")
scatter!(model.vortices.x,model.vortices.y,color=:blue,markersize=2,label="GridPotentialFlow.jl")
```

The vertical impulse and the vertical force (lift) can also be compared and show good agreement as well.


```julia
plot(xlabel="t",ylabel="Py")
plot!(T[1:end-1],imag.(imp),color=:blue,label="PotentialFlow.jl")
plot!(T[1:end-1],Py_hist,color=:red,label="GridPotentialFlow.jl")
```


```julia
plot(xlabel="t",ylabel="Fy")
plot!(T[2:end-1],imag.(force),color=:blue,label="PotentialFlow.jl")
plot!(T[2:end-1],Fy_hist,color=:red,label="GridPotentialFlow.jl")
```

## Added mass

The added mass tensor provides a measure of the inertial influence of the fluid on the body in response to changes in the body's translational or rotational motion. The coefficients of the added mass tensor of a body are obtained by computing the impulse components associated with a unit-valued component of motion. The motion's influence is both direct, via the surface velocity, and indirect, in the bound vortex sheet that develops on the surface.

The added mass for simple geometries can easiliy be calculated in an analytical way. For example, the entries of the translational added mass tensor for an ellipse with semi-major axis $a$ and semi-minor axis $b$ are $m_{xx} = \rho \pi b^2$ for motion in the $x$ direction, $m_{yy} = \rho \pi a^2$ for motion in the $y$-direction, with the diagonal entries $m_{xy}=m_{yx}=0$. By using `addedmass` from this package we can approximate these results numerically. In this case, we have one body and the method will return a matrix with the entries

$\begin{bmatrix}
m_{xx} & m_{xy} \\
m_{yx} & m_{yy}
\end{bmatrix}.$


```julia
Lx = 2.0
xlim = (-Lx/2,Lx/2)
ylim = (-Lx/2,Lx/2)
g = PhysicalGrid(xlim,ylim,0.01)
Δx = cellsize(g);


a = 0.5
b = 0.25
ellipse = Ellipse(a,b,Δx)
pfb = PotentialFlowBody(ellipse)
model = VortexModel(g,bodies=[pfb])
M = GridPotentialFlow.addedmass(model)
```




    2×2 Matrix{Float64}:
      0.209533     -2.32453e-16
     -2.01445e-16   0.804095



We can compare the values using the `@test` macro.


```julia
using Test
@test isapprox(M[1,1], π*b^2, rtol=1e-1)
```




    [32m[1mTest Passed[22m[39m




```julia
@test isapprox(M[2,2], π*a^2, atol=1e-1)
```




    [32m[1mTest Passed[22m[39m



In the case of $N$ bodies, `addedmass` returns the translational added mass tensor of size $2N$-by-$2N$. As an example, we will compute the added mass tensor for an array of cylinders. To compare it with an analytically obtained solution, we will calculate the added mass matrix tensor for the following gap-to-radius ratios.


```julia
GRratios = [0.1,0.2,0.3,0.4,0.5,0.6,0.8,1.0,1.2,1.4,1.6];


function rectangulararray(rows,columns,spacing)
     centers = VectorData(rows*columns)
     for j in 1:columns
        for i in 1:rows
            xc = (i-1)*spacing
            yc = (j-1)*spacing
            centers.u[(j-1)*rows+i] = xc
            centers.v[(j-1)*rows+i] = yc
        end
    end
    centers.u .-= mean(centers.u)
    centers.v .-= mean(centers.v)
    return centers
 end
```




    rectangulararray (generic function with 1 method)



For this case, the array consists of three rows of three cylinders.


```julia
n = 100
rows = 3
columns = 3
N = rows*columns
R = 0.20
𝒱 = π*R^2
bodies = fill(Circle(R,2*cellsize(g)),N);
Δs = minimum(dlength(bodies[1]));
```

We loop over the gap-to-radius ratios, position the bodies, and compute the ratio of the largest eigenvalue to the largest diagonal element of the translational added mass tensor. To position the bodies, we defined the method `rectangulararray` (see notebook in the examples folder).


```julia
using Statistics: mean
using LinearAlgebra: eigen
```


```julia
λoverMratios = zeros(length(GRratios))
for idx in 1:length(GRratios)
    global bodies
    gap = GRratios[idx]*R
    spacing = 2*R+gap
    bodycenters = rectangulararray(rows,columns,spacing)
    for i in 1:N
        Tf = RigidTransform((bodycenters.u[i],bodycenters.v[i]),0.0)
        global bodies[i] = Tf(deepcopy(bodies[i]))
    end
    body_list = PotentialFlowBody.(bodies)
    vm = VortexModel(g,bodies=body_list);
    M = GridPotentialFlow.addedmass(vm)/𝒱
    eigM = eigen(M);
    max_eig_value_coef = maximum(real(eigM.values))
    max_self_added_mass_coef = maximum(M)
    λoverMratios[idx] = max_eig_value_coef/max_self_added_mass_coef
end
```

    ┌ Warning: Assignment to `M` in soft scope is ambiguous because a global variable by the same name exists: `M` will be treated as a new local. Disambiguate by using `local M` to suppress this warning or `global M` to assign to the existing global variable.
    └ @ string:13


The same ratios, obtained in an analytical way, are available in literature [^1] and can be used to verify the numerical values.


```julia
nineRodArrayChen = [2.3637,2.2092,2.1007,2.0120,1.9350,1.8665,1.7494,1.6531,1.5732,1.5066,1.4508];
plot(xlabel="GR",ylabel="λ/max(Mij)")
plot!(GRratios,nineRodArrayChen,label="Chen1975")
plot!(GRratios,λoverMratios,label="GridPotentialFlow.jl")
```

[^1]: Chen, S. S. (1975) "Vibration of nuclear fuel bundles," *Nuclear Engineering and Design*, 35 (3), 399-–422.

---

*This notebook was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*

919d718850ed6d1747daff0f0cb791a630a520f4
ipynb
Jupyter Notebook
examples/6.-Force-and-the-added-mass.ipynb
JuliaIBPM/GridPotentialFlow.jl
b24127d036ed54a588a11ebedb35544acb8dea53
MIT
2021-02-25T02:10:16.000Z
2021-12-09T22:51:52.000Z
examples/6.-Force-and-the-added-mass.ipynb
JuliaIBPM/GridPotentialFlow.jl
b24127d036ed54a588a11ebedb35544acb8dea53
MIT
2021-04-01T20:39:56.000Z
2021-06-15T17:46:44.000Z
examples/6.-Force-and-the-added-mass.ipynb
JuliaIBPM/GridPotentialFlow.jl
b24127d036ed54a588a11ebedb35544acb8dea53
MIT
Qwen/Qwen-72B
1. YES
2. YES