Difference added mass force (moving domain/moving body)

[ddeboerfluid]

Dear all,

I have noticed a significant difference between the added mass force from WaterLily when using a moving domain, compared to a moving body (see figure). Interestingly, after the acceleration stops the forces are very similar again. I noticed it for a plate in 3D first, but it seems to happen in 2D and with other bodies as well. From theoretical added mass values I found for a plate, it seems like the moving body case is approximately correct, and the moving domain case overestimates the force.

What I did:

1. 2D circle, constant acceleration, then constant velocity. All by mapping the body.
2. 2D circle with same parameters. All by using a moving domain.

Below I added the code generating the data for the figure (csv file with t* and c_d).

Thank you for helping and I hope you can find out what is causing this.

Moving body code:

# Importing
    using WriteVTK, BiotSavartBCs, WaterLily, StaticArrays, CSV, DataFrames
    include("CSV_writer.jl")


    #Simulation settings 

        #Limit time step by overwritting default function
        WaterLily.CFL(a::Flow) = WaterLily.CFL(a; Δt_max=0.2)

    # CSV writer settings 

        #Set attributes 
        t_attrib(a::Simulation, forces) = sum(a.flow.Δt[1:end-1])*a.U/a.L;
        d_attrib(a::Simulation, forces) = forces[1];

        #Make attributes dict 
        attributes_dict         =   Dict("Time [-]" => t_attrib,
                                        "Drag [-]" => d_attrib) 

    #Circle SDF 
        function TwoD_circle_sdf(x, radius, center)
            return √sum(abs2, x .- center) - radius
        end
    
    #Motion function 
        function U_motion(x,t, U, L, acceleration_ND)

            #Get t_i 
            t_i                 =   t*U/L 

            #Get t_acc
            t_acc               =   U/acceleration_ND

            #If accelerating, define motion for accelerating part 
            if t_i<t_acc
                
                x - L*SA[0.5*acceleration_ND*t_i^2,0]
                
            else
                x - L*SA[0.5*acceleration_ND*t_acc^2 + U*(t_i-t_acc),0]
                
            end
        end

    #Function to simulate circle with constant velocity (2D)
    function Circle_acc_to_const(L,Re,domain,acceleration_ND;U=1,ϵ=0.5, f=Array)
        
        #Set radius and center 
        diameter, center  =   L, SA[domain[1]*L/2, domain[2]*L/2]

        #Set sdf 
        sdf     =   (x,t) -> TwoD_circle_sdf(x,diameter/2, center)

        map     =   (x,t) -> U_motion(x,t, U, L, acceleration_ND)

        #Write domain tuple
        domain_tuple    =   tuple(domain*L...)

        #Define simulation and return it
        Simulation(domain_tuple,(0,0),L;U,ν=U*L/Re,body=AutoBody(sdf, map),ϵ, mem=f)
    end

    #Run 
    function run(sim, case_name, t_max_star, wr_CSV, ω)

        #Set t0 to be zero
        t₀ = 0.0 

        #Do actual simulation and store fields and forces
        @time for tᵢ in range(t₀,t₀+t_max_star;step=0.05)

            #Get current time (sum of delta t until now)
            t = sum(sim.flow.Δt[1:end-1])

            #While time is lower than normalised ti, keep running
            while t < tᵢ*sim.L/sim.U

                #Measure 
                measure!(sim,t)
            
                #Calculate forces 
                forces = 2*WaterLily.pressure_force(sim.flow.p,sim.flow.f,sim.body,t)/(sim.L)
                
                #Append parameters to arrays in dict  
                CSV_write(wr_CSV, sim, forces)

                #Evolve flow
                biot_mom_step!(sim.flow,sim.pois,ω)

                #Evolve time
                t += sim.flow.Δt[end]
            end

            #Print progress
            println("tU/L=",round(tᵢ,digits=2),", Δt=",round(sim.flow.Δt[end],digits=3), ", -",round(t*sim.U/sim.L, digits=3))
        end

        #Close writer
        CSV_close(wr_CSV)
    end


#Set parameters 

    #Set Re 
    Re      =   80 #-

    #Set L (diameter)
    L       =   64 #-

    #Set domain 
    domain  =   [10,5]

    #Set nondimensional acceleration 
    acceleration_ND     =   5

