diff --git a/UI_runner.jl b/UI_runner.jl new file mode 100644 index 00000000..5cc07e9e --- /dev/null +++ b/UI_runner.jl @@ -0,0 +1,3 @@ +include("src/UnsteadyFlowSolvers.jl") + +UnsteadyFlowSolvers.runUI() diff --git a/bin/airfoil.dat b/bin/airfoil.dat new file mode 100644 index 00000000..5a5be3cc --- /dev/null +++ b/bin/airfoil.dat @@ -0,0 +1,33 @@ +1.0000 0.0020 +0.9500 0.0165 +0.9000 0.0300 +0.8000 0.0524 +0.7000 0.0732 +0.6000 0.0898 +0.5000 0.1026 +0.4000 0.1108 +0.3000 0.1128 +0.2000 0.1082 +0.1500 0.1008 +0.1000 0.0890 +0.0750 0.0800 +0.0500 0.0694 +0.0250 0.0540 +0.0125 0.0430 +0.0000 0.0230 +0.0125 0.0094 +0.0250 0.0060 +0.0500 0.0022 +0.0750 0.0006 +0.1000 0.0000 +0.1500 0.0002 +0.2000 0.0014 +0.3000 0.0034 +0.4000 0.0048 +0.5000 0.0056 +0.6000 0.0060 +0.7000 0.0050 +0.8000 0.0034 +0.9000 0.0016 +0.9500 0.0008 +1.0000 0.0000 diff --git a/bin/config.txt b/bin/config.txt new file mode 100644 index 00000000..280212d5 --- /dev/null +++ b/bin/config.txt @@ -0,0 +1,15 @@ +Airfoil: bin/airfoil.dat +Reynolds#: 1e5 +Mach: .5 + +ChangePanels: true +Panels: 200 +ChangeIterations: true +Iterations: 200 + +AlphaStart: 0 +AlphaEnd: 60 +AlphaStep: 1 + +XfoilPath: bin\xfoil.exe +PolarFile: bin\polar.txt diff --git a/bin/xfoil.exe b/bin/xfoil.exe new file mode 100644 index 00000000..6b13ab8c Binary files /dev/null and b/bin/xfoil.exe differ diff --git a/builddir/7z.dll b/builddir/7z.dll new file mode 100644 index 00000000..be29515b Binary files /dev/null and b/builddir/7z.dll differ diff --git a/builddir/BugpointPasses.dll b/builddir/BugpointPasses.dll new file mode 100644 index 00000000..90c3fc25 Binary files /dev/null and b/builddir/BugpointPasses.dll differ diff --git a/builddir/LLVMHello.dll b/builddir/LLVMHello.dll new file mode 100644 index 00000000..9a1f28a1 Binary files /dev/null and b/builddir/LLVMHello.dll differ diff --git a/builddir/UNSflow.exe b/builddir/UNSflow.exe new file mode 100644 index 00000000..af9504e2 Binary files /dev/null and b/builddir/UNSflow.exe differ diff --git a/builddir/libamd.dll b/builddir/libamd.dll new file mode 100644 index 00000000..e6074942 Binary files /dev/null and b/builddir/libamd.dll differ diff --git a/builddir/libatomic-1.dll b/builddir/libatomic-1.dll new file mode 100644 index 00000000..c296b794 Binary files /dev/null and b/builddir/libatomic-1.dll differ diff --git a/builddir/libcamd.dll b/builddir/libcamd.dll new file mode 100644 index 00000000..fe6ef424 Binary files /dev/null and b/builddir/libcamd.dll differ diff --git a/builddir/libccalltest.dll b/builddir/libccalltest.dll new file mode 100644 index 00000000..304de233 Binary files /dev/null and b/builddir/libccalltest.dll differ diff --git a/builddir/libccolamd.dll b/builddir/libccolamd.dll new file mode 100644 index 00000000..0983623f Binary files /dev/null and b/builddir/libccolamd.dll differ diff --git a/builddir/libcholmod.dll b/builddir/libcholmod.dll new file mode 100644 index 00000000..5b348b33 Binary files /dev/null and b/builddir/libcholmod.dll differ diff --git a/builddir/libcolamd.dll b/builddir/libcolamd.dll new file mode 100644 index 00000000..9100592a Binary files /dev/null and b/builddir/libcolamd.dll differ diff --git a/builddir/libdSFMT.dll b/builddir/libdSFMT.dll new file mode 100644 index 00000000..28455072 Binary files /dev/null and b/builddir/libdSFMT.dll differ diff --git a/builddir/libexpat-1.dll b/builddir/libexpat-1.dll new file mode 100644 index 00000000..e1a1a311 Binary files /dev/null and b/builddir/libexpat-1.dll differ diff --git a/builddir/libgcc_s_seh-1.dll b/builddir/libgcc_s_seh-1.dll new file mode 100644 index 00000000..518ffc18 Binary files /dev/null and b/builddir/libgcc_s_seh-1.dll differ diff --git a/builddir/libgfortran-3.dll b/builddir/libgfortran-3.dll new file mode 100644 index 00000000..b7670836 Binary files /dev/null and b/builddir/libgfortran-3.dll differ diff --git a/builddir/libgit2.dll b/builddir/libgit2.dll new file mode 100644 index 00000000..1ae7e4b8 Binary files /dev/null and b/builddir/libgit2.dll differ diff --git a/builddir/libgmp.dll b/builddir/libgmp.dll new file mode 100644 index 00000000..9726973f Binary files /dev/null and b/builddir/libgmp.dll differ diff --git a/builddir/libllvmcalltest.dll b/builddir/libllvmcalltest.dll new file mode 100644 index 00000000..bdd1658a Binary files /dev/null and b/builddir/libllvmcalltest.dll differ diff --git a/builddir/libmbedcrypto.dll b/builddir/libmbedcrypto.dll new file mode 100644 index 00000000..6f216825 Binary files /dev/null and b/builddir/libmbedcrypto.dll differ diff --git a/builddir/libmbedtls.dll b/builddir/libmbedtls.dll new file mode 100644 index 00000000..159cc9cb Binary files /dev/null and b/builddir/libmbedtls.dll differ diff --git a/builddir/libmbedx509.dll b/builddir/libmbedx509.dll new file mode 100644 index 00000000..38693f95 Binary files /dev/null and b/builddir/libmbedx509.dll differ diff --git a/builddir/libmpfr.dll b/builddir/libmpfr.dll new file mode 100644 index 00000000..ddd1ef73 Binary files /dev/null and b/builddir/libmpfr.dll differ diff --git a/builddir/libopenlibm.dll b/builddir/libopenlibm.dll new file mode 100644 index 00000000..3e48588c Binary files /dev/null and b/builddir/libopenlibm.dll differ diff --git a/builddir/libpcre2-8.dll b/builddir/libpcre2-8.dll new file mode 100644 index 00000000..27752917 Binary files /dev/null and b/builddir/libpcre2-8.dll differ diff --git a/builddir/libpcre2-posix.dll b/builddir/libpcre2-posix.dll new file mode 100644 index 00000000..86e6c276 Binary files /dev/null and b/builddir/libpcre2-posix.dll differ diff --git a/builddir/libquadmath-0.dll b/builddir/libquadmath-0.dll new file mode 100644 index 00000000..0325cecd Binary files /dev/null and b/builddir/libquadmath-0.dll differ diff --git a/builddir/libspqr.dll b/builddir/libspqr.dll new file mode 100644 index 00000000..f3aa3be9 Binary files /dev/null and b/builddir/libspqr.dll differ diff --git a/builddir/libssh2.dll b/builddir/libssh2.dll new file mode 100644 index 00000000..041637c7 Binary files /dev/null and b/builddir/libssh2.dll differ diff --git a/builddir/libssp-0.dll b/builddir/libssp-0.dll new file mode 100644 index 00000000..d3c4c314 Binary files /dev/null and b/builddir/libssp-0.dll differ diff --git a/builddir/libstdc++-6.dll b/builddir/libstdc++-6.dll new file mode 100644 index 00000000..73482372 Binary files /dev/null and b/builddir/libstdc++-6.dll differ diff --git a/builddir/libsuitesparse_wrapper.dll b/builddir/libsuitesparse_wrapper.dll new file mode 100644 index 00000000..07a55046 Binary files /dev/null and b/builddir/libsuitesparse_wrapper.dll differ diff --git a/builddir/libsuitesparseconfig.dll b/builddir/libsuitesparseconfig.dll new file mode 100644 index 00000000..2887f4b4 Binary files /dev/null and b/builddir/libsuitesparseconfig.dll differ diff --git a/builddir/libumfpack.dll b/builddir/libumfpack.dll new file mode 100644 index 00000000..a38c9e4f Binary files /dev/null and b/builddir/libumfpack.dll differ diff --git a/builddir/libwinpthread-1.dll b/builddir/libwinpthread-1.dll new file mode 100644 index 00000000..893d1be4 Binary files /dev/null and b/builddir/libwinpthread-1.dll differ diff --git a/builddir/zlib1.dll b/builddir/zlib1.dll new file mode 100644 index 00000000..08c6becb Binary files /dev/null and b/builddir/zlib1.dll differ diff --git a/main.jl b/main.jl new file mode 100644 index 00000000..abfede4b --- /dev/null +++ b/main.jl @@ -0,0 +1,98 @@ +include("src/UnsteadyFlowSolvers.jl") +## Define Geometry and Kinematics +# Kinematics +case = 3 # kinematics case to run +if case == 1 # ~ 8 hours to run, Ramshes @ 4 hrs + amp = 35 #degrees + k = .005 + tstart = 60 # non-dimensional time + t_tot = 277 # 138.525 for up only, 277 for full +elseif case == 2 + amp = 45 + k = .05 + tstart = 20 + t_tot = 60.15 # 29.5 for up only, 60.15 for full +elseif case == 3 + amp = 90 + k = .4 + tstart = 10 + t_tot = 25 # 12.5 for up only, 25 for full +elseif case == 4 + amp = 25 + k = .11 + a = 11 + tstart = 1 + t_tot = 7.09 # 3.55 for up only, 7.09 for full +elseif case == 5 + amp = 45 + k = 0.4 + a = 11 + tstart = 1 + t_tot = 4.52 # 2.26 for up only, 4.52 for full +elseif case == 6 + mean = 0 + amp = 30 + k = 0.1 + phi = 0 + t_tot = 31.4 #roughly +elseif case == 7 + mean = 0 + amp = 30 + k = 0.4 + phi = 0 + t_tot = 7.85 +end +#amp = amp * pi / 180 +if case < 4 +a = (pi^2 * k*180 )/(2*amp*pi *(1-0.1)) +end + +#alpha = 10 +if case < 6 + alphadef = UnsteadyFlowSolvers.EldRampReturntstartDef(amp*pi/180,k,a,tstart) # ConstDef(alpha*pi/180)# +else + alphadef = UnsteadyFlowSolvers.SinDef(mean,amp*pi/180,k,phi) +end +hdef = UnsteadyFlowSolvers.ConstDef(0) +udef = UnsteadyFlowSolvers.ConstDef(1) +full_kinem = UnsteadyFlowSolvers.KinemDef(alphadef, hdef, udef) +# Geometry +pvt = 0.25 ; +geometry = "bin/sd7003.dat"#"FlatPlate"# +lespcrit = [.2;] + +surf = UnsteadyFlowSolvers.TwoDSurf(geometry,pvt,full_kinem,lespcrit; ndiv = 101, camberType = "linear", rho = 1.225) +curfield = UnsteadyFlowSolvers.TwoDFlowField() +# Iteration Parameters +dtstar = .015 #UnsteadyFlowSolvers.find_tstep(alphadef) + +nsteps = Int(round(t_tot/dtstar)) +#nsteps = 1 #DEBUG + +println("Running LVE") +frames, IFR_field, mat, test1, test2 = UnsteadyFlowSolvers.LVE(surf,curfield,nsteps,dtstar,40,"longwake") +writedlm("myMat.txt",mat) + +# Velocity plot +using Plots +pyplot() +single = "false" ; +if single == "none" + mesh = frames[end] + + contour(mesh.x[1,:],mesh.z[:,1], mesh.velMag,fill=true,levels = 200,c = :lightrainbow_r) + scatter!(mesh.vorX, mesh.vorZ,markersize = 1.5, markerstrokestyle = :dot, markerstrokecolor = :white) + plot!(mesh.camX,mesh.camZ,color = :black, linewidth = 3, aspect_ratio = :equal,legend = :none,axis = false) +elseif single == "false" # Make full list of images in folder + dir = "case_$case steps_$nsteps" + mkdir(dir) + for i = 1:length(frames) + mesh = frames[i] + + contour(mesh.x[1,:],mesh.z[:,1], mesh.velMag,fill=true,levels = 200,c = :lightrainbow_r) + scatter!(mesh.vorX, mesh.vorZ,markersize = 1.5, markerstrokestyle = :dot, markerstrokecolor = :white) + plot!(mesh.camX,mesh.camZ,color = :black, linewidth = 3, aspect_ratio = :equal,legend = :none,axis = false) + + savefig("$dir/$i") + end +end diff --git a/src/UI/UI_main.jl b/src/UI/UI_main.jl new file mode 100644 index 00000000..c6bec7c8 --- /dev/null +++ b/src/UI/UI_main.jl @@ -0,0 +1,310 @@ +#= +UI_main.jl + Written by Matthew White + 11/14/2019 +=# +function runUI() + println("\33[2J") # clear command window + commandHist = [] # User command history log + pdef = ConstDef(0) # pre-define kinematics + udef = ConstDef(1) + alphadef = ConstDef(0) + alphadef = ConstDef(0) + curfield = UnsteadyFlowSolvers.TwoDFlowField() + delvort = UnsteadyFlowSolvers.delNone() + nsteps = 0 # pre-define global vars + trange = collect(1:10) + lockTime = false + mat = [] + frames = [] + + ## Dictonary of Cmdmands + # Stage 1 commands + list1 = ["airfoil","ndiv","lespcrit","alphadef","oper","pvt","runtime","pdef","veldef","dt","rho","cambercalc","motionfile"] + help1 = ["Define airfoil designation: coordinate file or FlatPlate","Define airfoil panel division count","Define critical LESP value for LEV shedding","Define pitch motion parameters","Change to operation phase","Define airfoil pivot location","Define simulation run time","Define plunge motion parameters","Define velocity profile parameters","Define time step","Define operating density","Define camber calculation method","Use a file to describe the motion"] + help1 = hcat(list1,help1) + help1[1,1] = "load" # overwrite to show command + help1 = sortslices(help1,dims = 1) + + # Stage 2 commands + list2 = ["wflag","sflag","wcount","run","save","anim","animstep","plot","vortplot","infoplot"] + help2 = ["If enabled, writes out iteration files to default save location","If enabled, start simulation from previous output files","Define number of output files","Run a simulation","Save the polar data to a file","If enabled, shows an animation of polars and motion during LVE simulation","Define step count between animation frames","Plot simulation output data","Plot vortex location data for each write step (Must have wflag enabled)","Plot step parameters for each write step (Must have wflag enabled)"] + help2 = sortslices(hcat(list2,help2),dims = 1) + + commandDict = commands(list1,list2) + commandDict = setDefaults(commandDict) + + stage = 1 # set preliminary stage + + # Run 'HELP' on startup + helpMe(stage,help1,help2) + + while true # Loop forever or unitl user exit + if stage == 1 # Set stage prompts + prompt = "UNSflow c> " + cmds = commandDict.s1 + elseif stage == 2 + prompt = "OPER c> " + cmds = commandDict.s2 + end + ## User Input + print(prompt) + usrCmd = ask() + commandHist = vcat(commandHist,usrCmd) + + ## check for multiple commands in one line + if occursin(" ",usrCmd) && length(usrCmd) > 2 # sub command + temp = split(usrCmd," ") + if length(temp) <= 2 # extra command added + usrCmd = lowercase(temp[1]) + subCmd = temp[2] + else + errorHandler(usrCmd,"many") + continue + end + else # no sub command + usrCmd = lowercase(usrCmd) + subCmd = "" + end + + ## Top level Commands + if stage == 1 && usrCmd == "exit" + close("all") + dirvec = readdir() # Clear step files if prompted + if "Step Files" in dirvec + if ynAsk("Clear step file cache?") + cleanWrite() + end + end + break # exit program + end + if stage == 2 && (usrCmd == "" || usrCmd == "back") + stage = 1 + close("all") + continue # return to stage 1 + end + if usrCmd == "?" + cmdValues(commandDict,stage) + continue # show variables + end + if usrCmd == "clc" + println("\33[2J") # clear command window + continue # clear screen + end + if usrCmd == "help" + helpMe(stage,help1,help2) + continue + end + + if stage == 1 && usrCmd == "load" # change from load command + usrCmd = "airfoil" + end + if stage == 2 && (usrCmd == "run") # Calculate surf and curfield + global surf = UnsteadyFlowSolvers.TwoDSurf(commandDict.s1["airfoil"],commandDict.s1["pvt"],kinem,[commandDict.s1["lespcrit"];]; ndiv = commandDict.s1["ndiv"], camberType = commandDict.s1["cambercalc"], rho = commandDict.s1["rho"]) + curfield = UnsteadyFlowSolvers.TwoDFlowField() + end + + ## General Commands + if haskey(cmds,usrCmd) # Test if command is available + ## Evaluate Cmdmands + if stage == 1 # Initilization Stage + if usrCmd == "airfoil" # load airfoil + o = loadGeo(subCmd) + if o != nothing + commandDict.s1[usrCmd] = o + makeInitPlot(commandDict.s1,trange,alphadef,pdef,udef) + end + + elseif usrCmd == "ndiv" # Set panel divsions + o = svPrompt("i", "Panel Count: ", subCmd) + if o != nothing + commandDict.s1[usrCmd] = o + makeInitPlot(commandDict.s1,trange,alphadef,pdef,udef) + end + + elseif usrCmd == "lespcrit" # Set panel divsions + o = svPrompt("f", "LESP Critical: ", subCmd) + if o != nothing ; commandDict.s1[usrCmd] = o end + + elseif usrCmd == "pvt" # Set pivot location + o = svPrompt("f", "Pivot Point: ", subCmd) + if o != nothing ; commandDict.s1[usrCmd] = o end + + elseif usrCmd == "dt" # Set time step + if !lockTime + o = svPrompt("f", "Time Step: ", subCmd) + if o != nothing + commandDict.s1[usrCmd] = o + if commandDict.s1["runtime"] != "---" + nsteps, trange = iterCalc(commandDict.s1["runtime"],o) + println(" The simulation will take $nsteps iterations.") + end + end + else + errorHandler(1,"The file '$(commandDict.s1["motionfile"])' prevents this from being changed.") + end + + elseif usrCmd == "runtime" # Set total time + if !lockTime + o = svPrompt("f", "Total Time: ", subCmd) + if o != nothing + commandDict.s1[usrCmd] = o + nsteps, trange = iterCalc(o,commandDict.s1["dt"]) + println(" The simulation will take $nsteps iterations.") + end + else + errorHandler(1,"The file '$(commandDict.s1["motionfile"])' prevents this from being changed.") + end + + elseif usrCmd == "rho" # Set density + o = svPrompt("f", "Density [kg/m^3]: ", subCmd) + if o != nothing ; commandDict.s1[usrCmd] = o end + + elseif usrCmd == "cambercalc" # Set camber calculation method + o = setCamberCalc(subCmd) + if o != nothing + commandDict.s1[usrCmd] = o + makeInitPlot(commandDict.s1,trange,alphadef,pdef,udef) + end + + elseif usrCmd == "motionfile" # load time, alphadef, h, and vel from file + t,a,h,u,fn = fileInput(subCmd) + if t != nothing + commandDict.s1["motionfile"] = fn + + alphadef = FileDef(t,a) + commandDict.s1["alphadef"] = "file" + + pdef = FileDef(t,h) + commandDict.s1["pdef"] = "file" + + udef = FileDef(t,u) + commandDict.s1["veldef"] = "file" + + commandDict.s1["dt"] = t[2]-t[1] + commandDict.s1["runtime"] = t[end] + lockTime = true # might need a way to toggle off if file is completely overrided + + trange = t + makeInitPlot(commandDict.s1,trange,alphadef,pdef,udef) + end + + elseif usrCmd == "alphadef" # Set alpha definition + alphavec = modef(subCmd) + if alphavec != nothing + alphadef = alphaCorr(alphavec) + param = join(alphavec[2:end], ", ") + commandDict.s1[usrCmd] = "$(alphavec[1])($param)" + rt = endTime(alphadef) + if rt != "---" + if commandDict.s1["runtime"] == "---" || rt > commandDict.s1["runtime"] && !lockTime + commandDict.s1["runtime"] = rt + nsteps, trange = iterCalc(commandDict.s1["runtime"],commandDict.s1["dt"]) + end + end + makeInitPlot(commandDict.s1,trange,alphadef,pdef,udef) + end + + elseif usrCmd == "pdef" || usrCmd == "veldef"# Set height and velocity definition + o = modef(subCmd) + if o != nothing + if o[1] == "const" + odef = ConstDef(o[2]) + elseif o[1] == "lin" + odef = LinearDef(o[2],o[3],o[4],o[5]) + elseif o[1] == "eld" + odef = EldRampReturntstartDef(o[2],o[3],o[4],o[5]) + elseif o[1] == "sin" + odef = SinDef(o[2],o[3],o[4],o[5]) + elseif o[1] == "cos" + odef = CosDef(o[2],o[3],o[4],o[5]) + elseif o[1] == "gust" + odef = StepGustDef(o[2],o[3],o[4]) + end + param = join(o[2:end], ", ") + commandDict.s1[usrCmd] = "$(o[1])($param)" + + rt = endTime(odef) + + if rt != "---" + if commandDict.s1["runtime"] == "---" || rt > commandDict.s1["runtime"] && !lockTime + commandDict.s1["runtime"] = rt + nsteps, trange = iterCalc(commandDict.s1["runtime"],commandDict.s1["dt"]) + end + end + + if usrCmd == "pdef" + pdef = odef + else + udef = odef + end + makeInitPlot(commandDict.s1,trange,alphadef,pdef,udef) + end + + elseif usrCmd == "oper" # Change to stage 2 + o = gotoOper(commandDict.s1) + if o + stage = 2 + global kinem = KinemDef(alphadef, pdef, udef) + nsteps = Int(ceil(commandDict.s1["runtime"]/commandDict.s1["dt"])) + end + close() + end + + elseif stage == 2 # Operation Stage + if usrCmd == "wcount" # Set step output file count + o = svPrompt("i","Step output file count: ", subCmd) + if o != nothing ; commandDict.s2[usrCmd] = o end + + elseif usrCmd == "wflag" # Write step output during sim + commandDict.s2[usrCmd] = ynAsk("Write step output?",subCmd) + + elseif usrCmd == "sflag" # Start from previous step output + commandDict.s2[usrCmd] = ynAsk("Start from previous step output?",subCmd) + + elseif usrCmd == "anim" # Set LVE animation state + commandDict.s2[usrCmd] = ynAsk("Display LVE Animation?", subCmd) + + elseif usrCmd == "animstep" # set animation step + o = svPrompt("i", "Animation Step: ", subCmd) + if o != nothing ; commandDict.s2[usrCmd] = o end + + elseif usrCmd == "run" # Run simulation + if commandDict.s2["wflag"] + cleanWrite() + end + o = runSim(commandDict,surf,curfield,nsteps,delvort,subCmd) + if o != nothing ; mat = o end + + elseif usrCmd == "plot" + plotResults(mat,subCmd) + + elseif usrCmd == "save" # Set step output file count + if mat != [] # mat has been defined + saveMat(subCmd,mat) + else + errorHandler(mat,"A simulation must be ran before saving the results.") + end + + elseif usrCmd == "infoplot" + if mat != [] + tt = plotTime(subCmd) + makeKinemClVortPlots2D(mat,tt) + else + errorHandler(1,"A simulation must be ran before infoplots can be generated.") + end + + elseif usrCmd == "vortplot" + tt = plotTime(subCmd) + makeVortPlots2D(tt) + + end + end + elseif usrCmd != "" # return error message + errorHandler(usrCmd) + continue + end + + end + #return commandDict +end diff --git a/src/UI/UI_utils.jl b/src/UI/UI_utils.jl new file mode 100644 index 00000000..f972a8dc --- /dev/null +++ b/src/UI/UI_utils.jl @@ -0,0 +1,271 @@ +#= +UI_utils.jl + Written by Matthew White + 11/14/2019 +=# +struct commands + s1::DataStructures.SortedDict + s2::DataStructures.SortedDict + + function commands(list1,list2) + # Stage 1 + s1 = DataStructures.SortedDict{String,Any}(list1[i] => "---" for i = 1:length(list1)) + + #Stage 2 + s2 = DataStructures.SortedDict{String,Any}(list2[i] => "---" for i = 1:length(list2)) + + new(s1,s2) + end +end + +function setDefaults(commandDict) + commandDict.s1["ndiv"] = 100 + commandDict.s1["pvt"] = 0.25 + commandDict.s1["lespcrit"] = 10000 + commandDict.s1["dt"] = .015 + commandDict.s1["veldef"] = "const(1)" + commandDict.s1["pdef"] = "const(0)" + commandDict.s1["rho"] = 1.225 + commandDict.s1["cambercalc"] = "radial" + commandDict.s1["oper"] = "N/A" + + commandDict.s2["wflag"] = false + commandDict.s2["sflag"] = false + commandDict.s2["wcount"] = 20 + commandDict.s2["anim"] = true + commandDict.s2["animstep"] = 50 + + return commandDict +end + +function helpMe(stage,help1,help2) + # Display help data based on stage + if stage == 1 # set relevant help results + help = help1 + elseif stage == 2 + help = help2 + end + c_offset = 12 + println() + println(" Available Commands:") + for i = 1:size(help,1) # print values + print("\t$(help[i,1])") + for j = 1:c_offset - length(help[i,1]) + print(" ") + end + print(":\t") + println("$(help[i,2])") + end + println() +end + +function ask(prompt::String = "") + # Ask a question + print(" $prompt") + usrCmd = strip(chomp(readline())) +end + +function cmdValues(commandDict,stage) + # display the values of all relevant variables + c_offset = 12 + + println() + if stage == 1 #decide paramters + dict = commandDict.s1 + exclude = ["oper"] + elseif stage == 2 + dict = commandDict.s2 + exclude = ["save","run","plot","vortplot","infoplot"] + end + + for entry in dict # print values + if sum(entry[1] .== exclude) == 0 # not on exclude list + print("\t$(entry[1])") + for i = 1:c_offset - length(entry[1]) + print(" ") + end + print(":\t") + println("$(entry[2])") + end + end + println() +end + +function errorHandler(usrCmd,type::String = "") + # Write out error messages depending on type of input expected + # "" if basic error message + + print(" ERROR: ") + if length(type) > 10 # custom error message + println(type) + else + if type == "" + println("'$usrCmd' not recognized. Type HELP for a list of available commands.") + elseif type == "many" + println("Too many arguments have been entered.") + elseif type == "valid" + println("'$usrCmd' is not a vaild input.") + elseif type == "list" + prompt = join(usrCmd[1:end-1], ", ") + println("Simulation options are $prompt and $(usrCmd[end]).") + else + print("'$usrCmd' must be ") + if type == "i" + println("an integer.") + elseif type == "f" + println("a float.") + elseif type == "s" + println("a string") + elseif tpye == "def" + println("either 'eld', 'const', 'sin', or 'cos'") + end + end + end +end + +function svPrompt(type, prompt, subCmd = "") + # Ask for single output variable + if type == "i" + fullType = Int64 + elseif type == "f" + fullType = Float64 + end + while true + if subCmd == "" + subCmd = ask(prompt) + if lowercase(subCmd) == "exit" || lowercase(subCmd) == "" + break + end + end + if type != "s" + outVar = tryparse(fullType,subCmd) + else + outVar = subCmd + end + if outVar != nothing + return outVar + else + errorHandler(subCmd,type) + end + subCmd = "" + end +end + +function ynAsk(prompt,state = "") + # Set a value to true or false + state = lowercase(state) + + if state == "on" + return true + elseif state == "off" + return false + elseif state != "" + errorHandler(state,"valid") + end + + while true + prompt = "$prompt [Y/N]: " + usrIn = lowercase(ask(prompt)) + if usrIn == "y" + return true + elseif usrIn == "n" + return false + else + errorHandler(usrIn,"Must be 'Y' or 'N'") + end + end +end + +function endTime(kindef) + # Returns the total time for a kinematic motion to run + + if typeof(kindef) == SinDef || typeof(kindef) == CosDef + rt = pi / kindef.k + + elseif typeof(kindef) == EldRampReturntstartDef + fr = kindef.K/(pi*abs(kindef.amp)); + t1 = kindef.tstart + t2 = t1 + (1. /(2*pi*fr)); + t3 = t2 + ((1/(4*fr)) - (1/(2*pi*fr))); + t4 = t3 + (1. /(2*pi*fr)); + rt = t4+t1; + + elseif typeof(kindef) == LinearDef + rt = kindef.tstart + kindef.len + + else + rt = "---" + end + return rt +end + +function iterCalc(rt,dt) + # Calculate iteration count and time vector + iter = Int(ceil(rt/dt)) + trange = collect(0:dt:iter*dt) + return iter, trange +end + +function kinrange(kindef,trange) + # calculate kinematic output at each time step + kinem = [kindef(i) for i in trange] + return kinem +end + +function airfoilGeometry(coord_file,ndiv,camberType) + # Extract geometry coordinates and calculate camber line + + theta = zeros(ndiv) + x = zeros(ndiv) + cam = zeros(ndiv) + cam_slope = zeros(ndiv) + + if (coord_file != "FlatPlate") + A = DelimitedFiles.readdlm(coord_file, Float64); # actual geometry + c = maximum(A[:,1]) - minimum(A[:,1]) + else + c = 1 + end + + + if camberType == "radial" + dtheta = pi/(ndiv-1) + for ib = 1:ndiv + theta[ib] = (ib-1.)*dtheta + x[ib] = c/2. *(1-cos(theta[ib])) + end + elseif camberType == "linear" + dx = c / (ndiv-1) + for ib = 2:ndiv + x[ib] = x[ib-1] + dx + end + end + + if (coord_file != "FlatPlate") + cam, cam_slope = camber_calc(x, coord_file) # camber line + else + A = hcat(x,cam) + end + C = hcat(x,cam) + return A,C +end + +function makeInitPlot(cmds,trange,alphadef,pdef,udef) + # Create all inputs for intitilization plot and plot results + + ## Airfoil + if cmds["airfoil"] != "---" + A, C = airfoilGeometry(cmds["airfoil"],cmds["ndiv"],cmds["cambercalc"]) + else + A = 0 + C = 0 + end + ## Kinem ranges + arange = kinrange(alphadef,trange) .* 180 ./ pi + hrange = kinrange(pdef,trange) + urange = kinrange(udef,trange) + + ## Plot + clf() + initPlot(A,C,trange,arange,hrange,urange) +end diff --git a/src/UI/s1Cmds.jl b/src/UI/s1Cmds.jl new file mode 100644 index 00000000..e4033101 --- /dev/null +++ b/src/UI/s1Cmds.jl @@ -0,0 +1,205 @@ +#= +s1Cmds.jl + Written by Matthew White + 11/14/2019 +=# +function loadGeo(geo = "") + # Assign airfoil file data + while true + if geo == "" # if no secondary command given + geo = ask("Enter an airfoil geometry file: ") + if lowercase(geo) == "exit" || lowercase(geo) == "" + break + end + end + if lowercase(geo) == "flatplate" + return "FlatPlate" + end + if geo[1] == '/' # remove any beginning slash + geo = geo[2:end] + end + if occursin("/",geo) # Check if file exsists + path = split(geo,"/") + fn = path[end] + path =join(path[1:end-1],"/") + try + files = cd(readdir, path) + catch + errorHandler(path,"'$path' is not a valid directory.") + geo = "" + continue + end + else + files = readdir() + fn = geo + end + if sum(files .== fn) == 0 # file does not exist in directory + errorHandler(fn,"File '$fn' does not exist.") + geo = "" + continue + end + return geo + end +end + +function fileInput(name = "") + exts = ["dat","txt"] # valid file extensions + while true + if name == "" + name = ask("File Name: ") + end + if lowercase(name) == "exit" || lowercase(name) == "" + break + end + parts = split(name,".") + if occursin("/",name) # can change in the future to look for available folders + errorHandler(name,"$name must not be inside a folder.") + + elseif length(parts) > 2 # extra periods + errorHandler(name,"$name is an invalid file name.") + + elseif sum(parts[end] .== exts) == 0 # not a valid extention + errorHandler(name,"$(parts[end]) is an invalid extension.") + + else + # potentially change for asking for which column is which + mat = DelimitedFiles.readdlm(name) + t = mat[:,1] + a = mat[:,2].*pi/180 + h = mat[:,3] + u = mat[:,4] + return t,a,h,u,name + end + name = "" + end +end + +function modef(motype = "") + # Create alphadef/pdef/udef motion vectors + + modefs = ["const","lin","eld","sin","cos"]#,"gust"] + outVar = [] + + while true # Motion type + if motype == "" # if no secondary command given + motype = lowercase(ask("Motion Type: ")) + if motype == "exit" || motype == "" + return nothing + end + end + if sum(motype .== modefs) > 0 # motion type matches + outVar = [outVar;motype] + break # next question + else + # Error message + errorHandler(modefs,"list") + end + motype = "" + end + + if motype == "const" + amp = svPrompt("f","Magnitude: ") + if amp == nothing ; return nothing end + outVar = [outVar;amp] + elseif motype == "lin" + sm = svPrompt("f","Start Magnitude: ") + if sm == nothing ; return nothing end + fm = svPrompt("f","Final Magnitude: ") + if fm == nothing ; return nothing end + l = svPrompt("f","Rise Time: ") + if l == nothing ; return nothing end + ynAsk("Change start time?") ? tstart = svPrompt("f","Start Time: ") : tstart = 1 + if tstart == nothing ; return nothing end + outVar = [outVar;tstart;sm;fm;l] + elseif motype == "eld" + amp = svPrompt("f","Amplitude: ") + if amp == nothing ; return nothing end + K = svPrompt("f","K: ") + if K == nothing ; return nothing end + ynAsk("Change a?") ? a = svPrompt("f","a: ") : a = (pi^2 * K*180 )/(2*amp*pi *(1-0.1)) + if a == nothing ; return nothing end + ynAsk("Change start time?") ? tstart = svPrompt("f","Start Time: ") : tstart = 1 + if tstart == nothing ; return nothing end + outVar = [outVar;amp;K;a;tstart] + elseif motype == "sin" || motype == "cos" + mean = svPrompt("f","Mean: ") + if mean == nothing ; return nothing end + amp = svPrompt("f","Amplitude: ") + if amp == nothing ; return nothing end + K = svPrompt("f","K: ") + if K == nothing ; return nothing end + phi = svPrompt("f","Phi: ") + if phi == nothing ; return nothing end + outVar = [outVar;mean;amp;K;phi] + elseif motype == "gust" + amp = svPrompt("f","Magnitude: ") + if amp == nothing ; return nothing end + ynAsk("Change start time?") ? tstart = svPrompt("f","Start Time: ") : tstart = 1 + if tstart == nothing ; return nothing end + tgust = svPrompt("f","Gust Start Time: ") + if tgust == nothing ; return nothing end + outVar = [outVar;amp;tstart;tgust] + end + return outVar +end + +function setCamberCalc(meth) + methods = ["radial","linear"] + eStr = join(methods[1:end-1], ",") + eStr = "Accepts either $eStr, or $(methods[end])" + while true + if meth == "" + prompt = "Camber calculation method: " + meth = lowercase(ask(prompt)) + end + if sum(meth .== methods) == 0 # not a method + errorHandler(meth,eStr) + else + return meth + end + meth = "" + end +end + +function alphaCorr(alphavec) + # Turn alpha into radians + alphavec[2] = alphavec[2]*pi/180 + + if alphavec[1] == "const" + alphadef = ConstDef(alphavec[2]) + elseif alphavec[1] == "lin" + alphavec[3:4] = alphavec[3:4].*pi./180 + alphadef = LinearDef(alphavec[2],alphavec[3],alphavec[4],alphavec[5]) + elseif alphavec[1] == "eld" + alphadef = EldRampReturntstartDef(alphavec[2],alphavec[3],alphavec[4],alphavec[5]) + elseif alphavec[1] == "sin" + alphavec[3] = alphavec[3]*pi/180 + alphavec[5] = alphavec[5]*pi/180 + alphadef = SinDef(alphavec[2],alphavec[3],alphavec[4],alphavec[5]) + elseif alphavec[1] == "cos" + alphavec[3] = alphavec[3]*pi/180 + alphavec[5] = alphavec[5]*pi/180 + alphadef = CosDef(alphavec[2],alphavec[3],alphavec[4],alphavec[5]) + elseif alphavec[1] == "gust" + alphadef = StepGustDef(alphavec[2],alphavec[3],alphavec[4]) + end + return alphadef +end + +function gotoOper(cmds) + # Cheack whether there are enough inputs to move to operation stage (2) + o = false + notdef = [k for (k,v) in cmds if v=="---" && k != "motionfile"] + if length(notdef) > 1 # incomplete intitilization + prompt = join(notdef[1:end-1], ", ") + prompt = "$prompt, and $(notdef[end]) have not been defined." + errorHandler(o,prompt) + elseif length(notdef) == 1 + prompt = "$(notdef[1]) has not been defined." + errorHandler(o,prompt) + else + println(" The simulation will take $(Int(ceil(cmds["runtime"]/cmds["dt"]))) iterations.") + o = true + end + return o +end diff --git a/src/UI/s2Cmds.jl b/src/UI/s2Cmds.jl new file mode 100644 index 00000000..a6ab23bf --- /dev/null +++ b/src/UI/s2Cmds.jl @@ -0,0 +1,158 @@ +#= +s2Cmds.jl + Written by Matthew White + 11/16/2019 +=# +using DelimitedFiles + +function saveMat(name,mat) + # Save the output matrix + exts = ["dat","txt"] # valid file extensions + while true + if name == "" + name = ask("File Name: ") + end + if lowercase(name) == "exit" || lowercase(name) == "" + break + end + parts = split(name,".") + if occursin("/",name) # can change in the future to look for available folders + errorHandler(name,"$name must not be inside a folder.") + + elseif length(parts) > 2 # extra periods + errorHandler(name,"$name is an invalid file name.") + + elseif sum(parts[end] .== exts) == 0 # not a valid extention + errorHandler(name,"$(parts[end]) is an invalid extension.") + + else + writedlm(name,mat) + break + end + name = "" + end +end + +function runSim(commandDict,surf,curfield,nsteps,delvort,sim) + # Run the specificed simulation, returns mat + sim = lowercase(sim) + simList = ["lve","ldvm"] + while true + if sim == "lve" + println() + println(" Running LVE Simulation...") + writeInterval = floor(nsteps/commandDict.s2["wcount"]) + frames, IFR_field, mat = LVE(surf,curfield,nsteps,commandDict.s1["dt"],commandDict.s2["wcount"],"UI Window",1000,commandDict.s2["anim"],30,commandDict.s2["wflag"],writeInterval) + println("\t Done!") + println() + return mat + elseif sim == "ldvm" + println() + println(" Running LDVM Simulation...") + writeInterval = commandDict.s1["runtime"]/commandDict.s2["wcount"] + commandDict.s2["wflag"] ? wflag = 1 : wflag = 0 + commandDict.s2["sflag"] ? sflag = 1 : sflag = 0 + mat, surf, curfield = ldvm(surf, curfield, nsteps, commandDict.s1["dt"],sflag, wflag, writeInterval, delvort) + println("\t Done!") + println() + return mat + end + + sim = ask("Simulation Type: ") + if sim == "" + break + elseif sum(simList .== sim) == 0 # if not on simulation list + errorHandler(simList,"lsit") + end + end +end + +function plotResults(mat,subCmd) + # Plot anything from sim results + ploty = [] + plotCnt = [] + plotx = [] + matCol = ["t","alpha","h","u","lesp","cl","cd","cm"] + matIdx = DataStructures.SortedDict{String,Any}(matCol[i] => i for i = 1:length(matCol)) + if haskey(matIdx,lowercase(subCmd)) # subcmd is plot type (infers only one plot to be calculated) + plotCnt = 1 + ploty = [subCmd] + else + plotCnt = tryparse(Int64,subCmd) + if plotCnt == nothing # subcmd is number of plots + plotCnt = [] + end + end + # Gather plot info + while true + if plotCnt == [] # number of plots + plotCnt = svPrompt("i","Number of plots: ") + if plotCnt == nothing + return nothing + end + end + if plotCnt > 7 # too many plots + plotCnt = [] + errorHandler(plotCnt,"$plotCnt exceeds the maximum number of plots.") + continue + else + break + end + end + if ploty == [] # determine y axis(es) + for i = 1:plotCnt + while true + ploty = [ploty;svPrompt("s","Plot $i y-axis: ")] + if ploty[i] == nothing + return nothing + end + if !haskey(matIdx,lowercase(ploty[i])) # not a vaild y axis + errorHandler(matCol,"list") + ploty = ploty[1:end-1] # reset ploty + continue + else + break + end + end + end + end + while true # determine x axis + plotx = svPrompt("s","X axis: ") + if plotx == nothing + return nothing + end + if !haskey(matIdx,lowercase(plotx)) # not a vaild x axis + errorHandler(matCol,"list") + continue + else + break + end + end + plotMat(mat,plotCnt,plotx,ploty,matIdx) + + if ynAsk("Save figure?") + fn = svPrompt("s","File name: ") + if occursin(".",fn) + temp = split(fn,".") # pyplot only supports .png + fn = temp[1] + end + fn = "$fn.png" + savefig(fn) + end +end + +function plotTime(tt = "") + # Asks for a specific time to plot step file results + ttparse = tryparse(Float64, tt) + if tt == "all" + return tt + elseif tt == "" || ttparse == nothing + if ynAsk("Plot all step files?") + return "all" + else + return svPrompt("f","Time to plot: ") + end + else + return ttparse + end +end diff --git a/src/UnsteadyFlowSolvers.jl b/src/UnsteadyFlowSolvers.jl index c0328e51..723f6f10 100755 --- a/src/UnsteadyFlowSolvers.jl +++ b/src/UnsteadyFlowSolvers.jl @@ -14,16 +14,19 @@ import NLsolve: nlsolve, not_in_place import Statistics: mean +import LinearAlgebra + +import Dates + import PyPlot: plot, scatter, figure, xlabel, ylabel, xlim, ylim, xticks, yticks, subplot, subplot2grid, legend, axis, savefig, - close, tight_layout + close, tight_layout, clf, title, show, pause import Plots: @layout import LaTeXStrings: @L_str -#For use in development and debugging -import Revise +import DataStructures export # kinematics types and funtions @@ -65,10 +68,18 @@ export ldvm, ldvmLin, ldvm2DOF, + LVE, # 2D plot output functions makeForcePlots2D, - makeVortPlots2D + makeVortPlots2D, + subPlot, + + # XFOIL Wrapper + xfoilWrapper, + + # User Interface + runUI ### source files @@ -85,9 +96,19 @@ include("delVort.jl") include("lowOrder2D/typedefs.jl") # type definitions include("lowOrder2D/calcs.jl") # calculation functions include("lowOrder2D/solvers.jl") # solver methods -include("lowOrder2D/postprocess.jl") # postprocessing functions +include("lowOrder2D/LVE.jl") # Lumped Vortex Element method +include("lowOrder2D/postprocess.jl") # 2D plotting functions include("plots/plots2D.jl") +# XFOIL Wrapper +include("xfoil/XfoilWrapper.jl") + +# User Interface +include("UI/UI_main.jl") +include("UI/UI_utils.jl") +include("UI/s1Cmds.jl") +include("UI/s2Cmds.jl") + end diff --git a/src/kinem.jl b/src/kinem.jl index b80fb4a3..ea8da5b6 100755 --- a/src/kinem.jl +++ b/src/kinem.jl @@ -74,6 +74,36 @@ function (eld::EldRampReturnDef)(tt) end +struct EldRampReturntstartDef <: MotionDef + amp :: Float64 + K :: Float64 + a :: Float64 + tstart :: Float64 +end + +function (eld::EldRampReturntstartDef)(tt) + fr = eld.K/(pi*abs(eld.amp)); + t1 = eld.tstart + t2 = t1 + (1. /(2*pi*fr)); + t3 = t2 + ((1/(4*fr)) - (1/(2*pi*fr))); + t4 = t3 + (1. /(2*pi*fr)); + t5 = t4+t1; + + nstep = round(Int,t5/0.015) + 1 + g = zeros(nstep) + t = zeros(nstep) + + for i = 1:nstep + t[i] = (i-1.)*0.015 + g[i] = log((cosh(eld.a*(t[i] - t1))*cosh(eld.a*(t[i] - t4)))/(cosh(eld.a*(t[i] - t2))*cosh(eld.a*(t[i] - t3)))) + end + maxg = maximum(g); + + gg = log((cosh(eld.a*(tt - t1))*cosh(eld.a*(tt - t4)))/(cosh(eld.a*(tt - t2))*cosh(eld.a*(tt - t3)))) + + return eld.amp*gg/(maxg); + +end struct ConstDef <: MotionDef amp :: Float64 @@ -217,3 +247,16 @@ function (eld::EldUpInttstartDef)(t) end amp end + +# Takes a vector of time-ordered data from an input file +struct FileDef <: MotionDef + tvec :: Vector{Float64} + var :: Vector{Float64} +end + +function (file::FileDef)(t) + # Find index of closest t to tvec + idx = argmin(abs.(file.tvec.- t)) + # return var at index + file.var[idx] +end diff --git a/src/lowOrder2D/LVE.jl b/src/lowOrder2D/LVE.jl new file mode 100644 index 00000000..cb15980b --- /dev/null +++ b/src/lowOrder2D/LVE.jl @@ -0,0 +1,302 @@ +#= +LVE.jl + Written by Matthew White + 5/21/2019 +=# +function LVE(surf::TwoDSurf,curfield::TwoDFlowField,nsteps::Int64 = 500,dtstar::Float64 = 0.015,frameCnt = 20,view = "square",LEV_Stop = 1000,anim = false,animStep = 30,wflag = false,writeInterval = 50) + ## Initialize variables + mat = zeros(8 ,0) # loads output matrix + prevCirc = zeros(surf.ndiv-1) + circChange = zeros(surf.ndiv-1) + vor_loc = zeros(surf.ndiv-1,2) # (x1,z1 ; x2,z2 ...) Inertial frame + coll_loc = zeros(surf.ndiv-1,2) # (x1,z1 ; x2,z2 ...) + global RHS = zeros(surf.ndiv) + global relVel = zeros(surf.ndiv-1,2) + IFR_vor = 0 .* vor_loc # Global frame + IFR_field = TwoDFlowField() + IFR_surf = TwoDSurf(surf.coord_file,surf.pvt,surf.kindef,surf.lespcrit; ndiv = surf.ndiv, camberType = "linear") + global dp = zeros(surf.ndiv-1,1) # change in panel pressure + frameStep = round((nsteps-1) / frameCnt) # Frames to capture animation on + frameStep < 1 ? frameStep = 1 : frameStep + frames = 0 + prevNow = 0 + LEVhist = hcat([],[]) + TEVhist = hcat([],[]) + animMat = hcat([],[],[],[],[],[],[],[]) + + # length of each panel + ds = sqrt.(( surf.x[1:end-1]-surf.x[2:end]).^2 + (surf.cam[1:end-1]-surf.cam[2:end]).^2) + # Surf cam_slope does not give the correct panel locations + cam_slope = asind.((surf.cam[2:end]-surf.cam[1:end-1])./ds) # [deg] + + # Normal vectors for each panel + n = hcat(-sind.(cam_slope),cosd.(cam_slope)) + # Tangential vectors for each panel + tau = hcat(cosd.(cam_slope),sind.(cam_slope)) + + ## Vortex Locations + # Located at 25% of each panel + # Vortex loactions (at 25% of panel length) + vor_loc[:,1] = surf.x[1:end-1] + .25 .* ds .* cosd.(cam_slope) + vor_loc[:,2] = surf.cam[1:end-1] + .25 .* ds .* sind.(cam_slope) + # Collocation points (at 75% panel chordwise) + coll_loc[:,1] = surf.x[1:end-1] + .75 .* ds .* cosd.(cam_slope) + coll_loc[:,2] = surf.cam[1:end-1] + .75 .* ds .* sind.(cam_slope) + + ## Critical gamma (LEV) + x = 1 - 2*ds[1] / surf.c + gamma_crit = surf.lespcrit / 1.13 * ( surf.kinem.u * surf.c * ( acos(x) + sin(acos(x) ))) + gamma_crit = gamma_crit[1] + + ## Refresh vortex values + refresh_vortex(surf,vor_loc) + + ## Time Stepping Loop + t = -dtstar + + for i = 0:nsteps + t += dtstar + + # Update Kinematics + update_kinem(surf, t) + + # Inertial reference frame + IFR_vor[:,1],IFR_vor[:,2] = IFR(surf,vor_loc[:,1],vor_loc[:,2],t) + IFR_surf = refresh_vortex(IFR_surf,IFR_vor) + + ## TEV setup + place_tev2(IFR_surf,IFR_field,dtstar,t) + + ## Influence Coeffcients + x_w, z_w = BFR(surf, IFR_field.tev[end].x, IFR_field.tev[end].z, t) + global a = influence_coeff(surf,curfield,coll_loc,n,dtstar,x_w,z_w) + + ## RHS Vector + RHS[end] = sum(prevCirc) # previous total circulation + # u_w, w_w + if length(curfield.tev) > 0 # Trailing edge vortexs exist + surf.uind[1:end-1], surf.wind[1:end-1] = ind_vel([curfield.tev; curfield.lev],coll_loc[:,1],coll_loc[:,2]) + end + for j = 1:surf.ndiv-1 + alpha = surf.kinem.alpha + # relVel = [U(t) ; W(t)] + global relVel[j,:] = [cos(alpha) -sin(alpha) ; sin(alpha) cos(alpha)] * [surf.kinem.u ; -surf.kinem.hdot] + [-surf.kinem.alphadot*coll_loc[j,2] ; surf.kinem.alphadot * (coll_loc[j,1] - surf.pvt)] + + # RHS = -[U_t + u_w , W_t + w_w] dot n + global RHS[j] = -((relVel[j,1] + surf.uind[j])*n[j,1] + (relVel[j,2] + surf.wind[j])*n[j,2]) + end + + ## Vortex Strenghths a*circ = RHS + global circ = a\RHS + + ## Circulation changes before each panel + for j = 1:surf.ndiv-1 + gamma1 = sum(prevCirc[1:j-1]) + gamma2 = sum(circ[1:j-1]) + circChange[j] = (gamma2-gamma1) ./ dtstar + end + + ## LEV Ejection + # Test for LEV criteria + if abs(circ[1]) > gamma_crit + if surf.levflag[1] == 0 && t < LEV_Stop # if lev is first in batch + gamma_crit_use = abs(circ[1])/circ[1] * gamma_crit # = gamma_crit w/ sign of circ[1] + global t_start = t + #println("LEV start batch:") + else # if lev sheds but is not first in batch + m_slp = -surf.lespcrit / (LEV_Stop - t_start) + c_slp = surf.lespcrit - (m_slp*t_start) + LESP_crit_use = (m_slp*t) + c_slp + x = 1 - 2*ds[1] / surf.c + gamma_crit_use = abs(circ[1])/circ[1] * LESP_crit_use / 1.13 * ( surf.kinem.u * surf.c * ( acos(x) + sin(acos(x) ))) # abs(circ[1])/circ[1] keeps sign of original LEV + end + #println("LEV Ejection: Step $i, time $t") + gamma_crit_use = gamma_crit_use[1] + + # Place LEV + place_lev2(surf,IFR_field,dtstar,t) + + surf.levflag[1] = 1 # Set LEV flag to track constant LEV streams + + # Convert coords to BFR + x_lev, z_lev = BFR(surf, IFR_field.lev[end].x, IFR_field.lev[end].z, t) + global a,a1 = mod_influence_coeff(surf,curfield,coll_loc,n,dtstar,x_w,z_w,x_lev,z_lev) + + # Recalculate RHS + # u_w, w_w + if length(curfield.tev) > 0 # Trailing edge vortexs exist + surf.uind[1:end-1], surf.wind[1:end-1] = ind_vel([curfield.tev; curfield.lev],coll_loc[:,1],coll_loc[:,2]) + end + for j = 1:surf.ndiv-1 + # RHS = -[U_t + u_w , W_t + w_w] dot n + global RHS[j] = -((relVel[j,1] + surf.uind[j])*n[j,1] + (relVel[j,2] + surf.wind[j])*n[j,2]) - gamma_crit_use*a1[j] + end + RHS[end] -= gamma_crit_use # first RHS - gamma_crit_use*a1[end] (which is always 1) + + global circ = a\RHS # Solve new circluations + + IFR_field.lev[end].s = circ[end] #store LEV strength + + circ[2:end] = circ[1:end-1] # shift circs + circ[1] = gamma_crit_use + + # LEV history data [t,s] + stamp = hcat(t,IFR_field.lev[end].s) + LEVhist = vcat(LEVhist,stamp) + else + surf.levflag[1] = 0 + end + + prevCirc = circ[1:end-1] + + for j = 1:surf.ndiv-1 # Set bv circs + surf.bv[j].s = circ[j] + IFR_surf.bv[j].s = circ[j] + end + if length(IFR_field.tev) > 1 # Set TEV circ + IFR_field.tev[end].s = circ[end] + end + + # LEV history data [t,s] + stamp = hcat(t,IFR_field.tev[end].s) + TEVhist = vcat(TEVhist,stamp) + + #= Position mesh + if i % frameStep == 0. && Int(i / frameStep) < frameCnt + + temp = Int(i / frameStep) + 1 + now = Dates.format(Dates.now(), "HH:MM") + if prevNow == now + println("$temp / $frameCnt") + else + println("$temp / $frameCnt, $now") + end + prevNow = now + # Velocities at mesh points + global mesh = meshgrid(surf,IFR_field.tev,IFR_field.lev,.25,t,100,view) + vorts = [IFR_field.tev ; IFR_field.lev ; IFR_surf.bv] + for j = 1:size(vorts,1) + u,w = ind_vel(vorts[j],mesh.x,mesh.z) + global mesh.uMat = mesh.uMat .+ u + global mesh.wMat = mesh.wMat .+ w + global mesh.velMag = sqrt.(mesh.uMat.^2 + mesh.wMat.^2) + global mesh.t = t + end + if frames == 0 + frames = [mesh;] + else + push!(frames,mesh) + end + end + =# + if i > 0 # Force Calculation + ## Pressures (Change in pressure for each panel) + for j = 1:length(dp) + global dp[j] = surf.rho * ( ( (relVel[j,1] + surf.uind[j])*tau[j,1] + (relVel[j,2] + surf.wind[j])*tau[j,2] ) * circ[j]/ds[j] + circChange[j]) + end + ## Loads + # Cn normal + ndem = .5 * surf.rho * surf.kinem.u^2 * surf.c # non-dimensionalization constant + cn = sum( dp[j]*ds[j] for j = 1:length(dp) ) / ndem + # LESP + x = 1 - 2*ds[1] / surf.c + surf.a0[1] = 1.13 * circ[1] / ( surf.kinem.u * surf.c * ( acos(x) + sin(acos(x) ))) + # Cs suction + cs = 2*pi*surf.a0[1]^2 + # Cl lift + alpha = surf.kinem.alpha + cl = cn*cos(alpha) + cs*sin(alpha) + # Cd drag + cd = cn*sin(alpha) - cs*cos(alpha) + # Cm moment + ndem = ndem * surf.c # change ndem to moment ndem + # Cm = Cmz + Cmx + cm = -sum( dp[j]*ds[j]*cosd(cam_slope[j])*(vor_loc[j,1]-surf.pvt) for j = 1:length(dp) ) / ndem + sum( dp[j]*ds[j]*sind(cam_slope[j])*vor_loc[j,2] for j = 1:length(dp) ) / ndem + + mat = hcat(mat,[t, surf.kinem.alpha*180/pi, surf.kinem.h, surf.kinem.u, surf.a0[1], cl, cd, cm]) # Pressures and loads + + end + ## Wake Rollup + IFR_field = wakeroll(IFR_surf,IFR_field,dtstar) + + # Convert IFR_field values to curfield (IFR -> BFR) + push!(curfield.tev,TwoDVort(0,0,IFR_field.tev[end].s,0.02*surf.c,0.,0.)) + if size(IFR_field.lev,1) > 0 && surf.levflag[1] == 1 + push!(curfield.lev,TwoDVort(0,0,IFR_field.lev[end].s,0.02*surf.c,0.,0.)) + end + vor_BFR(surf,t+dtstar,IFR_field,curfield) # Calculate IFR positions for next iteration + + + # Draw animation forces and vorticies alonglside calculation + if anim + if i % animStep == 0. && i > 0 + clf() + # Force Data + animMat = [animMat;mat[:,end]'] + cl = animMat[:,6] + cd = animMat[:,7] + cm = animMat[:,8] + tvec = animMat[:,1] + tmax = nsteps*dtstar + + # Vorticies + animMesh = meshgrid(surf,IFR_field.tev,IFR_field.lev,.25,t,100,view) + + # Calculate velocities + + # Plotting + subplot2grid((3,3),(0,0)) # Top left (Cl) + plot(tvec,cl,color = "black") + xlim(0,tmax) + ylabel("Cl") + + subplot2grid((3,3),(1,0)) # Middle left (Cd) + plot(tvec,cd,color = "black") + xlim(0,tmax) + ylabel("Cd") + + subplot2grid((3,3),(2,0)) # Bottom left (Cm) + plot(tvec,cm,color = "black") + xlim(0,tmax) + ylabel("Cm") + + subplot2grid((3,3),(0,1), colspan=2,rowspan=3,aspect=1) # Vortex plot + + if length(animMesh.tev) > 0 + scatter(animMesh.tev[:,1], animMesh.tev[:,2],c = "red", marker = "o", s = 1) + end + if length(animMesh.lev) > 0 + scatter(animMesh.lev[:,1], animMesh.lev[:,2],c = "blue", marker = "o", s = 1) + end + plot(animMesh.camX,animMesh.camZ,color = "black", linewidth = 3) + xmin = animMesh.x[1,1] + xmax = animMesh.x[1,end] + zmin = animMesh.z[end,1] + zmax = animMesh.z[1,1] + xlim(xmin,xmax) + ylim(zmin,zmax) + axis("off") + + tight_layout() + pause(0.0001) + show(block = false) + print("") # needed becasue for some reason the plot wont update without printing to terminal + end + end + + # Write iteration step variables out if wflag == true + if wflag + if i % writeInterval == 0 && i > 0 + dirname = "$(round(t, sigdigits=6))" + writeStamp(dirname, t, IFR_surf, IFR_field) + end + end + + + end + + mat = mat' + test1 = TEVhist + test2 = LEVhist + return frames,IFR_field,mat,test1,test2 +end diff --git a/src/lowOrder2D/calcs.jl b/src/lowOrder2D/calcs.jl index 92b05eba..420ac1bf 100755 --- a/src/lowOrder2D/calcs.jl +++ b/src/lowOrder2D/calcs.jl @@ -42,7 +42,7 @@ function add_indbound_b(surf::TwoDSurf, surfj::TwoDSurf) surf.uind[:] += uind surf.wind[:] += wind return surf - + end # Function for updating the downwash @@ -106,6 +106,9 @@ function update_kinem(surf::TwoDSurf, t) elseif (typeof(surf.kindef.alpha) == EldRampReturnDef) surf.kinem.alpha = surf.kindef.alpha(t) surf.kinem.alphadot = ForwardDiff.derivative(surf.kindef.alpha,t)*surf.uref/surf.c + elseif (typeof(surf.kindef.alpha) == EldRampReturntstartDef) + surf.kinem.alpha = surf.kindef.alpha(t) + surf.kinem.alphadot = ForwardDiff.derivative(surf.kindef.alpha,t)*surf.uref/surf.c elseif (typeof(surf.kindef.alpha) == ConstDef) surf.kinem.alpha = surf.kindef.alpha(t) surf.kinem.alphadot = 0. @@ -115,6 +118,9 @@ function update_kinem(surf::TwoDSurf, t) elseif (typeof(surf.kindef.alpha) == CosDef) surf.kinem.alpha = surf.kindef.alpha(t) surf.kinem.alphadot = ForwardDiff.derivative(surf.kindef.alpha,t)*surf.uref/surf.c + elseif (typeof(surf.kindef.alpha) == FileDef) + surf.kinem.alpha = surf.kindef.alpha(t) + surf.kinem.alphadot = ForwardDiff.derivative(surf.kindef.alpha,t)*surf.uref/surf.c elseif (typeof(surf.kindef.alpha) == VAWTalphaDef) surf.kinem.alpha = surf.kindef.alpha(t) surf.kinem.alphadot = ForwardDiff.derivative(surf.kindef.alpha,t)*surf.uref/surf.c @@ -137,6 +143,9 @@ function update_kinem(surf::TwoDSurf, t) elseif (typeof(surf.kindef.h) == EldRampReturnDef) surf.kinem.h = surf.kindef.h(t)*surf.c surf.kinem.hdot = ForwardDiff.derivative(surf.kindef.h,t)*surf.uref + elseif (typeof(surf.kindef.h) == EldRampReturntstartDef) + surf.kinem.h = surf.kindef.h(t)*surf.c + surf.kinem.hdot = ForwardDiff.derivative(surf.kindef.h,t)*surf.uref elseif (typeof(surf.kindef.h) == ConstDef) surf.kinem.h = surf.kindef.h(t)*surf.c surf.kinem.hdot = 0. @@ -146,6 +155,9 @@ function update_kinem(surf::TwoDSurf, t) elseif (typeof(surf.kindef.h) == CosDef) surf.kinem.h = surf.kindef.h(t)*surf.c surf.kinem.hdot = ForwardDiff.derivative(surf.kindef.h,t)*surf.uref + elseif (typeof(surf.kindef.h) == FileDef) + surf.kinem.h = surf.kindef.h(t) + surf.kinem.hdot = ForwardDiff.derivative(surf.kindef.h,t)*surf.uref elseif (typeof(surf.kindef.h) == VAWThDef) surf.kinem.h = surf.kindef.h(t)*surf.c surf.kinem.hdot = ForwardDiff.derivative(surf.kindef.h,t)*surf.uref @@ -159,6 +171,9 @@ function update_kinem(surf::TwoDSurf, t) elseif (typeof(surf.kindef.u) == EldRampReturnDef) surf.kinem.u = surf.kindef.u(t)*surf.uref surf.kinem.udot = ForwardDiff.derivative(surf.kindef.u,t)*surf.uref*surf.uref/surf.c + elseif (typeof(surf.kindef.u) == EldRampReturntstartDef) + surf.kinem.u = surf.kindef.u(t)*surf.uref + surf.kinem.udot = ForwardDiff.derivative(surf.kindef.u,t)*surf.uref*surf.uref/surf.c elseif (typeof(surf.kindef.u) == ConstDef) surf.kinem.u = surf.kindef.u(t)*surf.uref surf.kinem.udot = 0. @@ -171,6 +186,9 @@ function update_kinem(surf::TwoDSurf, t) elseif (typeof(surf.kindef.u) == LinearDef) surf.kinem.u = surf.kindef.u(t)*surf.uref surf.kinem.udot = ForwardDiff.derivative(surf.kindef.u,t)*surf.uref*surf.uref/surf.c + elseif (typeof(surf.kindef.u) == FileDef) + surf.kinem.u = surf.kindef.u(t) + surf.kinem.udot = ForwardDiff.derivative(surf.kindef.u,t)*surf.uref*surf.uref/surf.c elseif (typeof(surf.kindef.u) == VAWTuDef) surf.kinem.u = surf.kindef.u(t)*surf.uref surf.kinem.udot = ForwardDiff.derivative(surf.kindef.u,t)*surf.uref*surf.uref/surf.c @@ -408,7 +426,7 @@ function place_lev(surf::TwoDSurf,field::TwoDFlowField,dt) end function place_lev(surf::Vector{TwoDSurf},field::TwoDFlowField,dt,shed_ind::Vector{Int}) - nsurf = length(surf) + nsurf = length(surf) nlev = length(field.lev) for i = 1:nsurf @@ -422,7 +440,7 @@ function place_lev(surf::Vector{TwoDSurf},field::TwoDFlowField,dt,shed_ind::Vect else xloc = surf[i].bnd_x[1]+(1. /3.)*(field.lev[nlev-nsurf+i].x - surf[i].bnd_x[1]) zloc = surf[i].bnd_z[1]+(1. /3.)*(field.lev[nlev-nsurf+i].z - surf[i].bnd_z[1]) - push!(field.lev,TwoDVort(xloc,zloc,0.,0.02*surf[i].c,0.,0.)) + push!(field.lev,TwoDVort(xloc,zloc,0.,0.02*surf[i].c,0.,0.)) end else push!(field.lev, TwoDVort(0., 0., 0., 0., 0., 0.)) @@ -457,6 +475,17 @@ function ind_vel(src::Vector{TwoDVort},t_x,t_z) end return uind, wind end +# Same as above except allows for only one vortex and multiple collocation points +# Used for influence coefficients +function ind_vel(src::TwoDVort,t_x,t_z) + #'s' stands for source and 't' for target + xdist = src.x .- t_x + zdist = src.z .- t_z + distsq = xdist.*xdist + zdist.*zdist + uind = -src.s.*zdist./(2*pi*sqrt.(src.vc*src.vc*src.vc*src.vc .+ distsq.*distsq)) + wind = src.s.*xdist./(2*pi*sqrt.(src.vc*src.vc*src.vc*src.vc .+ distsq.*distsq)) + return uind, wind +end # Function determining the effects of interacting vortices - velocities induced on each other - classical n-body problem function mutual_ind(vorts::Vector{TwoDVort}) @@ -600,14 +629,14 @@ function update_kinem2DOF(surf::TwoDSurf, strpar :: TwoDOFPar, kinem :: KinemPar R1 = 4*strpar.kappa*surf.uref*surf.uref*cl/(pi*surf.c*surf.c) - 2*strpar.w_h*strpar.w_h*(strpar.cubic_h_1*kinem.h + strpar.cubic_h_3*kinem.h^3)/surf.c - strpar.x_alpha*sin(kinem.alpha)*kinem.alphadot*kinem.alphadot R2 = 8*strpar.kappa*surf.uref*surf.uref*cm/(pi*surf.c*surf.c) - strpar.w_alpha*strpar.w_alpha*strpar.r_alpha*strpar.r_alpha*(strpar.cubic_alpha_1*kinem.alpha + strpar.cubic_alpha_3*kinem.alpha^3) - + kinem.hddot_pr = (1/(m11*m22-m21*m12))*(m22*R1-m12*R2) kinem.alphaddot_pr = (1/(m11*m22-m21*m12))*(-m21*R1+m11*R2) #Kinematics are updated according to the 2DOF response kinem.alphadot = kinem.alphadot_pr + (dt/12.)*(23*kinem.alphaddot_pr - 16*kinem.alphaddot_pr2 + 5*kinem.alphaddot_pr3) kinem.hdot = kinem.hdot_pr + (dt/12)*(23*kinem.hddot_pr-16*kinem.hddot_pr2 + 5*kinem.hddot_pr3) - + kinem.alpha = kinem.alpha_pr + (dt/12)*(23*kinem.alphadot_pr-16*kinem.alphadot_pr2 + 5*kinem.alphadot_pr3) kinem.h = kinem.h_pr + (dt/12)*(23*kinem.hdot_pr - 16*kinem.hdot_pr2 + 5*kinem.hdot_pr3) @@ -616,9 +645,192 @@ function update_kinem2DOF(surf::TwoDSurf, strpar :: TwoDOFPar, kinem :: KinemPar surf.kinem.hdot = kinem.hdot surf.kinem.alpha = kinem.alpha surf.kinem.h = kinem.h - + return surf, kinem end +function IFR(surf::TwoDSurf,x,z,t::Float64) + # Given body frame and kinematics, find global positions + alpha = surf.kindef.alpha(t) + X0 = -surf.kindef.u(t)*t + Z0 = surf.kindef.h(t) + pvt = surf.pvt + + X = (x .- pvt).*cos(alpha) + z.*sin(alpha) .+ X0 .+ pvt + Z = -(x .- pvt).*sin(alpha) + z.*cos(alpha) .+ Z0 + + return X,Z +end + +function BFR(surf::TwoDSurf,X,Z,t::Float64) + # Given inertial frame and kinematics, find body positions + alpha = surf.kindef.alpha(t) + X0 = -surf.kindef.u(t)*t + Z0 = surf.kindef.h(t) + pvt = surf.pvt + + x = (X .- X0 .- pvt).*cos(alpha) - (Z .- Z0).*sin(alpha) .+ pvt + z = (X .- X0 .- pvt).*sin(alpha) + (Z .- Z0).*cos(alpha) + + return x,z +end + + +function refresh_vortex(surf::TwoDSurf,vor_loc) + # Updates vortex locations to vor_loc + for i=1:length(surf.bv) + surf.bv[i].x = vor_loc[i,1] + surf.bv[i].z = vor_loc[i,2] + end + return surf +end + +function place_tev2(surf::TwoDSurf,field::TwoDFlowField,dt,t) + TE_x, TE_z = IFR(surf,surf.x[end],surf.cam[end],t) + #= + if size(field.tev,1) == 0 + xloc = surf.bv[end].x + 0.5*surf.kinem.u*dt + zloc = surf.bv[end].z + else + xloc = surf.bv[end].x + (1. /3.)*(field.tev[end].x - surf.bv[end].x) + zloc = surf.bv[end].z + (1. /3.)*(field.tev[end].z - surf.bv[end].z) + end + =# + if size(field.tev,1) == 0 + xloc = TE_x + 0.5*surf.kinem.u*dt + zloc = TE_z + 0.5*surf.kinem.hdot*dt + else + xloc = TE_x + (1. /3.)*(field.tev[end].x - TE_x) + zloc = TE_z + (1. /3.)*(field.tev[end].z - TE_z) + end + push!(field.tev,TwoDVort(xloc,zloc,0.,0.02*surf.c,0.,0.)) + return field +end - +# Places a leading edge vortex +function place_lev2(surf::TwoDSurf,field::TwoDFlowField,dt,t) + LE_x, LE_z = IFR(surf,surf.x[1],surf.cam[1],t) + + le_vel_x = surf.kinem.u - surf.kinem.alphadot*sin(surf.kinem.alpha)*surf.pvt*surf.c + surf.uind[1] + le_vel_z = -surf.kinem.alphadot*cos(surf.kinem.alpha)*surf.pvt*surf.c- surf.kinem.hdot + surf.wind[1] + + if (surf.levflag[1] == 0) + xloc = LE_x + 0.5*le_vel_x*dt + zloc = LE_z + 0.5*le_vel_z*dt + else + xloc = LE_x + (1. /3.)*(field.lev[end].x - LE_x) + zloc = LE_z + (1. /3.)*(field.lev[end].z - LE_z) + end + + push!(field.lev,TwoDVort(xloc,zloc,0.,0.02*surf.c,0.,0.)) + + return field +end + +function vor_BFR(surf::TwoDSurf,t,IFR_field::TwoDFlowField,curfield::TwoDFlowField) + # finds and updates the BFR positions of the vortexes in IFR_field to + # curfield + ntev = size(IFR_field.tev,1) + nlev = size(IFR_field.lev,1) + vorX = [map(q -> q.x,IFR_field.tev); map(q -> q.x,IFR_field.lev)] + vorZ = [map(q -> q.z,IFR_field.tev); map(q -> q.z,IFR_field.lev)] + + # Body frame conversion + vorx, vorz = BFR(surf,vorX,vorZ,t) + + # Update bfr TEV + for i = 1:ntev + curfield.tev[i].x = vorx[i] + curfield.tev[i].z = vorz[i] + end + # Update bfr LEV + for i = 1:nlev + curfield.lev[i].x = vorx[i+ntev] + curfield.lev[i].z = vorz[i+ntev] + end + return curfield +end + +function influence_coeff(surf::TwoDSurf,curfield::TwoDFlowField,coll_loc,n,dt,x_w,z_w) + # With the surface and flowfield, determine the influence matrix "a" + a = zeros(surf.ndiv,surf.ndiv) + a[end,:] .= 1. # for wake portion + vc = surf.bv[1].vc + + for i = 1:surf.ndiv-1, j = 1:surf.ndiv + # i is test location, j is vortex source + t_x = coll_loc[i,1] + t_z = coll_loc[i,2] + + if j < surf.ndiv # Bound vortices (a_ij) + src = surf.bv[j] + xdist = src.x .- t_x + zdist = src.z .- t_z + + distsq = xdist.*xdist + zdist.*zdist # dist_type 1 + u = -zdist./(2*pi*distsq) + w = xdist./(2*pi*distsq) + else # Wake vorticy (a_iw) + xdist = x_w .- t_x + zdist = z_w .- t_z + + distsq = xdist.*xdist + zdist.*zdist # dist type 2 + u = -zdist./(2*pi*sqrt.(vc*vc*vc*vc .+ distsq.*distsq)) + w = xdist./(2*pi*sqrt.(vc*vc*vc*vc .+ distsq.*distsq)) + end + + a[i,j] = u*n[i,1] + w*n[i,2] + end + return a +end + +function mod_influence_coeff(surf::TwoDSurf,curfield::TwoDFlowField,coll_loc,n,dt,x_w,z_w,x_lev,z_lev) + # With the surface and flowfield, determine the influence matrix "a" including the + # modifications needed to calculate the LEV strength + a = zeros(surf.ndiv,surf.ndiv) + a[end,:] .= 1. # for wake portion + vc = surf.bv[1].vc + + for i = 1:surf.ndiv-1, j = 1:surf.ndiv + # i is test location, j is vortex source + t_x = coll_loc[i,1] + t_z = coll_loc[i,2] + + if j < surf.ndiv # Bound vortices (a_ij) + src = surf.bv[j] + xdist = src.x .- t_x + zdist = src.z .- t_z + + distsq = xdist.*xdist + zdist.*zdist # dist_type 1 + u = -zdist./(2*pi*distsq) + w = xdist./(2*pi*distsq) + else # Wake vorticy (a_iw) + xdist = x_w .- t_x + zdist = z_w .- t_z + + distsq = xdist.*xdist + zdist.*zdist # dist type 2 + u = -zdist./(2*pi*sqrt.(vc*vc*vc*vc .+ distsq.*distsq)) + w = xdist./(2*pi*sqrt.(vc*vc*vc*vc .+ distsq.*distsq)) + end + + a[i,j] = u*n[i,1] + w*n[i,2] + end + a1 = a[:,1] + a[:,1:end-1] = a[:,2:end] # shift over a matrix + + for i = 1:surf.ndiv-1 # calc LEV terms + t_x = coll_loc[i,1] + t_z = coll_loc[i,2] + + xdist = x_lev .- t_x + zdist = z_lev .- t_z + + distsq = xdist.*xdist + zdist.*zdist # dist type 2 + u = -zdist./(2*pi*sqrt.(vc*vc*vc*vc .+ distsq.*distsq)) + w = xdist./(2*pi*sqrt.(vc*vc*vc*vc .+ distsq.*distsq)) + + a[i,end] = u*n[i,1] + w*n[i,2] + end + + return a,a1 +end diff --git a/src/lowOrder2D/postprocess.jl b/src/lowOrder2D/postprocess.jl index 24236e1b..ccf08dd6 100755 --- a/src/lowOrder2D/postprocess.jl +++ b/src/lowOrder2D/postprocess.jl @@ -28,10 +28,18 @@ function calc_forces(surf::TwoDSurf, vels::Vector{Float64}) #Pitching moment is clockwise or nose up positive cm = cn*surf.pvt - 2*pi*(((surf.kinem.u + vels[1])*cos(surf.kinem.alpha)/surf.uref + (surf.kinem.hdot - vels[2])*sin(surf.kinem.alpha)/surf.uref)*(surf.a0[1]/4. + surf.aterm[1]/4. - surf.aterm[2]/8.) + (surf.c/surf.uref)*(7. *surf.a0dot[1]/16. + 3. *surf.adot[1]/16. + surf.adot[2]/16. - surf.adot[3]/64.)) - nonl_m - return cl, cd, cm + return cl, cd, cm, cn end function writeStamp(dirname::String, t::Float64, surf::TwoDSurf, curfield::TwoDFlowField) + + try + cd("Step Files") + catch + mkdir("Step Files") + cd("Step Files") + end + dirvec = readdir() if dirname in dirvec rm(dirname, recursive=true) @@ -99,7 +107,8 @@ function writeStamp(dirname::String, t::Float64, surf::TwoDSurf, curfield::TwoDF mat[:,6] = ql[:] DelimitedFiles.writedlm(f, mat) close(f) - + + cd("..") cd("..") end @@ -180,16 +189,16 @@ function calc_delcp(surf::TwoDSurf, vels::Vector{Float64}) gamint = zeros(surf.ndiv) p_com = zeros(surf.ndiv) gammod = zeros(surf.ndiv) - + for i = 1:surf.ndiv udash = surf.kinem.u*cos(surf.kinem.alpha) + surf.kinem.hdot*sin(surf.kinem.alpha) + surf.uind[i]*cos(surf.kinem.alpha) - surf.wind[i]*sin(surf.kinem.alpha) - + p_in[i] = 8*surf.uref*sqrt(surf.c)*surf.a0[1]*sqrt(surf.x[i])*udash/(surf.rho*surf.c + 2*surf.x[i]) + 4*surf.uref*sqrt(surf.c)*surf.a0dot[1]*sqrt(surf.x[i]) end for i = 2:surf.ndiv udash = surf.kinem.u*cos(surf.kinem.alpha) + surf.kinem.hdot*sin(surf.kinem.alpha) + surf.uind[i]*cos(surf.kinem.alpha) - surf.wind[i]*sin(surf.kinem.alpha) gam[i] = surf.a0[1]*cot(surf.theta[i]/2) - gamint[i] = surf.a0dot[1]*(surf.theta[i] + sin(surf.theta[i])) + gamint[i] = surf.a0dot[1]*(surf.theta[i] + sin(surf.theta[i])) for n = 1:surf.naterm gam[i] += surf.aterm[n]*sin(n*surf.theta[i]) gammod[i] += surf.aterm[n]*sin(n*surf.theta[i]) @@ -202,11 +211,11 @@ function calc_delcp(surf::TwoDSurf, vels::Vector{Float64}) gam[i] = gam[i]*surf.uref gammod[i] = gammod[i]*surf.uref gamint[i] = gamint[i]*surf.uref*surf.c - - p_out[i] = 2*udash*gam[i] + gamint[i] + + p_out[i] = 2*udash*gam[i] + gamint[i] p_com[i] = 2*udash*surf.uref*surf.a0[1]*(cos(surf.theta[i]/2)-1)/sin(surf.theta[i]/2) + 2*udash*gammod[i] + gamint[i] + 8*surf.uref*udash*surf.c*surf.a0[1]*sin(surf.theta[i]/2)/(surf.rho*surf.c + 2*surf.c*(sin(surf.theta[i]/2))^2) end - + return p_com, p_in, p_out end @@ -218,15 +227,15 @@ function calc_edgeVel(surf::TwoDSurf, vels::Vector{Float64}) q_com_u[1] = surf.uref*surf.a0[1]/sqrt(0.5*surf.rho) q_com_l[1] = surf.uref*surf.a0[1]/sqrt(0.5*surf.rho) - + for i = 2:surf.ndiv udash = (surf.kinem.u + vels[1])*cos(surf.kinem.alpha) + (surf.kinem.hdot - vels[2])*sin(surf.kinem.alpha) + surf.uind[i]*cos(surf.kinem.alpha) - surf.wind[i]*sin(surf.kinem.alpha) - + for n = 1:surf.naterm gammod[i] += surf.aterm[n]*sin(n*surf.theta[i]) end gammod[i] = gammod[i]*surf.uref - + q_com_u[i] = surf.uref*surf.a0[1]*(cos(surf.theta[i]/2) - 1)/sin(surf.theta[i]/2) + surf.uref*gammod[i] + (sqrt(surf.c)*surf.uref*surf.a0[1] + sqrt(surf.x[i])*udash)/sqrt(surf.x[i] + surf.rho*surf.c/2) q_com_l[i] = -surf.uref*surf.a0[1]*(cos(surf.theta[i]/2) - 1)/sin(surf.theta[i]/2) - surf.uref*gammod[i] + (-sqrt(surf.c)*surf.uref*surf.a0[1] + sqrt(surf.x[i])*udash)/sqrt(surf.x[i] + surf.rho*surf.c/2) end @@ -237,9 +246,9 @@ end function getEndCycle(mat, k) T = pi/k - + end_cycle = mat[end,1] - + #Find number of cycles ncyc = 0 for i = 1:1000 @@ -249,14 +258,14 @@ function getEndCycle(mat, k) end ncyc = ncyc + 1 end - + start_t = real(ncyc-1)*T end_t = real(ncyc)*T start_ind = argmin(abs.(mat[:,1] .- start_t)) end_ind = argmin(abs.(mat[:,1] .- end_t)) nlast = end_ind - start_ind + 1 - + newmat = zeros(nlast, 8) newmat[:,1] = (mat[start_ind:end_ind,1] .- start_t)/T for i = 2:7 diff --git a/src/lowOrder2D/solvers.jl b/src/lowOrder2D/solvers.jl index ee951f65..e61839c0 100755 --- a/src/lowOrder2D/solvers.jl +++ b/src/lowOrder2D/solvers.jl @@ -57,7 +57,7 @@ function lautat(surf::TwoDSurf, curfield::TwoDFlowField, nsteps::Int64 = 500, dt mat = mat' dt = dtstar*surf.c/surf.uref - + # if writeflag is on, determine the timesteps to write at if writeflag == 1 writeArray = Int64[] @@ -113,7 +113,7 @@ function lautat(surf::TwoDSurf, curfield::TwoDFlowField, nsteps::Int64 = 500, dt if wakerollup == 1 wakeroll(surf, curfield, dt) end - + # Calculate force and moment coefficients cl, cd, cm = calc_forces(surf, [curfield.u[1], curfield.w[1]]) @@ -159,7 +159,7 @@ function ldvm(surf::TwoDSurf, curfield::TwoDFlowField, nsteps::Int64 = 500, dtst mat = mat' dt = dtstar*surf.c/surf.uref - + # if writeflag is on, determine the timesteps to write at if writeflag == 1 writeArray = Int64[] @@ -190,18 +190,18 @@ function ldvm(surf::TwoDSurf, curfield::TwoDFlowField, nsteps::Int64 = 500, dtst #Add a TEV with dummy strength place_tev(surf,curfield,dt) - + #Update incduced velocities on airfoil update_indbound(surf, curfield) - + #Solve for TEV strength to satisfy Kelvin condition kelv = KelvinCondition(surf,curfield) soln = nlsolve(not_in_place(kelv), [-0.01]) curfield.tev[length(curfield.tev)].s = soln.zero[1] - + #Update adot update_a2a3adot(surf,dt) - + lesp = surf.a0[1] #Check for LESP condition @@ -215,7 +215,7 @@ function ldvm(surf::TwoDSurf, curfield::TwoDFlowField, nsteps::Int64 = 500, dtst kelvkutta = KelvinKutta(surf,curfield) soln = nlsolve(not_in_place(kelvkutta), [-0.01; 0.01]) (curfield.tev[length(curfield.tev)].s, curfield.lev[length(curfield.lev)].s) = soln.zero[1], soln.zero[2] - + #set flag for levshedding=on surf.levflag[1] = 1 else @@ -241,7 +241,7 @@ function ldvm(surf::TwoDSurf, curfield::TwoDFlowField, nsteps::Int64 = 500, dtst wakeroll(surf, curfield, dt) # Calculate force and moment coefficients - cl, cd, cm = calc_forces(surf, [curfield.u[1], curfield.w[1]]) + cl, cd, cm, cn = calc_forces(surf, [curfield.u[1], curfield.w[1]]) # write flow details if required if writeflag == 1 @@ -270,7 +270,7 @@ end function ldvm(surf::Vector{TwoDSurf}, curfield::TwoDFlowField, nsteps::Int64 = 500, dtstar::Float64 = 0.015, startflag = 0, writeflag = 0, writeInterval = 1000., delvort = delNone(); maxwrite = 50, nround=6) nsurf = length(surf) - + # If a restart directory is provided, read in the simulation data if startflag == 0 mat = zeros(0, 7*nsurf+1) @@ -287,7 +287,7 @@ function ldvm(surf::Vector{TwoDSurf}, curfield::TwoDFlowField, nsteps::Int64 = 5 mat = mat' dt = dtstar*minimum(map(q->q.c/q.uref, surf)) - + # if writeflag is on, determine the timesteps to write at if writeflag == 1 writeArray = Int64[] @@ -301,12 +301,12 @@ function ldvm(surf::Vector{TwoDSurf}, curfield::TwoDFlowField, nsteps::Int64 = 5 end end end - + lesp = zeros(nsurf) cl = zeros(nsurf) cd = zeros(nsurf) cm = zeros(nsurf) - + # time stepping for istep = 1:nsteps #Udpate current time @@ -318,11 +318,11 @@ function ldvm(surf::Vector{TwoDSurf}, curfield::TwoDFlowField, nsteps::Int64 = 5 for is = 1:nsurf #Update kinematic parameters update_kinem(surf[is], t) - + #Update bound vortex positions update_boundpos(surf[is], dt) end - + #Add TEVs with dummy strength place_tev(surf,curfield,dt) @@ -337,26 +337,26 @@ function ldvm(surf::Vector{TwoDSurf}, curfield::TwoDFlowField, nsteps::Int64 = 5 end end - + #Solve for TEV strength to satisfy Kelvin condition kelv = KelvinConditionMult(surf,curfield) soln = nlsolve(not_in_place(kelv), ones(nsurf)*-0.01) for i = 1:nsurf curfield.tev[end-nsurf+i].s = soln.zero[i] - + #Update adot update_a2a3adot(surf[i],dt) - + lesp[i] = surf[i].a0[1] end - + shed_ind = Int[] - + #Check for LESP condition for i = 1:nsurf - if (abs(lesp[i])>surf[i].lespcrit[1]) - + if (abs(lesp[i])>surf[i].lespcrit[1]) + push!(shed_ind, i) end end @@ -366,11 +366,11 @@ function ldvm(surf::Vector{TwoDSurf}, curfield::TwoDFlowField, nsteps::Int64 = 5 #2D iteration if LESP_crit is exceeded #Add LEVs with dummy strength place_lev(surf,curfield,dt,shed_ind) - + #Solve for TEV and LEV strengths to satisfy Kelvin condition and Kutta condition at leading edge kelvkutta = KelvinKuttaMult(surf,curfield,shed_ind) soln = nlsolve(not_in_place(kelvkutta), -0.01*ones(nshed+nsurf)) - + for i = 1:nsurf if i in shed_ind curfield.tev[end-nsurf+i].s = soln.zero[i] @@ -378,15 +378,15 @@ function ldvm(surf::Vector{TwoDSurf}, curfield::TwoDFlowField, nsteps::Int64 = 5 surf[i].levflag[1] = 1 else curfield.tev[end-nsurf+i].s = soln.zero[i] - surf[i].levflag[1] = 0 + surf[i].levflag[1] = 0 end end else for i = 1:nsurf surf[i].levflag[1] = 0 - end + end end - + for i = 1:nsurf #Set previous values of aterm to be used for derivatives in next time step surf[i].a0prev[1] = surf[i].a0[1] @@ -394,20 +394,20 @@ function ldvm(surf::Vector{TwoDSurf}, curfield::TwoDFlowField, nsteps::Int64 = 5 surf[i].aprev[ia] = surf[i].aterm[ia] end end - + # Delete or merge vortices if required mean_bndx = sum(map(q->q.bnd_x[Int(round(q.ndiv/2))], surf))/nsurf mean_bndz = sum(map(q->q.bnd_z[Int(round(q.ndiv/2))], surf))/nsurf controlVortCount(delvort, mean_bndx, mean_bndz, curfield) - + # free wake rollup wakeroll(surf, curfield, dt) - + # Calculate force and moment coefficients for i = 1:nsurf cl[i], cd[i], cm[i] = calc_forces(surf[i], [curfield.u[1], curfield.w[1]]) end - + # write flow details if required if writeflag == 1 if istep in writeArray @@ -415,7 +415,7 @@ function ldvm(surf::Vector{TwoDSurf}, curfield::TwoDFlowField, nsteps::Int64 = 5 writeStamp(dirname, t, surf, curfield) end end - + # for writing in resultsSummary matvect = [t;] for is = 1:nsurf @@ -423,7 +423,7 @@ function ldvm(surf::Vector{TwoDSurf}, curfield::TwoDFlowField, nsteps::Int64 = 5 surf[is].kinem.u, surf[is].a0[1], cl[is], cd[is], cm[is]]] end mat = hcat(mat,matvect) - + end mat = mat' @@ -456,7 +456,7 @@ function ldvmLin(surf::TwoDSurf, curfield::TwoDFlowField, nsteps::Int64 = 500, d mat = mat' dt = dtstar*surf.c/surf.uref - + # if writeflag is on, determine the timesteps to write at if writeflag == 1 writeArray = Int64[] @@ -665,7 +665,7 @@ function ldvm2DOF(surf::TwoDSurf, curfield::TwoDFlowField, strpar::TwoDOFPar, ki mat = mat' dt = dtstar*surf.c/surf.uref - + # if writeflag is on, determine the timesteps to write at if writeflag == 1 writeArray = Int64[] @@ -681,19 +681,19 @@ function ldvm2DOF(surf::TwoDSurf, curfield::TwoDFlowField, strpar::TwoDOFPar, ki end - + cl = 0. cm = 0. - + #Intialise flowfield - + for istep = 1:nsteps #Udpate current time t = t + dt #Update external flowfield update_externalvel(curfield, t) - + #Update kinematic parameters (based on 2DOF response) update_kinem2DOF(surf, strpar, kinem, dt, cl, cm) @@ -705,7 +705,7 @@ function ldvm2DOF(surf::TwoDSurf, curfield::TwoDFlowField, strpar::TwoDOFPar, ki #Update incduced velocities on airfoil update_indbound(surf, curfield) - + #Solve for TEV strength to satisfy Kelvin condition kelv = KelvinCondition(surf,curfield) soln = nlsolve(not_in_place(kelv), [-0.01]) @@ -714,7 +714,7 @@ function ldvm2DOF(surf::TwoDSurf, curfield::TwoDFlowField, strpar::TwoDOFPar, ki #Update adot #update_a2a3adot(surf,dt) update_atermdot(surf, dt) - + lesp = surf.a0[1] #Check for LESP condition @@ -729,7 +729,7 @@ function ldvm2DOF(surf::TwoDSurf, curfield::TwoDFlowField, strpar::TwoDOFPar, ki kelvkutta = KelvinKutta(surf,curfield) soln = nlsolve(not_in_place(kelvkutta), [-0.01; 0.01]) (curfield.tev[length(curfield.tev)].s, curfield.lev[length(curfield.lev)].s) = soln.zero[1], soln.zero[2] - + #set flag for levshedding=on surf.levflag[1] = 1 else @@ -738,7 +738,7 @@ function ldvm2DOF(surf::TwoDSurf, curfield::TwoDFlowField, strpar::TwoDOFPar, ki #Update rest of Fourier terms update_a2toan(surf) - + #Set previous values of aterm to be used for derivatives in next time step surf.a0prev[1] = surf.a0[1] for ia = 1:surf.naterm @@ -764,14 +764,14 @@ function ldvm2DOF(surf::TwoDSurf, curfield::TwoDFlowField, strpar::TwoDOFPar, ki writeStamp(dirname, t, surf, curfield) end end - + # for writing in resultsSummary mat = hcat(mat,[t, surf.kinem.alpha, surf.kinem.h, surf.kinem.u, surf.a0[1], cl, cd, cm]) - + end - + mat = mat' - + f = open("resultsSummary", "w") Serialization.serialize(f, ["#time \t", "alpha (deg) \t", "h/c \t", "u/uref \t", "A0 \t", "Cl \t", "Cd \t", "Cm \n"]) DelimitedFiles.writedlm(f, mat) @@ -779,4 +779,3 @@ function ldvm2DOF(surf::TwoDSurf, curfield::TwoDFlowField, strpar::TwoDOFPar, ki mat, surf, curfield end - diff --git a/src/lowOrder2D/typedefs.jl b/src/lowOrder2D/typedefs.jl index fa16359a..19fe0fb7 100755 --- a/src/lowOrder2D/typedefs.jl +++ b/src/lowOrder2D/typedefs.jl @@ -64,8 +64,8 @@ struct TwoDSurf levflag :: Vector{Int8} initpos :: Vector{Float64} rho :: Float64 - - function TwoDSurf(coord_file, pvt, kindef,lespcrit=zeros(1); c=1., uref=1., ndiv=70, naterm=35, initpos = [0.; 0.], rho = 0.04) + + function TwoDSurf(coord_file, pvt, kindef,lespcrit=zeros(1); c=1., uref=1., ndiv=70, naterm=35, initpos = [0.; 0.], rho = 0.04, camberType = "radial") theta = zeros(ndiv) x = zeros(ndiv) cam = zeros(ndiv) @@ -74,11 +74,19 @@ struct TwoDSurf bnd_z = zeros(ndiv) kinem = KinemPar(0, 0, 0, 0, 0, 0) - dtheta = pi/(ndiv-1) - for ib = 1:ndiv - theta[ib] = (ib-1.)*dtheta - x[ib] = c/2. *(1-cos(theta[ib])) + if camberType == "radial" + dtheta = pi/(ndiv-1) + for ib = 1:ndiv + theta[ib] = (ib-1.)*dtheta + x[ib] = c/2. *(1-cos(theta[ib])) + end + elseif camberType == "linear" + dx = c / (ndiv-1) + for ib = 2:ndiv + x[ib] = x[ib-1] + dx + end end + if (coord_file != "FlatPlate") cam, cam_slope = camber_calc(x, coord_file) end @@ -102,6 +110,9 @@ struct TwoDSurf elseif (typeof(kindef.alpha) == CosDef) kinem.alpha = kindef.alpha(0.) kinem.alphadot = ForwardDiff.derivative(kindef.alpha,0.)*uref/c + elseif (typeof(kindef.alpha) == FileDef) + kinem.alpha = kindef.alpha(0.) + kinem.alphadot = ForwardDiff.derivative(kindef.alpha,0.)*uref/c end if (typeof(kindef.h) == EldUpDef) @@ -122,6 +133,9 @@ struct TwoDSurf elseif (typeof(kindef.h) == CosDef) kinem.h = kindef.h(0.)*c kinem.hdot = ForwardDiff.derivative(kindef.h,0.)*uref + elseif (typeof(kindef.h) == FileDef) + kinem.h = kindef.h(0.) + kinem.hdot = ForwardDiff.derivative(kindef.h,0.)*uref end if (typeof(kindef.u) == EldUpDef) @@ -134,6 +148,9 @@ struct TwoDSurf elseif (typeof(kindef.u) == ConstDef) kinem.u = kindef.u(0.)*uref kinem.udot = 0. + elseif (typeof(kindef.alpha) == FileDef) + kinem.u = kindef.u(0.)*uref + kinem.udot = ForwardDiff.derivative(kindef.u,0.)*uref*uref/c end for i = 1:ndiv @@ -214,7 +231,7 @@ function (kelv::KelvinCondition)(tev_iter::Array{Float64}) ntev = length(kelv.field.tev) uprev, wprev = ind_vel([kelv.field.tev[ntev]], kelv.surf.bnd_x, kelv.surf.bnd_z) - + #Update the TEV strength kelv.field.tev[ntev].s = tev_iter[1] @@ -222,7 +239,7 @@ function (kelv::KelvinCondition)(tev_iter::Array{Float64}) kelv.surf.uind[:] = kelv.surf.uind[:] .- uprev .+ unow kelv.surf.wind[:] = kelv.surf.wind[:] .- wprev .+ wnow - + #Calculate downwash update_downwash(kelv.surf, [kelv.field.u[1],kelv.field.w[1]]) @@ -246,10 +263,10 @@ function (kelv::KelvinConditionMult)(tev_iter::Array{Float64}) nsurf = length(kelv.surf) val = zeros(nsurf) - + nlev = length(kelv.field.lev) ntev = length(kelv.field.tev) - + tev_list = kelv.field.tev[ntev-nsurf+1:ntev] for i = 1:nsurf @@ -260,7 +277,7 @@ function (kelv::KelvinConditionMult)(tev_iter::Array{Float64}) end end uprev, wprev = ind_vel([tev_list; bv_list], kelv.surf[i].bnd_x, kelv.surf[i].bnd_z) - + #Update the TEV strength kelv.field.tev[ntev-nsurf+i].s = tev_iter[i] @@ -268,18 +285,18 @@ function (kelv::KelvinConditionMult)(tev_iter::Array{Float64}) kelv.surf[i].uind[:] = kelv.surf[i].uind[:] .- uprev .+ unow kelv.surf[i].wind[:] = kelv.surf[i].wind[:] .- wprev .+ wnow - + #Calculate downwash update_downwash(kelv.surf[i], [kelv.field.u[1],kelv.field.w[1]]) end - + for i = 1:nsurf #Update Fourier coefficients and bv strength update_a0anda1(kelv.surf[i]) update_a2toan(kelv.surf[i]) update_bv(kelv.surf[i]) - - + + val[i] = kelv.surf[i].uref*kelv.surf[i].c*pi*(kelv.surf[i].a0[1] + kelv.surf[i].aterm[1]/2.) - kelv.surf[i].uref*kelv.surf[i].c*pi*(kelv.surf[i].a0prev[1] + kelv.surf[i].aprev[1]/2.) + kelv.field.tev[ntev-nsurf+i].s @@ -295,12 +312,12 @@ end function (kelv::KelvinKutta)(v_iter::Array{Float64}) val = zeros(2) - + nlev = length(kelv.field.lev) ntev = length(kelv.field.tev) uprev, wprev = ind_vel([kelv.field.tev[ntev]; kelv.field.lev[nlev]], kelv.surf.bnd_x, kelv.surf.bnd_z) - + #Update the TEV and LEV strengths kelv.field.tev[ntev].s = v_iter[1] kelv.field.lev[nlev].s = v_iter[2] @@ -309,17 +326,17 @@ function (kelv::KelvinKutta)(v_iter::Array{Float64}) kelv.surf.uind[:] = kelv.surf.uind[:] .- uprev .+ unow kelv.surf.wind[:] = kelv.surf.wind[:] .- wprev .+ wnow - + #Calculate downwash update_downwash(kelv.surf ,[kelv.field.u[1],kelv.field.w[1]]) #Calculate first two fourier coefficients update_a0anda1(kelv.surf) - val[1] = kelv.surf.uref*kelv.surf.c*pi*(kelv.surf.a0[1] + kelv.surf.aterm[1]/2.) - + val[1] = kelv.surf.uref*kelv.surf.c*pi*(kelv.surf.a0[1] + kelv.surf.aterm[1]/2.) - kelv.surf.uref*kelv.surf.c*pi*(kelv.surf.a0prev[1] + kelv.surf.aprev[1]/2.) + kelv.field.tev[ntev].s + kelv.field.lev[nlev].s - + if (kelv.surf.a0[1] > 0) lesp_cond = kelv.surf.lespcrit[1] else @@ -343,7 +360,7 @@ function (kelv::KelvinKuttaMult)(v_iter::Array{Float64}) nlev = length(kelv.field.lev) ntev = length(kelv.field.tev) - + tev_list = kelv.field.tev[ntev-nsurf+1:ntev] lev_list = kelv.field.lev[nlev-nsurf.+kelv.shed_ind] @@ -356,20 +373,20 @@ function (kelv::KelvinKuttaMult)(v_iter::Array{Float64}) end end uprev, wprev = ind_vel([tev_list; lev_list; bv_list], kelv.surf[i].bnd_x, kelv.surf[i].bnd_z) - + #Update the shed vortex strengths kelv.field.tev[ntev-nsurf+i].s = v_iter[i] - + if i in kelv.shed_ind levcount += 1 kelv.field.lev[nlev-nsurf+i].s = v_iter[nsurf+levcount] - end + end unow, wnow = ind_vel([tev_list; lev_list; bv_list], kelv.surf[i].bnd_x, kelv.surf[i].bnd_z) kelv.surf[i].uind[:] = kelv.surf[i].uind[:] .- uprev .+ unow kelv.surf[i].wind[:] = kelv.surf[i].wind[:] .- wprev .+ wnow - + #Calculate downwash update_downwash(kelv.surf[i], [kelv.field.u[1],kelv.field.w[1]]) end @@ -391,14 +408,105 @@ function (kelv::KelvinKuttaMult)(v_iter::Array{Float64}) else lesp_cond = -kelv.surf[i].lespcrit[1] end - + val[levcount+nsurf] = kelv.surf[i].a0[1] - lesp_cond else val[i] = kelv.surf[i].uref*kelv.surf[i].c*pi*(kelv.surf[i].a0[1] + kelv.surf[i].aterm[1]/2.) - kelv.surf[i].uref*kelv.surf[i].c*pi*(kelv.surf[i].a0prev[1] + kelv.surf[i].aprev[1]/2.) + kelv.field.tev[ntev-nsurf+i].s end - end + end return val end +mutable struct meshgrid + t :: Float64 + alpha :: Float64 + tev :: Array{Float64} + lev :: Array{Float64} + camX :: Vector{Float64} + camZ :: Vector{Float64} + x :: Array{Float64} + z :: Array{Float64} + uMat :: Array{Float64} + wMat :: Array{Float64} + velMag :: Array{Float64} + + function meshgrid(surf::TwoDSurf,tevs::Vector{TwoDVort},levs::Vector{TwoDVort},offset,t,width::Int64 = 100,view::String = "square") + farBnd = surf.x[end] + surf.c*offset + nearBnd = surf.x[1] - surf.c*offset + zBnd = ( farBnd - nearBnd ) / 2 + if view == "wake" + farBnd = 2*farBnd + elseif view == "largewake" + farBnd = 3*farBnd + zBnd = 2*zBnd + elseif view == "longwake" + farBnd = 5*farBnd + zBnd = 2*zBnd + elseif view == "UI Window" + farBnd = 2*farBnd + zBnd = 3/4*farBnd + end + + # Global frame translation + X0 = -surf.kindef.u(t)*t + Z0 = surf.kindef.h(t) + + lowX = nearBnd + X0 - surf.pvt*surf.c + uppX = farBnd + X0 - surf.pvt*surf.c + lowZ = -zBnd + Z0 + uppZ = zBnd + Z0 + + # Finding global camber line postion + camX, camZ = IFR(surf,surf.x,surf.cam,t) + + # Creating meshgrid + step = (farBnd-nearBnd)/(width-1) + range = lowX:step:uppX + height = zBnd*2/step + 1 + x = [ j for i = 1:height, j = range] + + range = uppZ:-step:lowZ + z = [ i for i = range , j = 1:size(x,2)] + + uMat = 0 .* x .+ surf.kinem.u + wMat = 0 .* x .+ surf.kinem.hdot + velMag = 0 .* x + + t = 0 + alpha = surf.kinem.alpha + circ = zeros(surf.ndiv-1) + + tevX = [] + tevZ = [] + for i = 1:length(tevs) + vx = tevs[i].x + vz = tevs[i].z + if vx >= lowX && vx <= uppX + if vz >= lowZ && vz <= uppZ + tevX = [tevX;vx] + tevZ = [tevZ;vz] + end + end + end + tev = hcat(tevX,tevZ) + levX = [] + levZ = [] + if length(levs) > 0 + for i = 1:length(levs) + vx = levs[i].x + vz = levs[i].z + if vx >= lowX && vx <= uppX + if vz >= lowZ && vz <= uppZ + levX = [levX;vx] + levZ = [levZ;vz] + end + end + end + end + lev = hcat(levX,levZ) + + new(t,alpha,tev,lev,camX, camZ, x, z, uMat, wMat, velMag) + end +end diff --git a/src/plots/plots2D.jl b/src/plots/plots2D.jl index 53f2bd0a..f25b6b01 100755 --- a/src/plots/plots2D.jl +++ b/src/plots/plots2D.jl @@ -10,7 +10,7 @@ function viewVort2D(tev::Vector{TwoDVort}, lev::Vector{TwoDVort}, bv::Vector{Two tev = hcat(map(q->q.s, tev), map(q->q.x, tev), map(q->q.z, tev)) lev = hcat(map(q->q.s, lev), map(q->q.x, lev), map(q->q.z, lev)) bv = hcat(map(q->q.s, bv), map(q->q.x, bv), map(q->q.z, bv)) - + scatter(tev[:,2], tev[:,3], s=5, c=tev[:,1], edgecolors="none") sc = scatter(lev[:,2], lev[:,3], s=5, c=lev[:,1], edgecolors="none") sc2 = scatter(bv[:,2], bv[:,3], s=5, c=bv[:,1], edgecolors="none") @@ -28,10 +28,21 @@ function viewVortConnect2D(tev::Matrix{Float64}, lev::Matrix{Float64}, bv::Matri plot(bv[:,2], bv[:,3], color = "black", linewidth=1.0) end -function makeVortPlots2D() +function makeVortPlots2D(tt = "all") +# tt is a target time to make the plot for. It will choose the closest time available + try + cd("Step Files") + catch + println(" Error: No 'Step Files' directory. Make sure 'wflag' is enabled before running simulation.") + return + end dirvec = readdir() dirresults = map(x->(v = tryparse(Float64,x); typeof(v) == Nothing ? 0.0 : v),dirvec) + if tt != "all" + figure() + ttidx = findmin(abs.(dirresults.-tt))[2] + end #Determine axis limits dirmax = maximum(dirresults) cd("$(dirmax)") @@ -40,7 +51,7 @@ function makeVortPlots2D() if isfile("boundv-1") == true #only 1 surface multsurfflag = 1 end - + if multsurfflag == 0 #single surface tev, _ = DelimitedFiles.readdlm("tev", '\t', Float64, header=true) @@ -56,13 +67,21 @@ function makeVortPlots2D() xmax = maximum([tev[:,2];lev[:,2];]) zmax = maximum([bv[:,2];tev[:,3];lev[:,3];bv[:,3];]) - cd("..") + cd("../..") - if "vortPlots" in dirvec - rm("vortPlots", recursive=true) + if tt == "all" + dirvec = readdir() + if "vortPlots" in dirvec + rm("vortPlots", recursive=true) + end + mkdir("vortPlots") end - mkdir("vortPlots") + cd("Step Files") + for i=1:length(dirresults) + if tt != "all" && ttidx != i + continue + end if dirresults[i] != 0 dirstr="$(dirresults[i])" cd(dirstr) @@ -76,11 +95,20 @@ function makeVortPlots2D() viewVort2D(tev, lev, bv) axis([xmin-1, xmax+1, zmin-1, zmax+1]) - savefig("../vortPlots/$(dirresults[i]).png") - close() + + if tt == "all" + savefig("../../vortPlots/$(dirresults[i]).png") + close() + end cd("..") end end + if tt == "all" + println(" All plots saved under the 'vortPlots' directory.") + else + pause(0.0001) + show(block = false) + end else #multiple surfaces @@ -138,16 +166,19 @@ function makeVortPlots2D() end end end + cd("..") end -function makeForcePlots2D() +function makeForcePlots2D(mat = "") dirvec = readdir() if "forcePlots" in dirvec rm("forcePlots", recursive=true) end mkdir("forcePlots") - mat, _ = DelimitedFiles.readdlm("resultsSummary", '\t', Float64, header=true) + if mat == "" + mat, _ = DelimitedFiles.readdlm("resultsSummary", '\t', Float64, header=true) + end if length(mat[1,:]) == 8 #only 1 surface t = mat[:,1] @@ -293,9 +324,9 @@ function makeForcePlots2D() end function checkConverge(k::Float64) - + mat, _ = DelimitedFiles.readdlm("resultsSummary", '\t', Float64, header=true) - + T = pi/k end_cycle = mat[end,1] @@ -310,7 +341,7 @@ function checkConverge(k::Float64) ncyc = ncyc + 1 end - #Lift convergence + #Lift convergence figure for i = 1:ncyc start_t = real(i-1)*T @@ -323,14 +354,14 @@ function checkConverge(k::Float64) xmax = 1. zmin = minimum(mat[:,6]) zmax = maximum(mat[:,6]) - axis([xmin, xmax, zmin, zmax]) + axis([xmin, xmax, zmin, zmax]) xlabel(L"$t^*$") ylabel(L"$C_l$") savefig("forcePlots/cl-convergence.png") close() - #Drag convergence + #Drag convergence figure for i = 1:ncyc start_t = real(i-1)*T @@ -343,14 +374,14 @@ function checkConverge(k::Float64) xmax = 1. zmin = minimum(mat[:,7]) zmax = maximum(mat[:,7]) - axis([xmin, xmax, zmin, zmax]) + axis([xmin, xmax, zmin, zmax]) xlabel(L"$t^*$") ylabel(L"$C_d$") savefig("forcePlots/cd-convergence.png") close() - #Pitching moment convergence + #Pitching moment convergence figure for i = 1:ncyc start_t = real(i-1)*T @@ -363,7 +394,7 @@ function checkConverge(k::Float64) xmax = 1. zmin = minimum(mat[:,8]) zmax = maximum(mat[:,8]) - axis([xmin, xmax, zmin, zmax]) + axis([xmin, xmax, zmin, zmax]) xlabel(L"$t^*$") ylabel(L"$C_m$") @@ -371,9 +402,12 @@ function checkConverge(k::Float64) close() end -function makeKinemClVortPlots2D() +function makeKinemClVortPlots2D(mat = "", tt = "all") +# tt is a target time to make the plot for. It will choose the closest time available - mat, _ = DelimitedFiles.readdlm("resultsSummary", '\t', Float64, header=true) + if mat == "" + mat, _ = DelimitedFiles.readdlm("resultsSummary", '\t', Float64, header=true) + end t = mat[:,1] len = length(t) @@ -382,19 +416,29 @@ function makeKinemClVortPlots2D() h = mat[:,3] u = mat[:,4] + try + cd("Step Files") + catch + println(" Error: No 'Step Files' directory. Make sure wflag is enabled.") + return + end + dirvec = readdir() dirresults = map(x->(v = tryparse(Float64,x); typeof(v) == Nothing ? 0.0 : v),dirvec) - + if tt != "all" + ttidx = findmin(abs.(dirresults.-tt))[2] + end + #Determine axis limits dirmax = maximum(dirresults) - + cd("$(dirmax)") - + multsurfflag = 0 if isfile("boundv-1") == true #only 1 surface error("this plot function is only written for single surface") end - + tev, _ = DelimitedFiles.readdlm("tev", '\t', Float64, header=true) lev = try DelimitedFiles.readdlm("lev", '\t', Float64, header=true)[1] @@ -402,28 +446,35 @@ function makeKinemClVortPlots2D() zeros(0,3) end bv, _ = DelimitedFiles.readdlm("boundv", '\t', Float64, header=true) - + xminv = minimum([tev[:,2];lev[:,2];bv[:,2];]) zminv = minimum([tev[:,3];lev[:,3];bv[:,3];]) xmaxv = maximum([tev[:,2];lev[:,2];]) zmaxv = maximum([bv[:,2];tev[:,3];lev[:,3];bv[:,3];]) - - cd("..") - - if "infoPlots" in dirvec - rm("infoPlots", recursive=true) + + cd("../..") + if tt == "all" + dirvec = readdir() + if "infoPlots" in dirvec + rm("infoPlots", recursive=true) + end + mkdir("infoPlots") end - mkdir("infoPlots") - + + cd("Step Files") for i=1:length(dirresults) + if tt != "all" && ttidx != i + continue + end if dirresults[i] != 0 + dirstr="$(dirresults[i])" cd(dirstr) t_cur = dirresults[i] - + figure(figsize=(8,8)) - + subplot2grid((6, 2), (0, 0)) plot(t, alpha) plot([t_cur; t_cur], [-10000; 10000], "k-") @@ -484,19 +535,30 @@ function makeKinemClVortPlots2D() bv, _ = DelimitedFiles.readdlm("boundv", '\t', Float64, header=true) viewVort2D(tev, lev, bv) axis([xminv-1, xmaxv+1, zminv-1, zmaxv+1]) - + tight_layout() - savefig("../infoPlots/$(dirresults[i]).png") - close() + if tt == "all" + savefig("../../infoPlots/$(dirresults[i]).png") + close() + end cd("..") end end + cd("..") + if tt == "all" + println(" All plots saved under the 'infoPlots' directory.") + else + pause(0.0001) + show(block = false) + end end -function makeKinemVelVortPlots2D() +function makeKinemVelVortPlots2D(mat = "",tt = "all") - mat, _ = DelimitedFiles.readdlm("resultsSummary", '\t', Float64, header=true) + if mat == "" + mat, _ = DelimitedFiles.readdlm("resultsSummary", '\t', Float64, header=true) + end t = mat[:,1] len = length(t) @@ -505,19 +567,30 @@ function makeKinemVelVortPlots2D() h = mat[:,3] u = mat[:,4] + try + cd("Step Files") + catch + println(" Error: No 'Step File' directory. Make sure wflag is enabled.") + return + end + dirvec = readdir() dirresults = map(x->(v = tryparse(Float64,x); typeof(v) == Nothing ? 0.0 : v),dirvec) - + if tt != "all" + ttidx = findmin(abs(dirresults-tt))[2] + end + + #Determine axis limits dirmax = maximum(dirresults) - + cd("$(dirmax)") - + multsurfflag = 0 if isfile("boundv-1") == true #only 1 surface error("this plot function is only written for single surface") end - + tev, _ = DelimitedFiles.readdlm("tev", '\t', Float64, header=true) lev = try DelimitedFiles.readdlm("lev", '\t', Float64, header=true)[1] @@ -525,28 +598,33 @@ function makeKinemVelVortPlots2D() zeros(0,3) end bv, _ = DelimitedFiles.readdlm("boundv", '\t', Float64, header=true) - + xminv = minimum([tev[:,2];lev[:,2];bv[:,2];]) zminv = minimum([tev[:,3];lev[:,3];bv[:,3];]) xmaxv = maximum([tev[:,2];lev[:,2];]) zmaxv = maximum([bv[:,2];tev[:,3];lev[:,3];bv[:,3];]) - - cd("..") - - if "infoPlots" in dirvec - rm("infoPlots", recursive=true) + + cd("../..") + + dirvec = readdir() + if "infoplots" in dirvec + rm("infoplots", recursive=true) end - mkdir("infoPlots") - + mkdir("infoplots") + + cd("Step Files") for i=1:length(dirresults) + if tt != "all" && ttidx != i + continue + end if dirresults[i] != 0 dirstr="$(dirresults[i])" cd(dirstr) t_cur = dirresults[i] - + figure(figsize=(8,8)) - + subplot2grid((6, 2), (0, 0)) plot(t, alpha) plot([t_cur; t_cur], [-10000; 10000], "k-") @@ -607,15 +685,105 @@ function makeKinemVelVortPlots2D() bv, _ = DelimitedFiles.readdlm("boundv", '\t', Float64, header=true) viewVort2D(tev, lev, bv) axis([xminv-1, xmaxv+1, zminv-1, zmaxv+1]) - + tight_layout() - savefig("../infoPlots/$(dirresults[i]).png") + savefig("../../infoplots/$(dirresults[i]).png") close() cd("..") end end + cd("..") +end + +function subPlotForces(mat, xVar :: Int64, xVarLabel :: String, saveName :: String) + # Create Cl,Cd,Cm,LESP vs xVar plot + # mat = [t , alpha , h , u , LESP , cl , cd , cm] for each timestep + xVar = mat[:,xVar] ; + cl = mat[:,6] ; + cd = mat[:,7] ; + cm = mat[:,8] ; + LESP = mat[:,5] ; + + subplot(221) + plot(xVar,cl) ; + xlabel(xVarLabel) + ylabel("Cl") + + subplot(222) + plot(xVar,cd) ; + xlabel(xVarLabel) + ylabel("Cd") + + subplot(223) + plot(xVar,cm) ; + xlabel(xVarLabel) + ylabel("Cm") + + subplot(224) + plot(xVar,LESP) ; + xlabel(xVarLabel) + ylabel("LESP") + + tight_layout() + savefig(saveName) + close() end +## UI Plots +function initPlot(A,C,t,a,h,u) + # Create intitilization plots that shows airfoil and camber and motion parameters + + # Camber Plot + if A != 0 # airfoil defined + subplot(121,aspect = 1) + plot(A[:,1],A[:,2],color = "black") # full airfoil + plot(C[:,1],C[:,2], color = "blue", marker = "|") # Camber + ylim(-.5,.5) + xlabel("x") + ylabel("z") + title("Geometry") + end + + # Motion plots + if sum(a) != 0 # alphadef defined + subplot(122) + plot(t,h, color = "black") # height + plot(t,u, color = "blue") # velocity + plot(t,a, color = "red") # alpha + xlabel(L"$t^*$") + legend(["Height","Velocity","Alpha [deg]"]) + title("Motion") + end + pause(0.0001) + show(block = false) +end - +function plotMat(mat,pltCnt,plotx,ploty,matIdx) + # plot subplots (number determined by plotCnt) with an x-axis of plotx + # and a y axis of ploty. plotx/y are indices of the columns of mat. + + figure() + # find subplot input dimensions + if pltCnt <= 2 # plot on one row + subRow = 1 + subCol = pltCnt + else + subRow = 2 + subCol = ceil(pltCnt/2) + end + subDim = subRow*100 + subCol*10 + + for i = 1:pltCnt + subplot(subDim + i) + x = matIdx[lowercase(plotx)] + y = matIdx[lowercase(ploty[i])] + plot(mat[:,x],mat[:,y],"b-") + xlabel(plotx) + ylabel(ploty[i]) + title("$(ploty[i]) vs $plotx") + end + tight_layout() + pause(0.0001) + show(block = false) +end diff --git a/src/utils.jl b/src/utils.jl index 9c98fb45..1f404e4a 100755 --- a/src/utils.jl +++ b/src/utils.jl @@ -4,12 +4,22 @@ Clears all timestamp directories in the current folder """ function cleanWrite() + dirvec = readdir() - dirresults = map(x->(v = tryparse(Float64,x); typeof(v) == Nothing ? 0.0 : v),dirvec) - for i =1:length(dirresults) - rm("$(dirresults[i])", force=true, recursive=true) + if "Step Files" in dirvec + try + rm("Step Files", recursive=true) + catch + println(" ERROR: Unable to reset 'Step Files' directory") + end + else + dirvec = readdir() + dirresults = map(x->(v = tryparse(Float64,x); typeof(v) == Nothing ? 0.0 : v),dirvec) + for i =1:length(dirresults) + rm("$(dirresults[i])", force=true, recursive=true) + end + # rm("*~", force=true) end - rm("*~", force=true) end @@ -98,7 +108,7 @@ function simpleTrapz(y::Vector{T}, x::Vector{T}) where {T<:Real} end # Aerofoil camber calculation from coordinate file -function camber_calc(x::Vector,airfoil::String) +function camber_calc(x::Vector,airfoil::Union{String,SubString{String}}) #Determine camber and camber slope on airfoil from airfoil input file ndiv = length(x); @@ -147,7 +157,7 @@ function camber_thick_calc(x::Vector,coord_file::String) m = parse(Int, coord_file[5])/100. p = parse(Int, coord_file[6])/10. th = parse(Int, coord_file[7:8])/100. - + b1 = 0.2969; b2 = -0.1260; b3 = -0.3516; b4 = 0.2843; b5 = -0.1015 for i = 2:ndiv thick[i] = 5*th*(b1*sqrt(x[i]) + b2*x[i] + b3*x[i]^2 + b4*x[i]^3 + b5*x[i]^4) @@ -226,4 +236,3 @@ function camber_thick_calc(x::Vector,coord_file::String) end return thick, thick_slope,rho, cam, cam_slope end - diff --git a/src/xfoil/XfoilWrapper.jl b/src/xfoil/XfoilWrapper.jl new file mode 100644 index 00000000..d5d71ac1 --- /dev/null +++ b/src/xfoil/XfoilWrapper.jl @@ -0,0 +1,79 @@ +#= +XfoilWrapper.jl + Written by Matthew White + 4/2/2019 + + Input: config file by name of Config.txt + Output: Cn,Cm and polar file +=# +function xfoilWrapper(cnfg::String = "bin/config.txt") + # import config file + cnfg = DelimitedFiles.readdlm(cnfg, String) + # airfoil + uppercase(cnfg[1,2]) == "NACA" ? airfoil = cnfg[1,3] : airfoil = cnfg[1,2] + re = cnfg[2,2] # Reynolds + mach = cnfg[3,2] # Mach # + changePanels = uppercase(cnfg[4,2]) .== "TRUE" # Manually change panels (otherwise runs PANE) + panels = cnfg[5,2] # Panels to change to + changeIter = uppercase(cnfg[6,2]) .== "TRUE" # change the number of iterations + iter = cnfg[7,2] # iteration count to change to + # AoA range + aStart = parse(Float32, cnfg[8,2]) + aEnd = parse(Float32, cnfg[9,2]) + aStep = parse(Float32, cnfg[10,2]) + # XFOIL path + xPath = cnfg[11,2] + # Polar file path + polarPath = cnfg[12,2] + + # --- Open XFOIL Pipe --- + p=open(pipeline(`$xPath`),"r+") + run(`rm -rf $polarPath`) # remove any old polar file + # Write to pipe + length(airfoil) <= 5 ? write(p,"naca $airfoil\n") : write(p,"load $airfoil\n") + changePanels ? write(p,"ppar $panels\n\n") : write(p,"pane\n") + write(p,"oper\n") + write(p,"visc $re\n") + if changeIter write(p,"iter $iter\n") end + write(p,"mach $mach\n") + write(p,"pacc\n") + write(p,"$polarPath\n\n") + write(p,"aseq $aStart $aEnd $aStep\n") + write(p,"\nquit") + + close(p) + + while !isfile(polarPath) end # wait unitl polar file has been created + while filesize(polarPath) < 1000 end # wait until at least one value has been written + + # --- Reading the Polar File --- + prevSize = 0 + sec = 0 + n = 10 + while true # wait unit the polar file is done changing + global polar = DelimitedFiles.readdlm(polarPath,Float64,skipstart = 12) + if polar[end,1] > aEnd - aStep + println("Done writing polar file") + break #File is fully loaded + else #check every n seconds for stable file size + sleep(1) + sec += 1 + if prevSize == filesize(polarPath) && sec == n + println("Method could not converge for the final alphas") + println(" Maximum angle of attack at $(polar[end,1]) degrees") + break # after n seconds if filesize is constant then use current polar + elseif sec >= n + sec = 0 # reset sec every n seconds + end + prevSize = filesize(polarPath) + end + end + + # polar: alpha | Cl | CD | CDp | Cm | Top_Xtr | Bot_Xtr + # Cn = Cl*cos(a) + Cd*sin(a) + Cn = polar[:,2].*cosd.(polar[:,1]) + polar[:,3].*sind.(polar[:,1]) + Cm = polar[:,5] + + Cn,Cm + +end