Sorry, my ignoramus questions about the Waterlily again

[tiinaroose]

In the moving sphere demo TwoD_cylinder_VIV.jl
there are lines:

initial condition FSI

p0=radius/3; v0=0; a0=0; t0=0

# motion function uses global var to adjust
posx(t) = p0 + (t-t0)*v0

# motion definition
map(x,t) = x - SA[0,posx(t)]

Why does this imply any movement? v0=0? I've made the ellipse move up and down just like the sphere, but I want to make it now rotate. However, I do not understand how the up and down movement is implemented in the first place in WaterLily.

[marinlauber]

Hi @tiinaroose,

Sorry for the late reply, I missed this issue.

This example is a bit advanced, it's perfectly normal not to get it on first read.

In WaterLily, if you are using an AutoBody, geometries are defined through signed-distance-function sdf(x,t) and the position through a mapping map(x,t). Internally, to impose the kinematic boundary condition at the wet interface between the fluid and solid, we need the solid's position, which we have from the map, and its velocity, that we must compute. WaterLily uses Automatic Differentiation (AD) to compute the velocity at the interface by differentiating the map function. This is all packed in the measure!(AutoBody, x, t) function that you can find here.

In the case of prescribed motion, this works fine because we know the motion and the velocity, and if we can express them in terms of a simple function, we are all set. For FSI, this is different. The map function is unknown and cannot, in general, be expressed simply.

The trick I am doing here is to compute the position of the cylinder by solving the 1D equation of motion, and then, I use the position to map the cylinder to the correct location through a map function that uses the global variable p0. This is what the lines you are referencing do. Solving the 1DoF equations is simple and is done every time step within the code

# compute motion and acceleration 1DOF
Δt = sim.flow.Δt[end]
accel = (force[2]- k*p0 + mₐ*a0)/(m + mₐ)
p0 += Δt*(v0+Δt*accel/2.) 
v0 += Δt*accel
a0 = accel

Now the trick is that if we were to differentiate this map function w.r.t to time to get the velocity

map(x,t) = x - SA[0,p0]

we would get a constant zero velocity. Which is not correct, and the kinematic BC would be violated. The trick to achieving this is to ensure that AutoDiff (which essentially uses a finite difference on the function map w.r.t time to get the velocity) gets the correct velocity. This is why I substituted the second part

map(x,t) = x - SA[0, p0 + (t-t0) * v0]

by setting t0=t at this line I obtain a function that evaluated at t gives me the position, and when AutoDiff evaluated it at t±Δt, I get the velocity v0 that I computed from the FSI

(x - [0, p0+(t+Δt-t0)*v0]) - (x - [0, p0+(t-Δt-t0)*v0]) / 2Δt = [0, (t+Δt-t0)*v0 - (t-Δt-t0)*v0] / 2Δt = [0, v0]

This is a trick to use AutoBody with non-prescribed motion.

Now, for your case with a rotating ellipse, you have to add rotation on top of the linear motion in much the same fashion

function map(x,t)
    α = α₀ + rot
    x₀ = x .- center; 
    R = SA[cos(α) sin(α); -sin(α) cos(α)]
    R*(x₀ - pos) - (t-t_init).*vel .- (t-t_init)*ω
end

Which moves and rotates the ellipse. To get the rotation angle rot and angular velocity ω, you need to solve the moment equation around the ellipse in addition to the force balance.

I actually have an example of a free-falling ellipse that does all that. I can send it if you want.

I hope this helps.