#Make writers and set folders 

    #Set case name 
    case_name   =   "Test_am_map64-10x5"

    #Make CSV writer 
    CSV_writer      =   get_CSVWriter("forces_moving_body.csv", attributes_dict)

#Construct sim 
    #Construct simulation 
    sim     =   Circle_acc_to_const(L, Re, domain, acceleration_ND)

    #Get ω for Biot-Savart
    ω       =   MLArray(sim.flow.σ)

#Run simulation 

    t_star_max  =   0.5

    #Run
    run(sim, case_name, t_star_max, CSV_writer, ω)

Moving domain code:

# Importing
    using WriteVTK, BiotSavartBCs, WaterLily, StaticArrays, CSV, DataFrames
    include("CSV_writer.jl")

#Simulation settings 

    #Limit time step by overwritting default function
    WaterLily.CFL(a::Flow) = WaterLily.CFL(a; Δt_max=0.2)

# CSV writer settings 

    #Set attributes 
    t_attrib(a::Simulation, forces) = sum(a.flow.Δt[1:end-1])*a.U/a.L;
    d_attrib(a::Simulation, forces) = forces[1];

    #Make attributes dict 
    attributes_dict         =   Dict("Time [-]" => t_attrib,
                                    "Drag [-]" => d_attrib) 

#Circle SDF 
    function TwoD_circle_sdf(x, radius, center)
        return √sum(abs2, x .- center) - radius
    end
    
#Motion function 
    function U_linear_to_constant(t, U, L, acceleration_ND)

        #Get t_i 
        t_i                 =   t*U/L 

        #Get t_acc
        t_acc               =   U/acceleration_ND

        #If accelerating, define motion for accelerating part 
        if t_i<t_acc
            updated_vel     =   acceleration_ND*t_i  
        
        else
            updated_vel     =   U

        end
    end

#Function to simulate circle with acceleration (2D)
    function Circle_acc_to_const(L,Re,domain,acceleration_ND;U=1,ϵ=0.5, f=Array)
        
        #Set radius and center 
        diameter, center  =   L, SA[domain[1]*L/2, domain[2]*L/2]

        #Set sdf 
        sdf     =   (x,t) -> TwoD_circle_sdf(x,diameter/2, center)

        #Get U of domain (horizontal velocity)
        U_domain    =   t -> U_linear_to_constant(t, U, L, acceleration_ND)

        #Get velocity vector of domain
        Ut(i,t::T) where T = i==1 ? convert(T,-U_domain(t)) : zero(T)

        #Write domain tuple
        domain_tuple    =   tuple(domain*L...)

        #Define simulation and return it
        Simulation(domain_tuple,Ut,L;U,ν=U*L/Re,body=AutoBody(sdf),ϵ, mem=f)
    end

    #Run 
    function run(sim, case_name, t_max_star, wr_CSV, ω)

        #Set t0 to be zero
        t₀ = 0.0 

        #Do actual simulation and store fields and forces
        @time for tᵢ in range(t₀,t₀+t_max_star;step=0.05)

            #Get current time (sum of delta t until now)
            t = sum(sim.flow.Δt[1:end-1])

            #While time is lower than normalised ti, keep running
            while t < tᵢ*sim.L/sim.U

                #Measure 
                measure!(sim,t)
            
                #Calculate forces 
                forces = 2*WaterLily.pressure_force(sim.flow.p,sim.flow.f,sim.body,t)/(sim.L)
                
                #Append parameters to arrays in dict  
                CSV_write(wr_CSV, sim, forces)

                #Evolve flow
                biot_mom_step!(sim.flow,sim.pois,ω)

                #Evolve time
                t += sim.flow.Δt[end]
            end

            #Print progress
            println("tU/L=",round(tᵢ,digits=2),", Δt=",round(sim.flow.Δt[end],digits=3), ", -",round(t*sim.U/sim.L, digits=3))
        end

        #Close writer
        CSV_close(wr_CSV)
    end


#Set parameters 

    #Set Re 
    Re      =   80 #-

    #Set L (diameter)
    L       =   64 #-

    #Set domain 
    domain  =   [10,5]

    #Set nondimensional acceleration 
    acceleration_ND     =   5

#Make writers and set folders 

    #Set case name 
    case_name   =   "Test_am64-10x5"

    #Make CSV writer 
    CSV_writer      =   get_CSVWriter("forces_moving_domain.csv", attributes_dict)

#Construct sim 
    #Construct simulation 
    sim     =   Circle_acc_to_const(L, Re, domain, acceleration_ND)

    #Get ω for Biot-Savart
    ω       =   MLArray(sim.flow.σ)

#Run simulation 

    t_star_max  =   0.5
    
    #Run
    run(sim, case_name, t_star_max, CSV_writer, ω)

CSV writer code:

#Code that defines CSV writer structure 

using CSV
using WaterLily 

#CSV writer structure 
struct CSVWriter
    fname         :: String
    directory     :: String
    attributes    :: Dict{String,Function}
    output        :: Dict{String,Array}
end 

#Define function to construct correct CSV Writer 
function get_CSVWriter(directory, attributes,fname="WaterLily")

    #Make empty output dict 
    output  =   Dict()

    #For all keys in attributes, set same keys in output dict and fill with empty arrays
    for key in collect(keys(attributes)) 
        output[key]   =   []    
    end

    CSVWriter(fname,directory,attributes,output)
end

#Define write csv function
function CSV_write(w::CSVWriter, a::Simulation, forces)
    #For all attributes, calculate and store in output 
    for (key,func) in w.attributes
        #Compute and store to output 
        push!(w.output[key], func(a,forces))
    end
end

#Define function to close csv (and save it)

function CSV_close(w::CSVWriter)
    CSV.write(w.directory, DataFrame(w.output), delim = ';')
end

[marinlauber]

Hey Dirk,

Thanks for picking this up. Maybe next time try to come up with a true MWE, without fileIO and this kind of stuff. It took me a while to make sure I was doing the same as you. Here is my MWE script

using WaterLily,StaticArrays
function make_sim_acc(; N=128, R=32, a=0.5, U=1, Re=1e3, mem=Array)
    circle(x,t) = √sum(abs2,x.-2N)-R # move the center
    Ut(i,t::T) where T = i==1 ? convert(T,ifelse(t≤U*R/a,a*t/R,U)) : zero(T) # velocity BC
    Simulation((4N,4N), Ut, R; U,ν=U*R/Re, body=AutoBody(circle), mem)
end
function make_sim(; N=128, R=32, a=0.5, U=1, Re=1e3, mem=Array)
    circle(x,t) = √sum(abs2,x.-2N)-R # move the center
    s(t) = ifelse(t≤U/a,0.5a*t^2,U*(t-0.5U/a)) # displacement
    move(x,t) = x + SA[R*s(t/R)-2R,0]
    Simulation((4N,4N), (0,0), R; U,ν=U*R/Re, body=AutoBody(circle,move), mem)
end
include("src/TwoD_plots.jl")

N = 2^7; R = N/4
acc = false
sim = acc ? make_sim_acc(mem=Array;N,R,Re=250) : make_sim(mem=Array;N,R,Re=250)
domain = inside(sim.flow.p)

forces = []
@gif for t in range(0,5;step=0.05)
    while sim_time(sim)<t
        measure!(sim,sum(@views(sim.flow.Δt))) # this is the same time as the force computation
        mom_step!(sim.flow,sim.pois)
        f = 2WaterLily.pressure_force(sim)
        push!(forces,[sim_time(sim),f[1]])
    end
    flood(sim.flow.p[domain],clims=(-1,1)); body_plot!(sim)
    @show t
end

if you use those you will realise the pressure field are different in the acceleration phase. This will directly translate into different forces and thus added-mass. I am not quite sure what the issue is, I suspect the pressure solver tolerances.

The forces match after the steady phase, but obviously not in the accelerating frame (y1=accelerating frame, y2 accelerating body)

[b-fg]

Thanks @marinlauber for the detailed reply. @ddeboerfluid, does this solve your issue?

[ddeboerfluid]

Hi @b-fg, I talked with @marinlauber about it and as I am interested mostly in the post-acceleration phase, it will have no influence on my results (the flow fields after the acceleration are identical). However, I think the issue should still be open as it is still present: the drag force during acceleration is overestimated (seems to be an integer factor) when accelerating the flow.