diff --git a/Project.toml b/Project.toml index 02e5175a..e2cf27b9 100644 --- a/Project.toml +++ b/Project.toml @@ -17,12 +17,13 @@ KinematicCoordinateTransformations = "730d3219-0a85-48f9-b699-9f31f8913d09" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Meshes = "eacbb407-ea5a-433e-ab97-5258b1ca43fa" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" +Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192" [compat] AcousticMetrics = "0.7.0" CCBlade = "0.2.4" -FLOWMath = "0.3.3" +FLOWMath = "0.4.1" FillArrays = "1.11.0" FlexiMaps = "0.1.18" Format = "1.3.7" @@ -30,5 +31,6 @@ JuliennedArrays = "0.4.0" KinematicCoordinateTransformations = "0.4.1" Meshes = "0.46.0" StaticArrays = "1.5.19" +Statistics = "1.9.0" WriteVTK = "1.17.1" julia = "1.9.1" diff --git a/docs/Project.toml b/docs/Project.toml index 28de79fa..2652f35e 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -7,6 +7,7 @@ DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" FLOWMath = "6cb5d3fb-0fe8-4cc2-bd89-9fe0b19a99d3" FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549" +FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b" GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a" Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" diff --git a/docs/make.jl b/docs/make.jl index ea4e2d8b..14a1bebb 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -5,6 +5,7 @@ function doit() IN_CI = get(ENV, "CI", nothing)=="true" makedocs(sitename="AcousticAnalogies.jl", modules=[AcousticAnalogies], doctest=false, + root=@__DIR__, format=Documenter.HTML(prettyurls=IN_CI), pages=["Introduction"=>"index.md", "Guided Example"=>"guided_example.md", diff --git a/docs/src/assets/openfast_example_with_obs-downstream_view-loading_y-0to6000-compressed.gif b/docs/src/assets/openfast_example_with_obs-downstream_view-loading_y-0to6000-compressed.gif new file mode 100644 index 00000000..15580936 Binary files /dev/null and b/docs/src/assets/openfast_example_with_obs-downstream_view-loading_y-0to6000-compressed.gif differ diff --git a/docs/src/assets/openfast_example_with_obs-iso_view-loading_y-0to6000-compressed.gif b/docs/src/assets/openfast_example_with_obs-iso_view-loading_y-0to6000-compressed.gif new file mode 100644 index 00000000..c00bd57a Binary files /dev/null and b/docs/src/assets/openfast_example_with_obs-iso_view-loading_y-0to6000-compressed.gif differ diff --git a/docs/src/openfast_example.md b/docs/src/openfast_example.md index 435133ae..eb7d9e19 100644 --- a/docs/src/openfast_example.md +++ b/docs/src/openfast_example.md @@ -1,9 +1,9 @@ - ```@meta CurrentModule = AADocs ``` # Compact Formulation 1A OpenFAST Example +## Introduction This example loads a .out file generated by the popular aeroserohydroelastic solver [OpenFAST](https://github.com/OpenFAST/openfast), which is released by the U.S. National Renewable Energy Laboratory to simulate wind turbines, @@ -12,32 +12,29 @@ The example simulates the acoustic emissions of the 3.4MW land-based reference w turbine released by the International Wind Energy Agency. The OpenFAST model is available at https://github.com/IEAWindTask37/IEA-3.4-130-RWT. - We start by loading Julia dependencies, which are available in the General registry ```@example first_example -using AcousticAnalogies -using AcousticMetrics +using AcousticAnalogies: AcousticAnalogies +using AcousticMetrics: AcousticMetrics using ColorSchemes: colorschemes +using FillArrays: FillArrays, getindex_value using GLMakie -using Interpolations: linear_interpolation -using FLOWMath: FLOWMath -using KinematicCoordinateTransformations -using LinearAlgebra: × +using KinematicCoordinateTransformations: SteadyRotYTransformation using StaticArrays: @SVector -using Statistics nothing # hide ``` +## Inputs Next, we set the user-defined inputs: * number of blades, usually 3 for modern wind turbines * hub radius in m, it is specified in the ElastoDyn main input file of OpenFAST * blade spanwise grid in m and the corresponding chord, also in m. The two arrays are specified in the AeroDyn15 blade input file -* Observer location in the global coordinate frame (located at the rotor center, x points downwind, z points vertically up, and y points sideways). In this case we picked the IEC-prescribed location (turbine height on the ground) by specifying the hub heigth of 110 m. +* Observer location in the global coordinate frame (located at the rotor center, x points downwind, z points vertically up, and y points sideways). In this case we picked the IEC-prescribed location (turbine height on the ground) by specifying the hub height of 110 m. * Air density and speed of sound * Path to the OpenFAST .out file. The file must contain these channels: Time (always available), Wind1VelX from InflowWind, RotSpeed from ElastoDyn, Nodal outputs Fxl and Fyl from AeroDyn15. the file is available in the repo under test/gen_test_data/openfast_data ```@example first_example -num_blades = 3 +# num_blades = 3 Rhub = 2. BlSpn = [0.0000e+00, 2.1692e+00, 4.3385e+00, 6.5077e+00, 8.6770e+00, 1.0846e+01, 1.3015e+01, 1.5184e+01, 1.7354e+01, 1.9523e+01, 2.1692e+01, 2.3861e+01, 2.6031e+01, 2.8200e+01, 3.0369e+01, 3.2538e+01, @@ -48,7 +45,6 @@ Chord = [2.600e+00, 2.645e+00, 3.020e+00, 3.437e+00, 3.781e+00, 4.036e+00, 4.201 2.986e+00, 2.770e+00, 2.581e+00, 2.412e+00, 2.266e+00, 2.142e+00, 2.042e+00, 1.964e+00, 1.909e+00, 1.870e+00, 1.807e+00, 1.666e+00, 1.387e+00, 9.172e-01, 1.999e-01] file_path = joinpath(@__DIR__, "..", "..", "test", "gen_test_data", "openfast_data", "IEA-3.4-130-RWT.out") -@show @__DIR__ HH = 110. # m RSpn = BlSpn .+ Rhub x0 = @SVector [HH .+ RSpn[end], 0.0, -HH] @@ -57,195 +53,55 @@ c0 = 340.0 # m/s nothing # hide ``` -After that, we start with the same computations specified in the other examples. Refer to those before jumping onto this example. +For the monopole/thickness noise, we need the cross-sectional area at each radial station. +If we know the cross-sectional area per chord squared, we can find the cross-sectional area this way: ```@example first_example -# Compute the mid-section locations in radial coordinates, m -radii = 0.5.*(RSpn[begin:end-1] .+ RSpn[begin+1:end]) -# Compute the length of each section along blade span, m -dradii = RSpn[begin+1:end] .- RSpn[begin:end-1] -# Compute the blade angles -θs = 2*pi/num_blades.*(0:(num_blades-1)) -# Create a linear interpolation object to interpolate chord onto the radial mid-points -itp_chord = linear_interpolation(RSpn, Chord) -# Perform interpolation -chord = itp_chord(radii) # Cross-sectional area of each element in m**2. This is taking a bit of a shortcut—the value of `cs_area_over_chord_squared` does not actually correspond to the IEAWindTask37 turbine blade. cs_area_over_chord_squared = 0.064 -cs_area = cs_area_over_chord_squared.*chord.^2 - -# # Code to parse the data from the OpenFAST .out file -# # Function to parse a line of data, converting strings to floats -# function parse_line(line) -# # Split the line by whitespace and filter out any empty strings -# elements = filter(x -> !isempty(x), split(line)) -# # Convert elements to Float64 -# return map(x -> parse(Float64, x), elements) -# end -# -# # Initialize an empty array to store the data -# data = [] -# -# # Open the file and read the data, skipping the first 8 lines -# open(file_path) do file -# # Skip the first 8 lines (header and description) -# for i in 1:8 -# readline(file) -# end -# -# # Read the rest of the lines and parse them -# for line in eachline(file) -# push!(data, parse_line(line)) -# end -# end -# -# # Convert the data to an array of arrays (matrix) -# data = reduce(hcat, data) -# time = data[1, :] -# avg_wind_speed = mean(data[2, :]) -# sim_length_s = time[end] - time[1] # s -# @show length(time) -# -# # Reopen the file and read the lines -# lines = open(file_path) do f -# readlines(f) -# end -# -# # Find the index of the line that contains the column headers -# header_index = findfirst(x -> startswith(x, "Time"), lines) -# -# # Extract the headers -# headers = split(lines[header_index], '\t') -# -# id_b1_Fx = findfirst(x -> x == "AB1N001Fxl", headers) -# id_b2_Fx = findfirst(x -> x == "AB2N001Fxl", headers) -# id_b3_Fx = findfirst(x -> x == "AB3N001Fxl", headers) -# id_b1_Fy = findfirst(x -> x == "AB1N001Fyl", headers) -# id_b2_Fy = findfirst(x -> x == "AB2N001Fyl", headers) -# id_b3_Fy = findfirst(x -> x == "AB3N001Fyl", headers) -# id_rot_speed = findfirst(x -> x == "RotSpeed", headers) -# n_elems = length(radii) -# Fx_b1_locs = data[id_b1_Fx:id_b1_Fx+n_elems,:] -# Fy_b1_locs = data[id_b1_Fy:id_b1_Fy+n_elems,:] -# Fx_b2_locs = data[id_b2_Fx:id_b2_Fx+n_elems,:] -# Fy_b2_locs = data[id_b2_Fy:id_b2_Fy+n_elems,:] -# Fx_b3_locs = data[id_b3_Fx:id_b3_Fx+n_elems,:] -# Fy_b3_locs = data[id_b3_Fy:id_b3_Fy+n_elems,:] -# -# # Reinterpolate onto the mid-sections -# Fx_b1 = Array{Float64}(undef, 29, 6001) -# Fx_b2 = Array{Float64}(undef, 29, 6001) -# Fx_b3 = Array{Float64}(undef, 29, 6001) -# Fy_b1 = Array{Float64}(undef, 29, 6001) -# Fy_b2 = Array{Float64}(undef, 29, 6001) -# Fy_b3 = Array{Float64}(undef, 29, 6001) -# for j in axes(Fx_b1_locs, 2) -# itp = LinearInterpolation(RSpn, Fx_b1_locs[:, j], extrapolation_bc=Line()) -# Fx_b1[:, j] = itp(radii) -# end -# for j in axes(Fx_b2_locs, 2) -# itp = LinearInterpolation(RSpn, Fx_b2_locs[:, j], extrapolation_bc=Line()) -# Fx_b2[:, j] = itp(radii) -# end -# for j in axes(Fx_b3_locs, 2) -# itp = LinearInterpolation(RSpn, Fx_b3_locs[:, j], extrapolation_bc=Line()) -# Fx_b3[:, j] = itp(radii) -# end -# for j in axes(Fy_b1_locs, 2) -# itp = LinearInterpolation(RSpn, Fy_b1_locs[:, j], extrapolation_bc=Line()) -# Fy_b1[:, j] = itp(radii) -# end -# for j in axes(Fy_b2_locs, 2) -# itp = LinearInterpolation(RSpn, Fy_b2_locs[:, j], extrapolation_bc=Line()) -# Fy_b2[:, j] = itp(radii) -# end -# for j in axes(Fy_b3_locs, 2) -# itp = LinearInterpolation(RSpn, Fy_b3_locs[:, j], extrapolation_bc=Line()) -# Fy_b3[:, j] = itp(radii) -# end +cs_area = cs_area_over_chord_squared .* Chord.^2 +nothing # hide +``` + +Next, we'll use the [`read_openfast_file`](@ref) function to read the OpenFAST output file. +This will read the data in the file, but also do a bit of processing necessary for an acoustic prediction. +Specifically, it will... + + * interpolate the cross-sectional area and loading from the blade element interfaces to the cell centers + * use second-order finite differences to differentiate the loading with respect to time. + * average the freestream velocity and RPM (if `average_freestream_vel` or `average_omega` keyword arguments are `true`) + +The [`read_openfast_file` doc string](@ref read_openfast_file) has more information. +The output of `read_openfast_file` is a `OpenFASTData` `struct` that has fields like `time`, `omega`, `axial_loading`, etc. that are read from the output file, and also fields like `radii_mid`, `circum_loading_mid_dot` that are created after the output file is read. +Check out the [`OpenFASTData` doc string](@ref OpenFASTData) for a list of all the fields. + +```@example first_example # Read the data from the file and create an `OpenFASTData` object, a simple struct with fields like `time`, `omega`, `axial_loading`, etc. -data = AcousticAnalogies.read_openfast_file(file_path) - -# Interpolate the loading data into the mid-sections. -# Start with the axial loading. -# data.axial_loading has size (num_times, num_radial, num_blades). -# Since we're interpolating onto the mid-sections, the interpolated data will have size (num_times, num_radial-1, num_blades). -axial_loading = data.axial_loading -axial_loading_mid = Array{eltype(axial_loading)}(undef, size(axial_loading, 1), size(axial_loading, 2)-1, size(axial_loading, 3)) - -# Now start doing the interpolation. -for b in 1:size(axial_loading, 3) - for t in 1:size(axial_loading, 1) - itp_loading = linear_interpolation(RSpn, axial_loading[t, :, b]) - axial_loading_mid[t, :, b] .= itp_loading(radii) - end -end +data = AcousticAnalogies.read_openfast_file(file_path, RSpn, cs_area; average_freestream_vel=true, average_omega=true) +``` -# Do the same thing for the circumferential loading. -circum_loading = data.circum_loading -circum_loading_mid = Array{eltype(circum_loading)}(undef, size(circum_loading, 1), size(circum_loading, 2)-1, size(circum_loading, 3)) -for bidx in 1:size(circum_loading, 3) - for tidx in 1:size(circum_loading, 1) - itp_loading = linear_interpolation(RSpn, circum_loading[tidx, :, bidx]) - circum_loading_mid[tidx, :, bidx] .= itp_loading(radii) - end -end +We can get the averaged rotation rate value from the `OpenFASTData` `struct` this way: -# New thought: F1A is a function of the time-derivative of the loading. -# How should we get that? -# We could use finite differences. -# I guess we'll do that. -# But is the time step constant? -# Let's find out. -# @show extrema((data.time[2:end] .- data.time[1:end-1]) .- (data.time[2] - data.time[1])) -dts = data.time[2:end] - data.time[1:end-1] -@show all(dts .≈ dts[1]) -dt = dts[1] - -# Looks like it is uniformly spaced. - -# Now use second-order finite differences to differentiate. -# We'll create a simple function to do that for us. -function finite_diff_2nd_order(f, h) - # Create an array the same size as `f`. - fdot = similar(f) - - @views begin - # These stencils are in Tannehill, Anderson, Pletcher, "Computational Fluid Mechanics and Heat Transfer," 2nd edition, page 50. - # First do the interior points. - fdot[begin+1:end-1, :, :] .= (f[begin+2:end, :, :] .- f[begin:end-2, :, :]) ./ (2*h) - - # Then the left boundary. - fdot[begin, :, :] .= (-3 .* f[begin] .+ 4 .* f[begin+1] .- f[begin+2]) ./ (2*h) - - # Then the right boundary. - fdot[end, :, :] .= (3 .* f[end] .- 4 .* f[end-1] .+ f[end-2]) ./ (2*h) - end - - return fdot -end +```@example first_example +omega_avg = getindex_value(data.omega) +@show omega_avg +nothing # hide +``` -# Now use that function to differentiate the axial and circumferential loading. -axial_loading_mid_dot = finite_diff_2nd_order(axial_loading_mid, dt) -circum_loading_mid_dot = finite_diff_2nd_order(circum_loading_mid, dt) +(When averaging rotation rate or freestream velocity, `read_openfast_file` uses a [`Fill`](https://juliaarrays.github.io/FillArrays.jl/stable/#FillArrays.Fill) `struct` from the [`FillArrays.jl`](https://github.com/JuliaArrays/FillArrays.jl) package to lazily represent the average `omega` value as a length-`num_times` `Vector`, and [`getindex_value`](https://juliaarrays.github.io/FillArrays.jl/stable/#FillArrays.getindex_value) is a function from `FillArrays.jl` that returns that single averaged value. +Could have also just indexed the `data.omega` array at the first value, or last, etc..) -# Extract the mean rotor speed -omega_rpm = mean(data.omega) -@show omega_rpm +```@example first_example +@show data.omega[1] data.omega[8] data.omega[end] nothing # hide ``` -Once we are here we have parsed the OpenFAST output file and are ready to run the F1A code. -Before that, let's plot the unsteady loading acting on the wind turbine blade. +Before we actually try an acoustic prediction, let's have a look at the loading. +We'll use the Makie plotting package to make the plots, and only plot 1 out of every 500 time steps (as seen in the `for tidx` line): ```@example first_example -# Let's plot the unsteady loading 1 of every 500 timesteps -# x-axis is the span position (mid-sections) -# times are indicated by the colorbar on the right of the plot. -# @assert size(Fx_b1) == size(Fx_b2) == size(Fx_b3) == size(Fy_b1) == size(Fy_b2) == size(Fy_b3) -# ntimes_loading = size(Fx_b1, 2) -ntimes_loading = size(axial_loading_mid, 1) +ntimes_loading = size(data.axial_loading_mid, 1) fig = Figure() ax11 = fig[1, 1] = Axis(fig, xlabel="Span Position (m)", ylabel="Fx (N/m)", title="blade 1") ax21 = fig[2, 1] = Axis(fig, xlabel="Span Position (m)", ylabel="Fy (N/m)") @@ -253,28 +109,29 @@ ax12 = fig[1, 2] = Axis(fig, xlabel="Span Position (m)", ylabel="Fx (N/m)", titl ax22 = fig[2, 2] = Axis(fig, xlabel="Span Position (m)", ylabel="Fy (N/m)") ax13 = fig[1, 3] = Axis(fig, xlabel="Span Position (m)", ylabel="Fx (N/m)", title="blade 3") ax23 = fig[2, 3] = Axis(fig, xlabel="Span Position (m)", ylabel="Fy (N/m)") -bpp = 60/omega_rpm/num_blades +num_blades = data.num_blades colormap = colorschemes[:viridis] time = data.time sim_length_s = time[end] - time[begin] for tidx in 1:500:ntimes_loading cidx = (time[tidx] - time[1])/sim_length_s - l1 = lines!(ax11, radii, axial_loading_mid[tidx,:,1], label ="b1", color=colormap[cidx]) - l1 = lines!(ax12, radii, axial_loading_mid[tidx,:,2], label ="b2", color=colormap[cidx]) - l1 = lines!(ax13, radii, axial_loading_mid[tidx,:,3], label ="b3", color=colormap[cidx]) - l2 = lines!(ax21, radii, circum_loading_mid[tidx,:,1], label ="b1", color=colormap[cidx]) - l2 = lines!(ax22, radii, circum_loading_mid[tidx,:,2], label ="b2", color=colormap[cidx]) - l2 = lines!(ax23, radii, circum_loading_mid[tidx,:,3], label ="b3", color=colormap[cidx]) + l1 = lines!(ax11, data.radii_mid, data.axial_loading_mid[tidx,:,1], label ="b1", color=colormap[cidx]) + l1 = lines!(ax12, data.radii_mid, data.axial_loading_mid[tidx,:,2], label ="b2", color=colormap[cidx]) + l1 = lines!(ax13, data.radii_mid, data.axial_loading_mid[tidx,:,3], label ="b3", color=colormap[cidx]) + l2 = lines!(ax21, data.radii_mid, data.circum_loading_mid[tidx,:,1], label ="b1", color=colormap[cidx]) + l2 = lines!(ax22, data.radii_mid, data.circum_loading_mid[tidx,:,2], label ="b2", color=colormap[cidx]) + l2 = lines!(ax23, data.radii_mid, data.circum_loading_mid[tidx,:,3], label ="b3", color=colormap[cidx]) end -linkxaxes!(ax21, ax11) linkxaxes!(ax21, ax11) linkxaxes!(ax12, ax11) linkxaxes!(ax22, ax11) linkxaxes!(ax13, ax11) linkxaxes!(ax23, ax11) + linkyaxes!(ax12, ax11) linkyaxes!(ax13, ax11) + linkyaxes!(ax22, ax21) linkyaxes!(ax23, ax21) @@ -288,18 +145,20 @@ hideydecorations!(ax23, grid=false) cbar = fig[:, 4] = Colorbar(fig; limits=(time[begin], time[end]), colormap=:viridis, label="time (sec)") -save(joinpath(@__DIR__, "Fx_t-all_time.png"), fig) +save(joinpath(@__DIR__, "openfast_example_loading.png"), fig) nothing # hide ``` -![Axial Loading](Fx_t-all_time.png) +![Loading](openfast_example_loading.png) + +The x axis of each subplot is the radial position along the blade, from hub to tip. +The top three plots show the axial loading, bottom three the circumferential, and there's one column for each blade. +And the colorbar indicates the simulation time. +The plot shows significant unsteadiness, which is cool to see. -Let's also see what the time derivative of the loading looks like. +We can also plot the loading time derivative in a similar form: ```@example first_example -# Let's plot the unsteady loading time derivative 1 of every 500 timesteps -# x-axis is the span position (mid-sections) -# times are indicated by the colorbar on the right of the plot. -ntimes_loading = size(axial_loading_mid_dot, 1) +ntimes_loading = size(data.axial_loading_mid_dot, 1) fig = Figure() ax11 = fig[1, 1] = Axis(fig, xlabel="Span Position (m)", ylabel="∂Fx/∂t (N/(m*s))", title="blade 1") ax21 = fig[2, 1] = Axis(fig, xlabel="Span Position (m)", ylabel="∂Fy/∂t (N/(m*s))") @@ -307,28 +166,29 @@ ax12 = fig[1, 2] = Axis(fig, xlabel="Span Position (m)", ylabel="∂Fx/∂t (N/( ax22 = fig[2, 2] = Axis(fig, xlabel="Span Position (m)", ylabel="∂Fy/∂t (N/(m*s))") ax13 = fig[1, 3] = Axis(fig, xlabel="Span Position (m)", ylabel="∂Fx/∂t (N/(m*s))", title="blade 3") ax23 = fig[2, 3] = Axis(fig, xlabel="Span Position (m)", ylabel="∂Fy/∂t (N/(m*s))") -bpp = 60/omega_rpm/num_blades +num_blades = data.num_blades colormap = colorschemes[:viridis] time = data.time sim_length_s = time[end] - time[begin] for tidx in 1:500:ntimes_loading cidx = (time[tidx] - time[1])/sim_length_s - l1 = lines!(ax11, radii, axial_loading_mid_dot[tidx,:,1], label ="b1", color=colormap[cidx]) - l1 = lines!(ax12, radii, axial_loading_mid_dot[tidx,:,2], label ="b2", color=colormap[cidx]) - l1 = lines!(ax13, radii, axial_loading_mid_dot[tidx,:,3], label ="b3", color=colormap[cidx]) - l2 = lines!(ax21, radii, circum_loading_mid_dot[tidx,:,1], label ="b1", color=colormap[cidx]) - l2 = lines!(ax22, radii, circum_loading_mid_dot[tidx,:,2], label ="b2", color=colormap[cidx]) - l2 = lines!(ax23, radii, circum_loading_mid_dot[tidx,:,3], label ="b3", color=colormap[cidx]) + l1 = lines!(ax11, data.radii_mid, data.axial_loading_mid_dot[tidx,:,1], label ="b1", color=colormap[cidx]) + l1 = lines!(ax12, data.radii_mid, data.axial_loading_mid_dot[tidx,:,2], label ="b2", color=colormap[cidx]) + l1 = lines!(ax13, data.radii_mid, data.axial_loading_mid_dot[tidx,:,3], label ="b3", color=colormap[cidx]) + l2 = lines!(ax21, data.radii_mid, data.circum_loading_mid_dot[tidx,:,1], label ="b1", color=colormap[cidx]) + l2 = lines!(ax22, data.radii_mid, data.circum_loading_mid_dot[tidx,:,2], label ="b2", color=colormap[cidx]) + l2 = lines!(ax23, data.radii_mid, data.circum_loading_mid_dot[tidx,:,3], label ="b3", color=colormap[cidx]) end -linkxaxes!(ax21, ax11) linkxaxes!(ax21, ax11) linkxaxes!(ax12, ax11) linkxaxes!(ax22, ax11) linkxaxes!(ax13, ax11) linkxaxes!(ax23, ax11) + linkyaxes!(ax12, ax11) linkyaxes!(ax13, ax11) + linkyaxes!(ax22, ax21) linkyaxes!(ax23, ax21) @@ -342,143 +202,227 @@ hideydecorations!(ax23, grid=false) cbar = fig[:, 4] = Colorbar(fig; limits=(time[begin], time[end]), colormap=:viridis, label="time (sec)") -save(joinpath(@__DIR__, "Fx_t-all_time-dot.png"), fig) +save(joinpath(@__DIR__, "openfast_example_loading_dot.png"), fig) +nothing # hide +``` +![Loading Time Derivative](openfast_example_loading_dot.png) + +## Constructing the Source Elements +Now, the next step is to turn the OpenFAST data into source elements. +This step is pretty easy, since there is a function called [`f1a_source_elements_openfast`](@ref f1a_source_elements_openfast) that takes the `OpenFASTData` `struct` and a few other parameters and will create the source elements for us. +But first we need to think about the coordinate system we'd like our source elements to be in. +Eventually, we want the turbine blades to be rotating about the positive x axis, with the freestream velocity pointing in the positive x axis. +But there are two things we need to account for to make that happen: + + * To do a proper noise prediction, AcousticAnalogies.jl needs the source elements' motion to be defined in a reference frame relative to the ambient fluid, not the ground. + Put another way, we need to manipulate the source elements in a way so that it appears that there is no freestream velocity—that the ambient fluid is stationary. + So instead of having blade elements that are only rotating about a fixed hub position relative to the ground in a freestream pointed in the positive x direction, we will have the blade elements translate in the *negative* x direction as they rotate, with no freestream velocity. + * The `f1a_source_elements_openfast` routine puts the source elements in the Standard AcousticAnalogies.jl Reference Frame™, where the source elements + + * begin with the hub (rotation center) at `x = 0` at source time `t = 0` + * rotate about either the positive x or negative x axis (depending on the value of the `positive_x_rotation` argument). + * translate in the positive x direction, + +So, to make all this work, we'll initially have the source elements translate in the positive x direction (as required by `f1a_source_elements_openfast`) and rotate about the *negative* x axis. +Then we'll rotate the source elements 180° about the y axis, which will mean they will be translating in the negative x axis, rotating about the positive y axis, just like what we intend. + +So, here's the first step: create the source elements from the `OpenFASTData` `struct`, where they'll be rotating about the negative x axis, translating along the positive x axis. + +```@example first_example +positive_x_rotation = false +ses_before_roty = AcousticAnalogies.f1a_source_elements_openfast(data, rho, c0, positive_x_rotation) +nothing # hide +``` + +The `f1a_source_elements_openfast` returns an array of [`CompactF1ASourceElement`] `structs`. +The array is of size `(num_times, num_radial, num_blades)`, where `ses[i, j, k]` refers to the source element of the `i`th time step, `j`th radial position, and `k`th blade: + +```@example first_example +@show size(ses_before_roty) +nothing # hide +``` + +Now we'll rotate each source element 180° about the positive y axis. + +```@example first_example +# Create the object from KinematicCoordinateTransformations.jl defining the 180° rotation about the y axis. +rot180degy = SteadyRotYTransformation(0, 0, pi) + +# Now rotate the source elements. +ses = rot180degy.(ses_before_roty) + +# Could have combined all that in one line, i.e., +# ses = rot180degy.(AcousticAnalogies.f1a_source_elements_openfast(data, rho, c0, positive_x_rotation)) +nothing # hide +``` + +## Defining the Observer +The last thing we need before we can perform the noise prediction is an acoustic observer. +The observer is just the computational equivalent of the microphone. +In this case we picked the IEC-prescribed location (turbine height on the ground) by specifying the hub height of 110 m. +So we need the observer to be 110 m below the hub. +We'll also have the observer positioned downstream of the turbine rotation plane by a certain amount. + +```@example first_example +x0_obs = @SVector [HH + RSpn[end], 0.0, -HH] +nothing # hide +``` + +(The `@SVector` macro creates a statically-size vector using the [StaticArrays.jl](https://github.com/JuliaArrays/StaticArrays.jl) package, which is good for performance but not required.) + +Now, just like with the source elements, we need to define the motion of the observer relative to the fluid, not the ground. +So, we'll use the same trick that we used with the source elements: have the observer translate in the negative x direction to account for the freestream velocity that's pointed in the positive x direction: + +```@example first_example +# Get the average freestream velocity from the OpenFAST data. +v_avg = getindex_value(data.v) + +# Create a vector defining the velocity of the observer. +v_obs = @SVector [-v_avg, 0, 0] +nothing # hide +``` + +Since the observer is moving, its position is obviously changing. +So the `x0_obs` will be the position of the observer at the start of the simulation, at the first source time level of the source elements. +We can get that first source time level this way: + +```@example first_example +t0_obs = data.time[1] +nothing # hide +``` + +Now we have enough information to create the observer object: + +```@example first_example +obs = AcousticAnalogies.ConstVelocityAcousticObserver(t0_obs, x0_obs, v_obs) +nothing # hide +``` + +That says that we want our observer to start at the location `x0_obs` at time `t0_obs`, and then move with constant velocity `v_obs` forever after. +After creating the observer, we can query its location at any time value after this way: + +```@example first_example +@show obs(t0_obs) # should be equal to `x0_obs` +@show obs(t0_obs + 1) +nothing # hide +``` + +## Visualization with VTK Files +That was a lot. +How will we know we did all that correctly? +The answer is: write out the source elements and observer we just created to VTK files, and then visualize them with our favorite visualization software (ParaView at the moment). + +The function we want is [`to_paraview_collection`](@ref). +Using it is simple: + +```@example first_example +AcousticAnalogies.to_paraview_collection("openfast_example_with_obs", (ses,); observers=(obs,)) +nothing # hide +``` + +(This form of `to_paraview_collection` expects multiple arrays of source elements and multiple observers. +But here we just have one array of source elements (`ses`) and one observer (`obs`), so we wrap each in a single-entry tuple, i.e., `(ses,)` and `(obs,)`.) + +That will write out a bunch of VTK files showing the motion of the source elements and the observer, all starting with the `name` argument to the function (`openfast_example_with_obs` here). +The one to focus on is `openfast_example_with_obs.pvd`, a [ParaView data file](https://www.paraview.org/Wiki/ParaView/Data_formats#PVD_File_Format) that describes how all the many VTK files that `to_paraview_collection` writes out fit together. + +The VTK files for the source elements will also contain all the data defined in the source element `struct`s (the loading, cross-sectional area, etc.). +That's really handy for checking that the loading is in the correct direction (remember, it needs to be the loading on the fluid, i.e. exactly opposite the loading on the blades). + +To that end, here's an animation of the blades and observer, with the blades colored by the loading per unit span in the y direction: + +![LoadingY, Iso](assets/openfast_example_with_obs-iso_view-loading_y-0to6000-compressed.gif) + +Things look pretty good: the observer (i.e. the gray sphere) and the blades are all translating in the negative x direction, and the blades are rotating about the positive x axis. +(The gray smearing along the path of the observer is an artifact of the compression process the `gif` went through to make the file smaller.) +The y-component of the loading also appears to be in the correction direction: for a wind turbine, we'd expect the loading on the fluid to oppose the motion of the blade in the circumferential direction, which is what the animation shows. + +One thing that is troubling about the previous animation is the location of the observer relative to the blades in the y direction. +Since the rotor hub starts at the origin and moves along the negative x axis, and since y component of the observer position is always zero, the observer should only be offset in the x and z directions relative to the hub path. +That's hard to see in the previous animation, but if we switch our perspective to be looking directly downstream (i.e., looking in the positive-x direction), everything appears as it should be: + +![LoadingY, Downstream](assets/openfast_example_with_obs-downstream_view-loading_y-0to6000-compressed.gif) + +## Noise Prediction +Now we're finally ready to do a noise prediction! +The relevant function for that is [`noise`](@ref), which takes in a source element and observer and returns an [`F1AOutput`](@ref) `struct`, representing the acoustic pressure experienced by the observer due to the source: + +```@example first_example +apth = AcousticAnalogies.noise.(ses, Ref(obs)) nothing # hide ``` -![Axial Loading Time Derivative](Fx_t-all_time-dot.png) -Perfect, we are now ready to run the core of AcousticAnalogies.jl +Notice that we used a `.` after the `noise` function, which [broadcasts](https://docs.julialang.org/en/v1/manual/arrays/#Broadcasting) the `noise` call over all source element-observer combinations. +(The `Ref(obs)` makes the single observer `struct` act as a scalar during broadcasting, meaning the same observer object is passed to each `noise` call.) +Because of the broadcasting, `apth` is an `Array` of `F1AOutput` `structs` with the same size as `ses`: + +```@example first_example +@show size(ses) size(apth) +nothing # hide +``` +We now have a noise prediction for each of the individual source elements in `ses` at the acoustic observer obs—specifically, `apth[i, j, k]` represents the acoustic pressure for the `i` time step, `j` radial location, and `k` blade. +What we ultimately want is the total noise prediction at `obs`—we want to add all the acoustic pressures for each time level in `apth` together. +But we can't add them directly, yet, since the observer times, the time at which each source's noise reaches the observer, are not all the same. +What we need to do is first interpolate the acoustic pressure time history of each source onto a common observer time grid, and then add them up. +We'll do this using the [`combine`](@ref) function. +First, we need to decide on the length of the observer time grid and how many points it will contain. +If the motion and loading of the blades was steady, then one blade pass would be sufficient, but for this example that is not the case, so we'll use a longer observer time: ```@example first_example -# To do F1A correctly, we need to put all source elements in a coordinate system that -# moves with the fluid, i.e. one in which the fluid appears stationary. -# So, to do that, we have the blades translating in the -# negative x direction at the average horizontal wind speed. -avg_wind_speed = mean(data.v) -v = -avg_wind_speed # m/s -omega = omega_rpm * 2*pi/60 # rad/s - -# some reshaping, ses[i, j, k] holds the CompactSourceElement at src_time[i], radii[j], and blade number k -θs = reshape(θs, 1, 1, :) -radii = reshape(radii, 1, :, 1) -dradii = reshape(dradii, 1, :, 1) -cs_area = reshape(cs_area, 1, :, 1) -src_times = reshape(time, :, 1, 1) # This isn't really necessary. -# fx = cat(transpose(Fx_b1), transpose(Fx_b2), transpose(Fx_b3), dims=3) -# fc = cat(transpose(Fy_b1), transpose(Fy_b2), transpose(Fy_b3), dims=3) -# The `axial_loading_mid` and `circum_loading_mid` have shape `(num_times, num_radial-1, num_blades)`, so no reshaping necessary. -fx = axial_loading_mid -fc = circum_loading_mid - -# AcousticAnalogies.jl requires the loading *on the fluid*, not on the blade/lifting surface/body that is commonly calculated by aerodynamic prediction codes. -# The freestream velocity is in the positive x direction, so, for a wind turbine, we'd expect the axial loading on the blades to also be in the positive x direction. -# The plots of the loading above show that the axial loading is positive. -# The axial loading on the fluid is opposite that, so we'll need to switch the sign on the axial loading. - -# Next, we'll think about the circumferential loading. -# We will be rotating the blades about the positive x axis. -# The `CompactSourceElement` constructor assumes that the blade is initially aligned with the y axis (and then rotated by the angle `θ` about the x axis). -# So the blade will be moving in the positive z direction, and the circumferential loading on the blade will be in negative z direction, matching the plot of `circum_loading` above. -# But we want the loading on the fluid, so we'll also need to switch the sign on the `circum_loading` array. -ses = CompactSourceElement.(rho, c0, radii, θs, dradii, cs_area, -fx, 0.0, -fc, src_times) - -# t0 = 0.0 # Time at which the angle between the source and target coordinate systems is equal to offest. -# offset = 0.0 # Angular offset between the source and target coordinate systems at t0. -# # steady rotation around the x axis -# rot_trans = SteadyRotXTransformation(t0, omega, offset) -# -# # orient the rotation axis of the blades as it is the global frame -# rot_axis = @SVector [1.0, 0.0, 0.0] # rotation axis aligned along global x-axis -# blade_axis = @SVector [0.0, 0.0, 1.0] # blade 1 pointing up, along z-axis -# global_trans = ConstantLinearMap(hcat(rot_axis, blade_axis, rot_axis×blade_axis)) -# -# # blade to move with the appropriate forward velocity, and -# # start from the desired location in the global reference frame -# y0_hub = @SVector [0.0, 0.0, 0.0] # Position of the hub at time t0 -# v0_hub = SVector{3}(v.*rot_axis) # Constant velocity of the hub in the global reference frame -# const_vel_trans = ConstantVelocityTransformation(t0, y0_hub, v0_hub) -# -# # combine these three transformations into one, and then use that on the SourceElements -# trans = compose.(src_times, Ref(const_vel_trans), compose.(src_times, Ref(global_trans), Ref(rot_trans))) -# -# # trans will perform the three transformations from right to left (rot_trans, global_trans, const_vel_trans) -# ses = ses .|> trans -# -# # The ses object now describes how each blade element source is moving through the global reference -# # frame over the time src_time. As it does this, it will emit acoustics that can be sensed by an acoustic observer -# # (a human, or a microphone). The exact "amount" of acoustics the observer will experience depends -# # on the relative location and motion between each source and the observer. -# -# # This creates an acoustic observer moving with constant velocity v0_hub that is at location `x0` at time `t0`. -# obs = ConstVelocityAcousticObserver(t0, x0, v0_hub) -# -# # Now, in order to perform the F1A calculation, -# # we need to know when each acoustic disturbance emitted -# # by the source arrives at the observer. This is referred -# # to an advanced time calculation, and is done this way: -# apth = f1a.(ses, Ref(obs)) -# -# # We now have a noise prediction for each of the individual source elements in ses at the acoustic observer obs. -# # What we ultimately want is the total noise prediction at obs—we want to add all the acoustic pressures in apth together. -# # But we can't add them directly, yet, since the observer times are not all the same. What we need to do -# # is first interpolate the apth of each source onto a common observer time grid, and then add them up. -# # We'll do this using the AcousticAnalogies.combine function. -# period = 2*pi/omega -# bpp = period/num_blades # blade passing period -# obs_time_range = sim_length_s/60*omega_rpm*bpp -# -# # Note that we need to be careful to avoid extrapolation in the `combine` calculation. -# # That won't happen in this case, since obs_time_range/sim_length_s is 1/3, so the observer time -# # range is much less than the source time range. -# # The observer time range is 1/3 of the source time range, and we're using the same number of -# # simulation times, so that means the observer time step is 1/3 that of the source time step. -# num_obs_times = length(time) -# apth_total = combine(apth, obs_time_range, num_obs_times, 1) -# # The loading data is unsteady, so we may need to be careful to window the time history -# # to avoid problems with discontinuities going from the begining/end of the pressure time history. -# nothing # hide +rev_period = 2*pi/omega_avg +bpp = rev_period/num_blades # blade passing period +omega_rpm = omega_avg * 60/(2*pi) +obs_time_range = sim_length_s/60*omega_rpm*bpp +num_obs_times = length(data.time) +nothing # hide ``` -We can now have a look at the total acoustic pressure time history at the observer: +So that says that we'll have an output observer time length of `obs_time_range` with `num_obs_times` points. +Note that we need to be careful to avoid extrapolation in the `combine` calculation, which will happen if the observer time specified via the `obs_time_range` and `num_obs_times` arguments to `combine` extends past the times contained in the `apth` array. +That won't happen in this case, since `obs_time_range/sim_length_s` is 1/3, so the observer time range is much less than the source time range. + +Now we call `combine` to get the total acoustic pressure time history: ```@example first_example -# fig = Figure() -# ax1 = fig[1, 1] = Axis(fig, xlabel="time, s", ylabel="monopole, Pa") -# ax2 = fig[2, 1] = Axis(fig, xlabel="time, s", ylabel="dipole, Pa") -# ax3 = fig[3, 1] = Axis(fig, xlabel="time, s", ylabel="total, Pa") -# l1 = lines!(ax1, time, apth_total.p_m) -# l2 = lines!(ax2, time, apth_total.p_d) -# l3 = lines!(ax3, time, apth_total.p_m.+apth_total.p_d) -# hidexdecorations!(ax1, grid=false) -# hidexdecorations!(ax2, grid=false) -# save(joinpath(@__DIR__, "openfast-apth_total.png"), fig) -# nothing # hide +time_axis = 1 +apth_total = AcousticAnalogies.combine(apth, obs_time_range, num_obs_times, time_axis) +nothing # hide ``` +With that, we're finally able to plot the acoustic pressure time history: + +```@example first_example +fig = Figure() +ax1 = fig[1, 1] = Axis(fig, xlabel="time, s", ylabel="monopole, Pa") +ax2 = fig[2, 1] = Axis(fig, xlabel="time, s", ylabel="dipole, Pa") +ax3 = fig[3, 1] = Axis(fig, xlabel="time, s", ylabel="total, Pa") +l1 = lines!(ax1, time, apth_total.p_m) +l2 = lines!(ax2, time, apth_total.p_d) +l3 = lines!(ax3, time, apth_total.p_m.+apth_total.p_d) +hidexdecorations!(ax1, grid=false) +hidexdecorations!(ax2, grid=false) +save(joinpath(@__DIR__, "openfast-apth_total.png"), fig) +nothing # hide +``` +![Acoustic Pressure Time History](openfast-apth_total.png) + The plot shows that the monopole/thickness noise is much lower than the dipole/loading noise. Wind turbine blades are relatively slender, which would tend to reduce thickness noise. Also the observer is downstream of the rotation plane, which is where loading noise is traditionally thought to dominate (monopole/thickness noise is more significant in the rotor rotation plane, usually). +(Although we didn't use the actual cross-sectional area for the blades, which directly affects the monopole/thickness noise.) We now calculate the overall sound pressure level from the acoustic pressure time history. Next, we will calculate the narrowband spectrum. -Finally, we will calculate the overall A-weighted sound pressure level from the narrowband spectrum. +Finally, we will calculate the overall sound pressure level from the narrowband spectrum. ```@example first_example -# oaspl_from_apth = AcousticMetrics.OASPL(apth_total) -# nbs = AcousticMetrics.MSPSpectrumAmplitude(apth_total) -# oaspl_from_nbs = AcousticMetrics.OASPL(nbs) -# (oaspl_from_apth, oaspl_from_nbs) -# nothing # hide +oaspl_from_apth = AcousticMetrics.OASPL(apth_total) +nbs = AcousticMetrics.MSPSpectrumAmplitude(apth_total) +oaspl_from_nbs = AcousticMetrics.OASPL(nbs) +@show oaspl_from_apth oaspl_from_nbs +nothing # hide ``` -As a last step, we create VTK files that one can visualize in Paraview (or similar software) - -```@example first_example -# name = joinpath(@__DIR__, "vtk", "iea3p4_vtk") -# mkpath(dirname(name)) -# outfiles = AcousticAnalogies.to_paraview_collection(name, ses) -# end # module -# nothing # hide -``` +The OASPL values calculated from the acoustic pressure time history and the narrowband spectrum are the same, as they should be according to [Parseval's theorem](https://en.wikipedia.org/wiki/Parseval%27s_theorem). diff --git a/src/AcousticAnalogies.jl b/src/AcousticAnalogies.jl index b7f926e7..6031b5f0 100644 --- a/src/AcousticAnalogies.jl +++ b/src/AcousticAnalogies.jl @@ -5,14 +5,15 @@ using CCBlade: CCBlade using CSV: CSV using DataFrames: DataFrames using FillArrays: Fill -using FLOWMath: akima, linear, ksmax, norm_cs_safe, dot_cs_safe, atan_cs_safe, abs_cs_safe +using FLOWMath: FLOWMath, akima, linear, ksmax, norm_cs_safe, dot_cs_safe, atan_cs_safe, abs_cs_safe using FlexiMaps: mapview using Format: format, FormatExpr using JuliennedArrays: JuliennedArrays using KinematicCoordinateTransformations: KinematicTransformation, SteadyRotXTransformation, ConstantVelocityTransformation, compose using LinearAlgebra: cross, norm, mul! using Meshes: Meshes -using StaticArrays: @SVector +using StaticArrays: @SVector, SVector +using Statistics: mean using WriteVTK: WriteVTK include("utils.jl") @@ -63,9 +64,12 @@ export f1a_source_elements_ccblade, tblte_source_elements_ccblade, lblvs_source_ include("bpm_test_utils.jl") include("openfast_helpers.jl") -export read_openfast_file, OpenFASTData +export AbstractTimeDerivMethod, NoTimeDerivMethod, SecondOrderFiniteDiff, calculate_loading_dot! +export AbstractRadialInterpMethod, FLOWLinearInterp, FLOWAkimaInterp, interpolate_to_cell_centers! +export OpenFASTData, read_openfast_file, f1a_source_elements_openfast include("writevtk.jl") +export to_paraview_collection include("deprecated.jl") diff --git a/src/f1a.jl b/src/f1a.jl index 7a407479..4e17f190 100644 --- a/src/f1a.jl +++ b/src/f1a.jl @@ -68,6 +68,51 @@ function CompactF1ASourceElement(ρ0, c0, r, θ, Δr, Λ, fn, fr, fc, τ) return CompactF1ASourceElement(ρ0, c0, Δr, Λ, y0dot, y1dot, y2dot, y3dot, f0dot, f1dot, τ, span_uvec) end +""" + CompactF1ASourceElement(ρ0, c0, r, θ, Δr, Λ, fn, fndot, fr, frdot, fc, fcdot, τ) + +Construct a source element to be used with the compact form of Farassat's formulation 1A, using position and loading data expressed in a cylindrical coordinate system. + +The `r` and `θ` arguments are used to define the radial and circumferential position of the source element in a cylindrical coordinate system. +Likewise, the `fn`, `fr`, and `fc` arguments are used to define the normal, radial, and circumferential loading per unit span *on the fluid* (in a reference frame moving with the element) in the same cylindrical coordinate system. +The `fndot`, `frdot`, and `fcdot` arguments are the time-derivative of the normal, radial, and circumferential loading per unit span, again *on the fluid* and in a reference frame moving with the element, in the cylindrical coordinate system. +The cylindrical coordinate system is defined as follows: + + * The normal axial direction is in the positive x axis + * The circumferential/azimuth angle `θ` is defined such that `θ = 0` means the radial direction is aligned with the positive y axis, and a positive `θ` indicates a right-handed rotation around the positive x axis. + +Note that, for a proper noise prediction, the source element needs to be transformed into the "global" frame, aka, the reference frame of the fluid. +This can be done easily with the transformations provided by the `KinematicCoordinateTransformations` package, or manually by modifying the components of the source element struct. + +# Arguments +- ρ0: Ambient air density (kg/m^3) +- c0: Ambient speed of sound (m/s) +- r: radial coordinate of the element in the blade-fixed coordinate system (m) +- θ: angular offest of the element in the blade-fixed coordinate system (rad) +- Δr: length of the element (m) +- Λ: cross-sectional area of the element (m^2) +- fn: normal load per unit span *on the fluid* (N/m) +- fndot: time derivative of the normal load per unit span *on the fluid* (N/(m*s)) +- fr: radial load *on the fluid* (N/m) +- frdot: time derivative of the radial load *on the fluid* (N/(m*s)) +- fc: circumferential load *on the fluid* (N/m) +- fcdot: time derivative of the circumferential load *on the fluid* (N/(m*s)) +- τ: source time (s) +""" +function CompactF1ASourceElement(ρ0, c0, r, θ, Δr, Λ, fn, fndot, fr, frdot, fc, fcdot, τ) + s, c = sincos(θ) + y0dot = @SVector [0, r*c, r*s] + T = eltype(y0dot) + y1dot = @SVector zeros(T, 3) + y2dot = @SVector zeros(T, 3) + y3dot = @SVector zeros(T, 3) + f0dot = @SVector [fn, c*fr - s*fc, s*fr + c*fc] + f1dot = @SVector [fndot, c*frdot - s*fcdot, s*frdot + c*fcdot] + span_uvec = @SVector [0, c, s] + + return CompactF1ASourceElement(ρ0, c0, Δr, Λ, y0dot, y1dot, y2dot, y3dot, f0dot, f1dot, τ, span_uvec) +end + """ (trans::KinematicTransformation)(se::CompactF1ASourceElement) diff --git a/src/openfast_helpers.jl b/src/openfast_helpers.jl index d1fd7234..f4dc3fe5 100644 --- a/src/openfast_helpers.jl +++ b/src/openfast_helpers.jl @@ -1,91 +1,230 @@ +abstract type AbstractTimeDerivMethod end +struct NoTimeDerivMethod <: AbstractTimeDerivMethod end +struct SecondOrderFiniteDiff <: AbstractTimeDerivMethod end + +abstract type AbstractRadialInterpMethod end +struct FLOWLinearInterp <: AbstractRadialInterpMethod end +struct FLOWAkimaInterp <: AbstractRadialInterpMethod end + """ Struct for holding data from an OpenFAST (AeroDyn?) output file. # Fields: * `time`: vector of simulation times with size `(num_times,)` +* `dtime_dtau`: vector of derivative of simulation times with respect to non-dimensional/computational time with size `(num_times,)` * `v`: vector of freestream velocity time history with size `(num_times,)` * `azimuth`: vector of azimuth angle time history with size `(num_times,)` * `omega`: vector of rotation rate time history with size `(num_times,)` * `pitch`: array of pitch angle time history with size `(num_times, num_blades)` +* `radii`: vector of blade radial locations with size `(num_radial,)` +* `radii_mid`: vector of cell-centered/midpoint blade radial locations with size `(num_radial-1,)` +* `cs_area`: vector of cross-sectional areas with size `(num_radial,)`. +* `cs_area_mid`: vector of cell-centered/midpoint cross-sectional areas with size `(num_radial-1,)`. * `axial_loading`: array of axial loading time history with size `(num_times, num_radial, num_blades)` +* `axial_loading_mid`: array of axial loading time history at cell-centered blade radial locations with size `(num_times, num_radial-1, num_blades)` +* `axial_loading_mid_dot`: array of axial loading temporal derivative time history at cell-centered blade radial locations with size `(num_times, num_radial-1, num_blades)` * `circum_loading`: array of circumferential loading time history with size `(num_times, num_radial, num_blades)` +* `circum_loading_mid`: array of circum loading time history at cell-centered blade radial locations with size `(num_times, num_radial-1, num_blades)` +* `circum_loading_mid_dot`: array of circum loading temporal derivative time history at cell-centered blade radial locations with size `(num_times, num_radial-1, num_blades)` """ -struct OpenFASTData{TTime,TV,TAzimuth,TOmega,TPitch,TAxialLoading,TTangentialLoading} +struct OpenFASTData{TRadialInterpMethod,TTimeDerivMethod, + TTime,TdTimedTau,TV,TAzimuth,TOmega,TPitch, + TRadii,TRadiiMid,TDRadii, + TCSArea,TCSAreaMid, + TAxialLoading,TAxialLoadingMid,TAxialLoadingMidDot, + TCircumLoading,TCircumLoadingMid,TCircumLoadingMidDot} time::TTime + dtime_dtau::TdTimedTau v::TV azimuth::TAzimuth omega::TOmega pitch::TPitch + radii::TRadii + radii_mid::TRadiiMid + dradii::TDRadii + cs_area::TCSArea + cs_area_mid::TCSAreaMid axial_loading::TAxialLoading - circum_loading::TTangentialLoading - - function OpenFASTData(time, v, azimuth, omega, pitch, axial_loading, circum_loading) - # Figure out what num_times is. - if time !== nothing - num_times = length(time) - elseif v !== nothing - num_times = length(v) - elseif azimuth !== nothing - num_times = length(azimuth) - elseif omega !== nothing - num_times = length(omega) - elseif pitch !== nothing - num_times = size(pitch, 1) - elseif axial_loading !== nothing - num_times = size(axial_loading, 1) - elseif circum_loading !== nothing - num_times = size(circum_loading, 1) - end + axial_loading_mid::TAxialLoadingMid + axial_loading_mid_dot::TAxialLoadingMidDot + circum_loading::TCircumLoading + circum_loading_mid::TCircumLoadingMid + circum_loading_mid_dot::TCircumLoadingMidDot + num_blades::Int + + function OpenFASTData{ + TRadialInterpMethod,TTimeDerivMethod, + TTime,TdTimedTau,TV,TAzimuth,TOmega,TPitch, + TRadii,TRadiiMid,TDRadii, + TCSArea,TCSAreaMid, + TAxialLoading,TAxialLoadingMid,TAxialLoadingMidDot, + TCircumLoading,TCircumLoadingMid,TCircumLoadingMidDot}( + time, dtime_dtau, v, azimuth, omega, pitch, + radii, radii_mid, dradii, + cs_area, cs_area_mid, + axial_loading, axial_loading_mid, axial_loading_mid_dot, + circum_loading, circum_loading_mid, circum_loading_mid_dot) where { + TRadialInterpMethod,TTimeDerivMethod, + TTime,TdTimedTau,TV,TAzimuth,TOmega,TPitch, + TRadii,TRadiiMid,TDRadii, + TCSArea,TCSAreaMid, + TAxialLoading,TAxialLoadingMid,TAxialLoadingMidDot, + TCircumLoading,TCircumLoadingMid,TCircumLoadingMidDot} + + num_times = length(time) + num_radial = length(radii) + num_radial_mid = num_radial - 1 # Figure out what num_blades is. if pitch !== nothing num_blades = size(pitch, 2) elseif axial_loading !== nothing num_blades = size(axial_loading, 3) + elseif axial_loading_mid !== nothing + num_blades = size(axial_loading_mid, 3) + elseif axial_loading_mid_dot !== nothing + num_blades = size(axial_loading_mid_dot, 3) elseif circum_loading !== nothing num_blades = size(circum_loading, 3) + elseif circum_loading_mid !== nothing + num_blades = size(circum_loading_mid, 3) + elseif circum_loading_mid_dot !== nothing + num_blades = size(circum_loading_mid_dot, 3) + else + num_blades = 0 end - # Figure out what num_radial is. - if axial_loading !== nothing - num_radial = size(axial_loading, 2) - elseif circum_loading !== nothing - num_radial = size(circum_loading, 2) - end - - if time !== nothing - # I don't think this can ever happen since the check for time !== nothing at the beginning is first, but whatever. - size(time) == (num_times,) || raise(ArgumentError("size(time) = $(size(time)) inconsistent with other inputs")) + if dtime_dtau !== nothing + size(dtime_dtau) == (num_times,) || throw(ArgumentError("size(dtime_dtau) = $(size(dtime_dtau)) does not match size(time) = $(size(time))")) end if v !== nothing - size(v) == (num_times,) || raise(ArgumentError("size(v) = $(size(v)) does not match size(time) = $(size(time))")) + size(v) == (num_times,) || throw(ArgumentError("size(v) = $(size(v)) does not match size(time) = $(size(time))")) end if azimuth !== nothing - size(azimuth) == (num_times,) || raise(ArgumentError("size(azimuth) = $(size(azimuth)) does not match size(time) = $(size(time))")) + size(azimuth) == (num_times,) || throw(ArgumentError("size(azimuth) = $(size(azimuth)) does not match size(time) = $(size(time))")) end if omega !== nothing - size(omega) == (num_times,) || raise(ArgumentError("size(omega) = $(size(omega)) does not match size(time) = $(size(time))")) + size(omega) == (num_times,) || throw(ArgumentError("size(omega) = $(size(omega)) does not match size(time) = $(size(time))")) end if pitch !== nothing - size(pitch) == (num_times, num_blades) || raise(ArgumentError("size(pitch) = $(size(pitch)) not consistent with size(time) = $(size(time)) and/or size(axial_loading) = $(size(axial_loading))")) + size(pitch) == (num_times, num_blades) || throw(ArgumentError("size(pitch) = $(size(pitch)) not consistent with size(time) = $(size(time)) and/or size(axial_loading) = $(size(axial_loading))")) + end + + if radii_mid !== nothing + size(radii_mid) == (num_radial_mid,) || throw(ArgumentError("size(radii_mid) = $(size(radii_mid)) not consistent with length(radii)-1 = $(num_radial_mid)")) + end + + if dradii !== nothing + size(dradii) == (num_radial_mid,) || throw(ArgumentError("size(dradii) = $(size(dradii)) not consistent with length(radii)-1 = $(num_radial_mid)")) + end + + if cs_area !== nothing + size(cs_area) == (num_radial,) || throw(ArgumentError("size(cs_area) = $(size(cs_area)) not consistent with length(radii) = $(num_radial)")) + end + + if cs_area_mid !== nothing + size(cs_area_mid) == (num_radial_mid,) || throw(ArgumentError("size(cs_area_mid) = $(size(cs_area_mid)) not consistent with length(radii)-1 = $(num_radial_mid)")) end if axial_loading !== nothing - size(axial_loading) == (num_times, num_radial, num_blades) || raise(ArgumentError("size(axial_loading) = $(size(axial_loading)) not consistent with size(time) = $(size(time))")) + size(axial_loading) == (num_times, num_radial, num_blades) || throw(ArgumentError("size(axial_loading) = $(size(axial_loading)) not consistent with size(time) = $(size(time)), length(radii) = $(num_radial), blade count = $(num_blades)")) + end + + if axial_loading_mid !== nothing + size(axial_loading_mid) == (num_times, num_radial_mid, num_blades) || throw(ArgumentError("size(axial_loading_mid) = $(size(axial_loading_mid)) not consistent with size(time) = $(size(time)), length(radii)-1 = $(num_radial_mid), blade count = $(num_blades)")) + end + + if axial_loading_mid_dot !== nothing + size(axial_loading_mid_dot) == (num_times, num_radial_mid, num_blades) || throw(ArgumentError("size(axial_loading_mid_dot) = $(size(axial_loading_mid_dot)) not consistent with size(time) = $(size(time)), length(radii)-1 = $(num_radial_mid), blade count = $(num_blades)")) end if circum_loading !== nothing - size(circum_loading) == (num_times, num_radial, num_blades) || raise(ArgumentError("size(circum_loading) = $(size(circum_loading)) not consistent with size(time) = $(size(time)) and/or size(axial_loading) = $(size(axial_loading))")) + size(circum_loading) == (num_times, num_radial, num_blades) || throw(ArgumentError("size(circum_loading) = $(size(circum_loading)) not consistent with size(time) = $(size(time)), length(radii) = $(num_radial), blade count = $(num_blades)")) + end + + if circum_loading_mid !== nothing + size(circum_loading_mid) == (num_times, num_radial_mid, num_blades) || throw(ArgumentError("size(circum_loading_mid) = $(size(circum_loading_mid)) not consistent with size(time) = $(size(time)), length(radii)-1 = $(num_radial_mid), blade count = $(num_blades)")) + end + + if circum_loading_mid_dot !== nothing + size(circum_loading_mid_dot) == (num_times, num_radial_mid, num_blades) || throw(ArgumentError("size(circum_loading_mid_dot) = $(size(circum_loading_mid_dot)) not consistent with size(time) = $(size(time)), length(radii)-1 = $(num_radial_mid), blade count = $(num_blades)")) end - return new{typeof(time),typeof(v),typeof(azimuth),typeof(omega),typeof(pitch),typeof(axial_loading),typeof(circum_loading)}(time, v, azimuth, omega, pitch, axial_loading, circum_loading) + return new{TRadialInterpMethod,TTimeDerivMethod, + typeof(time),typeof(dtime_dtau),typeof(v),typeof(azimuth),typeof(omega),typeof(pitch), + typeof(radii),typeof(radii_mid),typeof(dradii), + typeof(cs_area),typeof(cs_area_mid), + typeof(axial_loading),typeof(axial_loading_mid),typeof(axial_loading_mid_dot), + typeof(circum_loading),typeof(circum_loading_mid),typeof(circum_loading_mid_dot)}( + time, dtime_dtau, v, azimuth, omega, pitch, + radii, radii_mid, dradii, + cs_area, cs_area_mid, + axial_loading, axial_loading_mid, axial_loading_mid_dot, + circum_loading, circum_loading_mid, circum_loading_mid_dot, + num_blades) end end +function OpenFASTData{TRadialInterpMethod,TTimeDerivMethod}(time, v, azimuth, omega, pitch, radii, cs_area, axial_loading, circum_loading) where {TRadialInterpMethod<:AbstractRadialInterpMethod,TTimeDerivMethod<:AbstractTimeDerivMethod} + dtime_dtau = similar(time) + + # Find the blade element midpoint locations. + radii_mid = 0.5 .* (@view(radii[begin:end-1]) .+ @view(radii[begin+1:end])) + + # Find the radial spacing. + dradii = @view(radii[begin+1:end]) .- @view(radii[begin:end-1]) + + if cs_area !== nothing + cs_area_mid = similar(cs_area, length(cs_area)-1) + else + cs_area_mid = nothing + end + + if axial_loading !== nothing + axial_loading_mid = similar(axial_loading, size(axial_loading, 1), size(axial_loading, 2)-1, size(axial_loading, 3)) + axial_loading_mid_dot = zero(axial_loading_mid) + else + axial_loading_mid = axial_loading_mid_dot = nothing + end + + if circum_loading !== nothing + circum_loading_mid = similar(circum_loading, size(circum_loading, 1), size(circum_loading, 2)-1, size(circum_loading, 3)) + circum_loading_mid_dot = zero(circum_loading_mid) + else + circum_loading_mid = circum_loading_mid_dot = nothing + end + + return OpenFASTData{TRadialInterpMethod,TTimeDerivMethod}( + time, dtime_dtau, v, azimuth, omega, pitch, + radii, radii_mid, dradii, + cs_area, cs_area_mid, + axial_loading, axial_loading_mid, axial_loading_mid_dot, + circum_loading, circum_loading_mid, circum_loading_mid_dot) +end + +function OpenFASTData{TRadialInterpMethod,TTimeDerivMethod}( + time, dtime_dtau, v, azimuth, omega, pitch, + radii, radii_mid, dradii, + cs_area, cs_area_mid, + axial_loading, axial_loading_mid, axial_loading_mid_dot, + circum_loading, circum_loading_mid, circum_loading_mid_dot) where {TRadialInterpMethod,TTimeDerivMethod} + return OpenFASTData{TRadialInterpMethod,TTimeDerivMethod, + typeof(time),typeof(dtime_dtau),typeof(v),typeof(azimuth),typeof(omega),typeof(pitch), + typeof(radii),typeof(radii_mid),typeof(dradii), + typeof(cs_area),typeof(cs_area_mid), + typeof(axial_loading),typeof(axial_loading_mid),typeof(axial_loading_mid_dot), + typeof(circum_loading),typeof(circum_loading_mid),typeof(circum_loading_mid_dot)}( + time, dtime_dtau, v, azimuth, omega, pitch, + radii, radii_mid, dradii, + cs_area, cs_area_mid, + axial_loading, axial_loading_mid, axial_loading_mid_dot, + circum_loading, circum_loading_mid, circum_loading_mid_dot) +end + function _get_num_blades(fmt, column_names) # Match the format against the column names. ms = match.(fmt, column_names) @@ -130,8 +269,113 @@ function _get_num_blades_num_radial(fmt, column_names) return num_blades, num_radial, ms_only_matches end +function interpolate_to_cell_centers!(data::OpenFASTData{FLOWLinearInterp}) + if (data.cs_area !== nothing) && (data.cs_area_mid !== nothing) + data.cs_area_mid .= FLOWMath.linear.(Ref(data.radii), Ref(data.cs_area), data.radii_mid) + end + + if (data.axial_loading !== nothing) && (data.axial_loading_mid !== nothing) + for bidx in 1:size(data.axial_loading_mid, 3) + for tidx in 1:size(data.axial_loading_mid, 1) + data.axial_loading_mid[tidx, :, bidx] .= FLOWMath.linear.(Ref(data.radii), Ref(@view(data.axial_loading[tidx, :, bidx])), data.radii_mid) + end + end + end + + if (data.circum_loading !== nothing) && (data.circum_loading_mid !== nothing) + for bidx in 1:size(data.circum_loading_mid, 3) + for tidx in 1:size(data.circum_loading_mid, 1) + data.circum_loading_mid[tidx, :, bidx] .= FLOWMath.linear.(Ref(data.radii), Ref(@view(data.circum_loading[tidx, :, bidx])), data.radii_mid) + end + end + end + + return nothing +end + +function interpolate_to_cell_centers!(data::OpenFASTData{FLOWAkimaInterp}) + if (data.cs_area !== nothing) && (data.cs_area_mid !== nothing) + spline_cs_area = FLOWMath.Akima(data.radii, data.cs_area) + data.cs_area_mid .= spline_cs_area.(data.radii_mid) + end + + if (data.axial_loading !== nothing) && (data.axial_loading_mid !== nothing) + for bidx in 1:size(data.axial_loading_mid, 3) + for tidx in 1:size(data.axial_loading_mid, 1) + spline_axial = FLOWMath.Akima(data.radii, @view(data.axial_loading[tidx, :, bidx])) + data.axial_loading_mid[tidx, :, bidx] .= spline_axial.(data.radii_mid) + end + end + end + + if (data.circum_loading !== nothing) && (data.circum_loading_mid !== nothing) + for bidx in 1:size(data.circum_loading_mid, 3) + for tidx in 1:size(data.circum_loading_mid, 1) + spline_circum = FLOWMath.Akima(data.radii, @view(data.circum_loading[tidx, :, bidx])) + data.circum_loading_mid[tidx, :, bidx] .= spline_circum.(data.radii_mid) + end + end + end + + return nothing +end + +function calculate_loading_dot!(data::OpenFASTData{TRadialInterpMethod,NoTimeDerivMethod}) where {TRadialInterpMethod} + if data.axial_loading_mid_dot !== nothing + fill!(data.axial_loading_mid_dot, 0) + end + + if data.circum_loading_mid_dot !== nothing + fill!(data.circum_loading_mid_dot, 0) + end + + return nothing +end + +function _finite_diff_2nd_order!(df_dtau, f) + # `f` is the input array, which we assume is of size `(num_times, num_radial, num_blades)`. + # `df_dtau` is the derivative wrt `tau`, the non-dimensional time. + + @views begin + # These stencils are in Tannehill, Anderson, Pletcher, "Computational Fluid Mechanics and Heat Transfer," 2nd edition, page 50. + # First do the interior points. + df_dtau[begin+1:end-1, :, :] .= 0.5 .* (f[begin+2:end, :, :] .- f[begin:end-2, :, :]) + + # Then the left boundary. + df_dtau[begin, :, :] .= 0.5 .* (-3 .* f[begin, :, :] .+ 4 .* f[begin+1, :, :] .- f[begin+2, :, :]) + + # Then the right boundary. + df_dtau[end, :, :] .= 0.5 .* (3 .* f[end, :, :] .- 4 .* f[end-1, :, :] .+ f[end-2, :, :]) + end + + return nothing +end + +function calculate_loading_dot!(data::OpenFASTData{TRadialInterpMethod,SecondOrderFiniteDiff}) where {TRadialInterpMethod} + # First get dt/dτ. + _finite_diff_2nd_order!(data.dtime_dtau, data.time) + + if (data.axial_loading_mid !== nothing) && (data.axial_loading_mid_dot !== nothing) + # Now get the derivatitve of the axial loading wrt tau, the non-dimensional time. + _finite_diff_2nd_order!(data.axial_loading_mid_dot, data.axial_loading_mid) + + # Now get the derivative of the axial loading with respect to the dimensional time via `dfdt = df/dτ*dτ/dt = (df/dτ)/(dt/dτ) + data.axial_loading_mid_dot ./= data.dtime_dtau + end + + if (data.circum_loading_mid !== nothing) && (data.circum_loading_mid_dot !== nothing) + # Now get the derivatitve of the circum loading wrt tau, the non-dimensional time. + _finite_diff_2nd_order!(data.circum_loading_mid_dot, data.circum_loading_mid) + + # Now get the derivative of the circum loading with respect to the dimensional time via `dfdt = df/dτ*dτ/dt = (df/dτ)/(dt/dτ) + data.circum_loading_mid_dot ./= data.dtime_dtau + end + + return nothing +end + """ - read_openfast_file(fname; + read_openfast_file(fname, radii, cs_area=nothing; header_keyword="Time", has_units_header=true, time_column_name="Time", @@ -140,33 +384,49 @@ end omega_column_name="RotSpeed", pitch_fmt=r"BlPitch(?[[:digit:]]+)", axial_loading_fmt=r"AB(?[[:digit:]]+)N(?[[:digit:]]+)Fxl", - circum_loading_fmt=r"AB(?[[:digit:]]+)N(?[[:digit:]]+)Fyl") + circum_loading_fmt=r"AB(?[[:digit:]]+)N(?[[:digit:]]+)Fyl", + radial_interp_method=FLOWLinearInterp, + time_deriv_method=SecondOrderFiniteDiff) Read an OpenFAST output file and return a [`OpenFASTData`](@ref) object. +The `Azimuth` and `BlPitch` columns are assumed to be in degrees and will be converted to radians. +Likewise, the `RotSpeed` column is assumed to be in revolutions per minute and will be converted to radians per second. + # Arguments * `fname`: name of the OpenFAST output file to read +* `radii`: `Vector` of blade radial coordinates +* `cs_area`: `Vector` of radial distribution of cross-sectional areas, or `nothing` to ignore * `header_keyword="Time"`: string at the beginning of the header line (maybe always "Time"?) * `has_units_header=true`: if true, assume the file has a line directly after the header line with the units of each column * `time_column_name=header_keyword`: name of time column in file. Set to `nothing` to skip. * `freestream_vel_column_name`: name of the freestream velocity column in the file. Set to `nothing` to skip. * `azimuth_column_name`: name of the azimuth column in the file. Set to `nothing` to skip. -* `omega_column_name`: name of the omega column in the file. Set to `nothing` to skip. +* `omega_column_name`: name of the omega (rotation rate) Set to `nothing` to skip. * `pitch_fmt`: Format for finding all pitch columns in the file. Should be a regex with a capture group named `blade` for the blade index, or `nothing` to skip. * `axial_loading_fmt`: Format for finding all axial loading columns in the file. Should be a regex with a captures groups named `blade` and `radial` for the blade and radial indices, or `nothing` to skip. * `circum_loading_fmt`: Format for finding all radial loading columns in the file. Should be a regex with a captures groups named `blade` and `radial` for the blade and radial indices, or `nothing` to skip. +* `radial_interp_method`: `<:AbstractRadialInterpMethod` indicating method used to interpolate loading from blade element "interfaces" to midpoints. +* `time_deriv_method`: `<:AbstractTimeDerivMethod` indicating the method used to calculate the loading time derivatives. +* `average_freestream_vel=false`: Store possibily unsteady freestream velocity in the `OpenFASTData` object if `false`, store average value otherwise. +* `average_omega=false`: Store possibily unsteady omega (rotation rate) in the `OpenFASTData` object if `false`, store average value otherwise. """ -function read_openfast_file(fname; +function read_openfast_file(fname, radii, cs_area=nothing; header_keyword="Time", has_units_header=true, - # time_column_name="Time", time_column_name=header_keyword, freestream_vel_column_name="Wind1VelX", azimuth_column_name="Azimuth", omega_column_name="RotSpeed", pitch_fmt=r"BlPitch(?[[:digit:]]+)", axial_loading_fmt=r"AB(?[[:digit:]]+)N(?[[:digit:]]+)Fxl", - circum_loading_fmt=r"AB(?[[:digit:]]+)N(?[[:digit:]]+)Fyl") + circum_loading_fmt=r"AB(?[[:digit:]]+)N(?[[:digit:]]+)Fyl", + radial_interp_method=FLOWLinearInterp, + time_deriv_method=SecondOrderFiniteDiff, + average_freestream_vel=true, + average_omega=true) + + num_radial = length(radii) # Remove leading and trailing whitespace from header keyword. header_keyword_s = strip(header_keyword) @@ -214,19 +474,29 @@ function read_openfast_file(fname; if freestream_vel_column_name === nothing v = nothing else - v = df[!, freestream_vel_column_name] + v_tmp = df[!, freestream_vel_column_name] + if average_freestream_vel + v = Fill(mean(v_tmp), length(v_tmp)) + else + v = v_tmp + end end if azimuth_column_name === nothing azimuth = nothing else - azimuth = df[!, azimuth_column_name] + azimuth = df[!, azimuth_column_name] .* (pi/180) end if omega_column_name === nothing omega = nothing else - omega = df[!, omega_column_name] + omega_tmp = df[!, omega_column_name] .* (2*pi/60) + if average_omega + omega = Fill(mean(omega_tmp), length(omega_tmp)) + else + omega = omega_tmp + end end if pitch_fmt === nothing @@ -240,7 +510,7 @@ function read_openfast_file(fname; pitch = Array{TF_pitch, 2}(undef, num_times, pitch_num_blades) for m in pitch_matches b = parse(Int, m[:blade]) - pitch[:, b] .= df[!, m.match] + pitch[:, b] .= df[!, m.match] .* (pi/180) end end @@ -263,10 +533,10 @@ function read_openfast_file(fname; if circum_loading_fmt === nothing circum_loading = nothing else - # Make sure we get the same thing if we use the circumferential loading format. + # Get the number of blades and radial stations according to the circumferential loading format. circum_num_blades, circum_num_radial, circum_matches = _get_num_blades_num_radial(circum_loading_fmt, colnames) - # Decide on an element type for the axial loading. + # Decide on an element type for the circumferential loading. TF_circum = promote_type(eltype.(getproperty.(Ref(df), getproperty.(circum_matches, :match)))...) circum_loading = Array{TF_circum, 3}(undef, num_times, circum_num_radial, circum_num_blades) for m in circum_matches @@ -276,5 +546,109 @@ function read_openfast_file(fname; end end - return OpenFASTData(time, v, azimuth, omega, pitch, axial_loading, circum_loading) + # Create the openfast data struct. + data = OpenFASTData{radial_interp_method,time_deriv_method}(time, v, azimuth, omega, pitch, radii, cs_area, axial_loading, circum_loading) + + # Interpolate the loading to the cell centers. + interpolate_to_cell_centers!(data) + + # Calculate the loading time derivatives. + calculate_loading_dot!(data) + + return data +end + +""" + f1a_source_elements_openfast(data::OpenFASTData, rho0, c0, area_per_chord2::Vector, positive_x_rotation::Bool=true) + +Construct and return an array of `CompactF1ASourceElement` objects from OpenFAST data. + +# Arguments +- `data`: OpenFAST data object. +- `rho0`: Ambient air density (kg/m^3) +- `c0`: Ambient speed of sound (m/s) +- `positive_x_rotation`: rotate blade around the positive-x axis if `true`, negative-x axis otherwise +""" +function f1a_source_elements_openfast(data::OpenFASTData, rho0, c0, positive_x_rotation::Bool=true) + + # if length(area_per_chord2) != length(data.radii_mid) + # throw(ArgumentError("length of area_per_chord2 = $(length(area_per_chord2)) should be equal to length(data.radii_mid) = $(length(data.radii_mid))")) + # end + + # OK, what's the coordinate system here? + # Well, I guess I can decide. + # For CCBlade, I've been assuming that the blade elements are translating in the positive x direction, rotating about the positive x axis if `positive_x_rotation` is `true`, negative x axis otherwise. + # So for consistency let's say I do the same here. + # Then that would mean the freestream velocity would be pointed in the negative x direction. + # The velocity in the example OpenFAST file is always positive, so, I'll need to keep that in mind. + # For the position of the blade elements, we'll initially be aligned with the y axis, and I'm assuming that the radial locations are all positive, so no sign switch necessary there. + + # But what to do about the loading? + # What is the loading sign convention? + # Well, for the axial loading, I would expect that the axial loading on the fluid would oppose the freestream velocity for a wind turbine. + # So, if the blades are translating in the positive x direction, the axial loading on the fluid would also be in the positive x direction. + # I can see that the axial loading in the openfast file is all positive, so no sign switch is necessary. + # For the circumferential loading, for a wind turbine, I would expect the loading on the blade to be in the same direction as the blade rotation, so the loading on the fluid would be opposite the blade's rotation. + # So, if the first blade is initially aligned with the y axis, rotating around the positive x axis, then the blade would initially be moving in the positive z direction, and so the loading on the fluid should be in the negative z direction. + # And it looks like the circumferential loading in the OpenFAST file is negative, so no sign switch necessary. + # But if the blade is rotating about the negative x axis, then the loading on the blade would be in the negative direction, and thus the loading on the fluid would be in the positive direction, and so we'd need to switch the sign. + + # Now, let's get some transformation stuff going. + # We're going to be translating in the positive x direction. + # Oh, but what do we do about the fact that the axial velocity isn't necessarily constant? + # Well, what I need is the position and velocity in the axial direction. + # It would also be nice to get the acceleration and jerk, too. + # OK, I have the velocity, but the position is the tricky part. + # But I should be able to just integrate it, I suppose. + # So let's do that. + # We'll assume the position at the first time is at the origin. + # Now we need to integrate the velocity. + # (This assumes the hub starts at the origin. + # If that's not the case, just add `x0` to it, where `x0` is the axial position at the first time level.) + x = FLOWMath.cumtrapz(data.time, data.v) + + # So now we have the axial position and axial velocity as a function of data.time. + # Now we should be able to create transformations that take that into account. + # There will be one for each data.time level. + # For each entry in `data.time`, i.e. for `data.time[i]`, it will start at `[x[i], 0, 0]` and have velocity `[v[i], 0, 0]`. + # This won't take into account the effect the possibily non-constant velocity has on the acceleration or jerk. + # const_vel_trans = ConstantVelocityTransformation.(data.time, SVector{3,typeof(x)}.(x, 0, 0), SVector{3,typeof(x)}.(data.v, 0, 0)) + + # Next, we need to figure out the rotation stuff. + # We actually have both a time history of omega and azimuthal angles. + # So we could create a bunch of rotational transformations that have the correct azimuth offset and omega. + # But it wouldn't take into account the effect the time derivative of omega has on the things we care about: the derivatives of position and loading. + # I would need to do some work on that. + # Also, I'm going to assume that omega and azimuth are always positive, and so I'll switch their signs if we are rotating about the negative x axis. + # rot_trans = SteadyRotXTransformation.(data.time, data.omega.*ifelse(positive_x_rotation, 1, -1), data.azimuth.*ifelse(positive_x_rotation, 1, -1)) + + # So, want everything to have shape (num_times, num_radial_mid, num_blades). + radii_mid_rs = reshape(data.radii_mid, 1, :) + dradii_rs = reshape(data.dradii, 1, :) + + num_blades = data.num_blades + thetas = 2*pi/num_blades.*(0:(num_blades-1)) .* ifelse(positive_x_rotation, 1, -1) + thetas_rs = reshape(thetas, 1, 1, :) + # cs_area_rs = reshape(cs_area, 1, :) + cs_area_rs = reshape(data.cs_area_mid, 1, :) + + # All the loading arrays are in the right shape, and no sign switch is necessary. + fn = data.axial_loading_mid + fndot = data.axial_loading_mid_dot + # The OpenFAST file doesn't have any data for the loading in the radial direction (which should be quite small of course). + fr = zero(eltype(fn)) + frdot = zero(eltype(fn)) + fc = data.circum_loading_mid + fcdot = data.circum_loading_mid_dot + + # Now we should be able to do all this. + # This is a bit too cute, but I can't help myself. + ses = (CompactF1ASourceElement.(rho0, c0, radii_mid_rs, thetas_rs, dradii_rs, cs_area_rs, fn, fndot, fr, frdot, fc.*ifelse(positive_x_rotation, 1, -1), fcdot.*ifelse(positive_x_rotation, 1, -1), data.time) .|> + compose.(data.time, + ConstantVelocityTransformation.(data.time, SVector{3,eltype(x)}.(x, 0, 0), SVector{3,eltype(data.v)}.(data.v, 0, 0)), + SteadyRotXTransformation.(data.time, data.omega.*ifelse(positive_x_rotation, 1, -1), data.azimuth.*ifelse(positive_x_rotation, 1, -1)) + ) + ) + + return ses end diff --git a/test/openfast_helper_tests.jl b/test/openfast_helper_tests.jl index fee261fc..44974d5c 100644 --- a/test/openfast_helper_tests.jl +++ b/test/openfast_helper_tests.jl @@ -1,81 +1,47 @@ module OpenFASTHelperTests -using AcousticAnalogies -using Statistics: mean -using Test +using SafeTestsets: @safetestset +using Test: @testset -@testset "OpenFAST reader test" begin - @testset "default" begin - fname = joinpath(@__DIR__, "gen_test_data", "openfast_data", "IEA-3.4-130-RWT-small.out") +@testset "OpenFAST reader tests" begin - data = read_openfast_file(fname) - num_times = length(data.time) - num_blades = size(data.pitch, 2) - num_radial = size(data.axial_loading, 2) - @test size(data.time) == (num_times,) - @test size(data.v) == (num_times,) - @test size(data.azimuth) == (num_times,) - @test size(data.omega) == (num_times,) - @test size(data.pitch) == (num_times, num_blades) - @test size(data.axial_loading) == (num_times, num_radial, num_blades) - @test size(data.circum_loading) == (num_times, num_radial, num_blades) + @testset "reading" begin + @safetestset "default" begin + include("openfast_reading_default.jl") + end - @test all(data.time .≈ (60.0:0.01:60.1)) - @test all(data.v .≈ 7) - @test all(data.pitch .≈ 0) - # Just some coarse tests. - @test mean(data.omega) ≈ 8.126818181818182 - @test mean(data.axial_loading) ≈ 1632.1274949494953 - @test mean(data.circum_loading) ≈ -217.45748949494939 - end - - @testset "different time column name" begin - fname = joinpath(@__DIR__, "gen_test_data", "openfast_data", "IEA-3.4-130-RWT-small-FooTime.out") + @safetestset "different time column name" begin + include("openfast_reading_different_time_column_name.jl") + end - data = read_openfast_file(fname; header_keyword="FooTime") - num_times = length(data.time) - num_blades = size(data.pitch, 2) - num_radial = size(data.axial_loading, 2) - @test size(data.time) == (num_times,) - @test size(data.v) == (num_times,) - @test size(data.azimuth) == (num_times,) - @test size(data.omega) == (num_times,) - @test size(data.pitch) == (num_times, num_blades) - @test size(data.axial_loading) == (num_times, num_radial, num_blades) - @test size(data.circum_loading) == (num_times, num_radial, num_blades) - - @test all(data.time .≈ (60.0:0.01:60.1)) - @test all(data.v .≈ 7) - @test all(data.pitch .≈ 0) - # Just some coarse tests. - @test mean(data.omega) ≈ 8.126818181818182 - @test mean(data.axial_loading) ≈ 1632.1274949494953 - @test mean(data.circum_loading) ≈ -217.45748949494939 + @safetestset "no units header" begin + include("openfast_reading_no_units_header.jl") + end end - @testset "nu units header" begin - fname = joinpath(@__DIR__, "gen_test_data", "openfast_data", "IEA-3.4-130-RWT-small-no_units.out") + @testset "radial interpolation" begin + @safetestset "linear" begin + include("openfast_radial_interpolation_linear.jl") + end - data = read_openfast_file(fname; has_units_header=false) - num_times = length(data.time) - num_blades = size(data.pitch, 2) - num_radial = size(data.axial_loading, 2) - @test size(data.time) == (num_times,) - @test size(data.v) == (num_times,) - @test size(data.azimuth) == (num_times,) - @test size(data.omega) == (num_times,) - @test size(data.pitch) == (num_times, num_blades) - @test size(data.axial_loading) == (num_times, num_radial, num_blades) - @test size(data.circum_loading) == (num_times, num_radial, num_blades) + @safetestset "akima" begin + include("openfast_radial_interpolation_akima.jl") + end + end - @test all(data.time .≈ (60.0:0.01:60.1)) - @test all(data.v .≈ 7) - @test all(data.pitch .≈ 0) - # Just some coarse tests. - @test mean(data.omega) ≈ 8.126818181818182 - @test mean(data.axial_loading) ≈ 1632.1274949494953 - @test mean(data.circum_loading) ≈ -217.45748949494939 + @testset "time derivatives" begin + @safetestset "constant time step" begin + include("openfast_time_derivatives_constant_time_step.jl") + end + + @safetestset "non-constant time step" begin + include("openfast_time_derivatives_nonconstant_time_step.jl") + end + + @safetestset "no time diff" begin + include("openfast_time_derivatives_nonconstant_time_step_no_diff.jl") + end end end diff --git a/test/openfast_radial_interpolation_akima.jl b/test/openfast_radial_interpolation_akima.jl new file mode 100644 index 00000000..5c606d68 --- /dev/null +++ b/test/openfast_radial_interpolation_akima.jl @@ -0,0 +1,62 @@ +using AcousticAnalogies +using Test + +T = 2.5 +t0 = 0.3 +R = 2.3 +num_times = 64 +num_radial = 30 +num_blades = 3 + +cs(r) = 0.1*sin(2*pi*r/R + 0.1*pi) + 0.2*cos(4*pi*r/R + 0.2*pi) +fn(t, r, b) = 0.2*sin(2*pi/T*t*r/R*b + 0.1*pi) + 0.3*cos(4*pi/T*t*r/R*b + 0.2*pi) +fc(t, r, b) = 0.4*sin(2*pi/T*t*r/R*b + 0.3*pi) + 0.5*cos(4*pi/T*t*r/R*b + 0.4*pi) + +dt = T/num_times +time = t0 .+ (0:(num_times-1)) .* dt + +radii = range(0.2*R, R; length=num_radial) +radii_mid = 0.5 .* (@view(radii[1:end-1]) .+ @view(radii[2:end])) +dradii = nothing + +dtime_dtau = v = azimuth = omega = pitch = nothing +axial_loading_mid_dot = circum_loading_mid_dot = nothing + +r = reshape(radii, 1, :) +b = reshape(1:num_blades, 1, 1, :) +cs_area = @. cs(radii) +axial_loading = @. fn(time, r, b) +circum_loading = @. fc(time, r, b) +cs_area_mid = similar(cs_area, num_radial-1) +axial_loading_mid = similar(axial_loading, num_times, num_radial-1, num_blades) +circum_loading_mid = similar(circum_loading, num_times, num_radial-1, num_blades) + +data = OpenFASTData{FLOWAkimaInterp,SecondOrderFiniteDiff}( + time, dtime_dtau, v, azimuth, omega, pitch, + radii, radii_mid, dradii, + cs_area, cs_area_mid, + axial_loading, axial_loading_mid, axial_loading_mid_dot, + circum_loading, circum_loading_mid, circum_loading_mid_dot) + +interpolate_to_cell_centers!(data) + +r_mid = reshape(radii_mid, 1, :) +cs_area_mid_exact = @. cs(radii_mid) +axial_loading_mid_exact = @. fn(time, r_mid, b) +circum_loading_mid_exact = @. fc(time, r_mid, b) + +csm_min, csm_max = extrema(cs_area_mid_exact) +csm_err = maximum(abs.((data.cs_area_mid .- cs_area_mid_exact)./(csm_max - csm_min))) +@test csm_err < 0.00105 + +alm_min = reshape(minimum(axial_loading_mid_exact; dims=2), num_times, 1, num_blades) +alm_max = reshape(maximum(axial_loading_mid_exact; dims=2), num_times, 1, num_blades) +# @show maximum(abs.((data.axial_loading_mid .- axial_loading_mid_exact)./(alm_max .- alm_min))) +alm_err = maximum(abs.((data.axial_loading_mid .- axial_loading_mid_exact)./(alm_max .- alm_min))) +@test alm_err < 0.033 + +clm_min = reshape(minimum(circum_loading_mid_exact; dims=2), num_times, 1, num_blades) +clm_max = reshape(maximum(circum_loading_mid_exact; dims=2), num_times, 1, num_blades) +# @show maximum(abs.((data.circum_loading_mid .- circum_loading_mid_exact)./(clm_max .- clm_min))) +clm_err = maximum(abs.((data.circum_loading_mid .- circum_loading_mid_exact)./(clm_max .- clm_min))) +@test clm_err < 0.033 diff --git a/test/openfast_radial_interpolation_linear.jl b/test/openfast_radial_interpolation_linear.jl new file mode 100644 index 00000000..fab25e4c --- /dev/null +++ b/test/openfast_radial_interpolation_linear.jl @@ -0,0 +1,60 @@ +using AcousticAnalogies +using Test + +T = 2.5 +t0 = 0.3 +R = 2.3 +num_times = 64 +num_radial = 30 +num_blades = 3 + +cs(r) = 0.1*sin(2*pi*r/R + 0.1*pi) + 0.2*cos(4*pi*r/R + 0.2*pi) +fn(t, r, b) = 0.2*sin(2*pi/T*t*r/R*b + 0.1*pi) + 0.3*cos(4*pi/T*t*r/R*b + 0.2*pi) +fc(t, r, b) = 0.4*sin(2*pi/T*t*r/R*b + 0.3*pi) + 0.5*cos(4*pi/T*t*r/R*b + 0.4*pi) + +dt = T/num_times +time = t0 .+ (0:(num_times-1)) .* dt + +radii = range(0.2*R, R; length=num_radial) +radii_mid = 0.5 .* (@view(radii[1:end-1]) .+ @view(radii[2:end])) +dradii = nothing + +dtime_dtau = v = azimuth = omega = pitch = nothing +axial_loading_mid_dot = circum_loading_mid_dot = nothing + +r = reshape(radii, 1, :) +b = reshape(1:num_blades, 1, 1, :) +cs_area = @. cs(radii) +axial_loading = @. fn(time, r, b) +circum_loading = @. fc(time, r, b) +cs_area_mid = similar(cs_area, num_radial-1) +axial_loading_mid = similar(axial_loading, num_times, num_radial-1, num_blades) +circum_loading_mid = similar(circum_loading, num_times, num_radial-1, num_blades) + +data = OpenFASTData{FLOWLinearInterp,SecondOrderFiniteDiff}( + time, dtime_dtau, v, azimuth, omega, pitch, + radii, radii_mid, dradii, + cs_area, cs_area_mid, + axial_loading, axial_loading_mid, axial_loading_mid_dot, + circum_loading, circum_loading_mid, circum_loading_mid_dot) + +interpolate_to_cell_centers!(data) + +r_mid = reshape(radii_mid, 1, :) +cs_area_mid_exact = @. cs(radii_mid) +axial_loading_mid_exact = @. fn(time, r_mid, b) +circum_loading_mid_exact = @. fc(time, r_mid, b) + +csm_min, csm_max = extrema(cs_area_mid_exact) +csm_err = maximum(abs.((data.cs_area_mid .- cs_area_mid_exact)./(csm_max - csm_min))) +@test csm_err < 0.00664 + +alm_min = reshape(minimum(axial_loading_mid_exact; dims=2), num_times, 1, num_blades) +alm_max = reshape(maximum(axial_loading_mid_exact; dims=2), num_times, 1, num_blades) +alm_err = maximum(abs.((data.axial_loading_mid .- axial_loading_mid_exact)./(alm_max .- alm_min))) +@test alm_err < 0.069 + +clm_min = reshape(minimum(circum_loading_mid_exact; dims=2), num_times, 1, num_blades) +clm_max = reshape(maximum(circum_loading_mid_exact; dims=2), num_times, 1, num_blades) +clm_err = maximum(abs.((data.circum_loading_mid .- circum_loading_mid_exact)./(clm_max .- clm_min))) +@test clm_err < 0.063 diff --git a/test/openfast_reading_default.jl b/test/openfast_reading_default.jl new file mode 100644 index 00000000..a0b2c254 --- /dev/null +++ b/test/openfast_reading_default.jl @@ -0,0 +1,32 @@ +using AcousticAnalogies +using Statistics: mean +using Test + +fname = joinpath(@__DIR__, "gen_test_data", "openfast_data", "IEA-3.4-130-RWT-small.out") + +num_radial = 30 +radii = range(0.2, 1.0; length=num_radial) +data = read_openfast_file(fname, radii) +num_times = length(data.time) +num_blades = size(data.pitch, 2) +num_radial = size(data.axial_loading, 2) +@test size(data.time) == (num_times,) +@test size(data.v) == (num_times,) +@test size(data.azimuth) == (num_times,) +@test size(data.omega) == (num_times,) +@test size(data.pitch) == (num_times, num_blades) +@test size(data.axial_loading) == (num_times, num_radial, num_blades) +@test size(data.circum_loading) == (num_times, num_radial, num_blades) + +@test all(data.time .≈ (60.0:0.01:60.1)) +@test all(data.v .≈ 7) +@test all(data.pitch .≈ 0*pi/180) +# Just some coarse tests. +@test mean(data.omega) ≈ 8.126818181818182*(2*pi/60) +@test mean(data.axial_loading) ≈ 1632.1274949494953 +@test mean(data.circum_loading) ≈ -217.45748949494939 + +# Make sure the averaging of the freestream velocity and omega works. +data2 = read_openfast_file(fname, radii; average_freestream_vel=true, average_omega=true) +@test all(data2.v .≈ 7) +@test all(data2.omega .≈ 8.126818181818182*(2*pi/60)) diff --git a/test/openfast_reading_different_time_column_name.jl b/test/openfast_reading_different_time_column_name.jl new file mode 100644 index 00000000..1460b39b --- /dev/null +++ b/test/openfast_reading_different_time_column_name.jl @@ -0,0 +1,32 @@ +using AcousticAnalogies +using Statistics: mean +using Test + +fname = joinpath(@__DIR__, "gen_test_data", "openfast_data", "IEA-3.4-130-RWT-small-FooTime.out") + +num_radial = 30 +radii = range(0.2, 1.0; length=num_radial) +data = read_openfast_file(fname, radii; header_keyword="FooTime") +num_times = length(data.time) +num_blades = size(data.pitch, 2) +num_radial = size(data.axial_loading, 2) +@test size(data.time) == (num_times,) +@test size(data.v) == (num_times,) +@test size(data.azimuth) == (num_times,) +@test size(data.omega) == (num_times,) +@test size(data.pitch) == (num_times, num_blades) +@test size(data.axial_loading) == (num_times, num_radial, num_blades) +@test size(data.circum_loading) == (num_times, num_radial, num_blades) + +@test all(data.time .≈ (60.0:0.01:60.1)) +@test all(data.v .≈ 7) +@test all(data.pitch .≈ 0*pi/180) +# Just some coarse tests. +@test mean(data.omega) ≈ 8.126818181818182*(2*pi/60) +@test mean(data.axial_loading) ≈ 1632.1274949494953 +@test mean(data.circum_loading) ≈ -217.45748949494939 + +# Make sure the averaging of the freestream velocity and omega works. +data2 = read_openfast_file(fname, radii; header_keyword="FooTime", average_freestream_vel=true, average_omega=true) +@test all(data2.v .≈ 7) +@test all(data2.omega .≈ 8.126818181818182*(2*pi/60)) diff --git a/test/openfast_reading_no_units_header.jl b/test/openfast_reading_no_units_header.jl new file mode 100644 index 00000000..5c2dcf2f --- /dev/null +++ b/test/openfast_reading_no_units_header.jl @@ -0,0 +1,33 @@ +using AcousticAnalogies +using Statistics: mean +using Test + +fname = joinpath(@__DIR__, "gen_test_data", "openfast_data", "IEA-3.4-130-RWT-small-no_units.out") + +num_radial = 30 +radii = range(0.2, 1.0; length=num_radial) +data = read_openfast_file(fname, radii; has_units_header=false) +num_times = length(data.time) +num_blades = size(data.pitch, 2) +num_radial = size(data.axial_loading, 2) +@test size(data.time) == (num_times,) +@test size(data.v) == (num_times,) +@test size(data.azimuth) == (num_times,) +@test size(data.omega) == (num_times,) +@test size(data.pitch) == (num_times, num_blades) +@test size(data.axial_loading) == (num_times, num_radial, num_blades) +@test size(data.circum_loading) == (num_times, num_radial, num_blades) + +@test all(data.time .≈ (60.0:0.01:60.1)) +@test all(data.v .≈ 7) +@test all(data.pitch .≈ 0) +@test all(data.pitch .≈ 0*pi/180) +# Just some coarse tests. +@test mean(data.omega) ≈ 8.126818181818182*(2*pi/60) +@test mean(data.axial_loading) ≈ 1632.1274949494953 +@test mean(data.circum_loading) ≈ -217.45748949494939 + +# Make sure the averaging of the freestream velocity and omega works. +data2 = read_openfast_file(fname, radii; has_units_header=false, average_freestream_vel=true, average_omega=true) +@test all(data2.v .≈ 7) +@test all(data2.omega .≈ 8.126818181818182*(2*pi/60)) diff --git a/test/openfast_time_derivatives_constant_time_step.jl b/test/openfast_time_derivatives_constant_time_step.jl new file mode 100644 index 00000000..6fb5cee4 --- /dev/null +++ b/test/openfast_time_derivatives_constant_time_step.jl @@ -0,0 +1,58 @@ +using AcousticAnalogies +using Polynomials: fit +using Test + +v = azimuth = omega = pitch = nothing +dradii = nothing +cs_area = cs_area_mid = nothing +axial_loading = circum_loading = nothing +T = 2.5 +t0 = 0.3 +fn(t, r, b) = 0.2*sin(2*pi/T*t*r*b + 0.1*pi) + 0.3*cos(4*pi/T*t*r*b + 0.2*pi) +fndot(t, r, b) = 0.2*2*pi/T*r*b*cos(2*pi/T*t*r*b + 0.1*pi) - 0.3*4*pi/T*r*b*sin(4*pi/T*t*r*b + 0.2*pi) +fc(t, r, b) = 0.4*sin(2*pi/T*t*r*b + 0.3*pi) + 0.5*cos(4*pi/T*t*r*b + 0.4*pi) +fcdot(t, r, b) = 0.4*2*pi/T*r*b*cos(2*pi/T*t*r*b + 0.3*pi) - 0.5*4*pi/T*r*b*sin(4*pi/T*t*r*b + 0.4*pi) +# r = reshape(range(1.0, 2.0; length=5), 1, :) +radii = range(1.0, 2.0; length=5) +radii_mid = 0.5 .* (@view(radii[1:end-1]) .+ @view(radii[2:end])) +r = reshape(radii_mid, 1, :) +b = reshape(1:3, 1, 1, :) +errs_fn_l2 = Vector{Float64}() +errs_fc_l2 = Vector{Float64}() +dts = Vector{Float64}() +for N in 120:10:150 + dt = T/N + push!(dts, dt) + time = t0 .+ (0:(N-1)) .* dt + dtime_dtau = similar(time) + axial_loading_mid = fn.(time, r, b) + circum_loading_mid = fc.(time, r, b) + axial_loading_mid_dot = similar(axial_loading_mid) + circum_loading_mid_dot = similar(circum_loading_mid) + # data = OpenFASTData{SecondOrderFiniteDiff}(time, v, azimuth, omega, pitch, axial_loading_mid, circum_loading_mid) + data = OpenFASTData{FLOWLinearInterp,SecondOrderFiniteDiff}( + time, dtime_dtau, v, azimuth, omega, pitch, + radii, radii_mid, dradii, + cs_area, cs_area_mid, + axial_loading, axial_loading_mid, axial_loading_mid_dot, + circum_loading, circum_loading_mid, circum_loading_mid_dot) + AcousticAnalogies.calculate_loading_dot!(data) + @test all(data.dtime_dtau .≈ dt) + err_fn_l2 = sqrt(sum((data.axial_loading_mid_dot .- fndot.(time, r, b)).^2) / length(data.axial_loading_mid_dot)) + push!(errs_fn_l2, err_fn_l2) + err_fc_l2 = sqrt(sum((data.circum_loading_mid_dot .- fcdot.(time, r, b)).^2) / length(data.circum_loading_mid_dot)) + push!(errs_fc_l2, err_fc_l2) +end + +# err ≈ dt^p +# err2/err1 ≈ (dt2^p)/(dt1^p) = (dt2/dt1)^p +# log(err2/err1) ≈ log((dt2/dt1)^p) = p*log(dt2/dt1) +# p ≈ log(err2/err1)/log(dt2/dt1) ≈ (log(err2) - log(err1))/(log(dt2) - log(dt1)) + +# Fit a line through the errors on a log-log plot, then check that the slope is 2 (second-order). +l_fn = fit(log.(dts), log.(errs_fn_l2), 1) +@test isapprox(l_fn.coeffs[2], 2, atol=0.02) + +# Fit a line through the errors on a log-log plot, then check that the slope is 2 (second-order). +l_fc = fit(log.(dts), log.(errs_fc_l2), 1) +@test isapprox(l_fc.coeffs[2], 2, atol=0.02) diff --git a/test/openfast_time_derivatives_nonconstant_time_step.jl b/test/openfast_time_derivatives_nonconstant_time_step.jl new file mode 100644 index 00000000..8bc75e1c --- /dev/null +++ b/test/openfast_time_derivatives_nonconstant_time_step.jl @@ -0,0 +1,58 @@ +using AcousticAnalogies +using Polynomials: fit +using Test + +v = azimuth = omega = pitch = nothing +cs_area = cs_area_mid = nothing +axial_loading = circum_loading = nothing +dradii = nothing +T = 2.5 +t0 = 0.3 +fn(t, r, b) = 0.2*sin(2*pi/T*t*r*b + 0.1*pi) + 0.3*cos(4*pi/T*t*r*b + 0.2*pi) +fndot(t, r, b) = 0.2*2*pi/T*r*b*cos(2*pi/T*t*r*b + 0.1*pi) - 0.3*4*pi/T*r*b*sin(4*pi/T*t*r*b + 0.2*pi) +fc(t, r, b) = 0.4*sin(2*pi/T*t*r*b + 0.3*pi) + 0.5*cos(4*pi/T*t*r*b + 0.4*pi) +fcdot(t, r, b) = 0.4*2*pi/T*r*b*cos(2*pi/T*t*r*b + 0.3*pi) - 0.5*4*pi/T*r*b*sin(4*pi/T*t*r*b + 0.4*pi) +# r = reshape(range(1.0, 2.0; length=5), 1, :) +radii = range(1.0, 2.0; length=5) +radii_mid = 0.5 .* (@view(radii[1:end-1]) .+ @view(radii[2:end])) +r = reshape(radii_mid, 1, :) +b = reshape(1:3, 1, 1, :) +errs_fn_l2 = Vector{Float64}() +errs_fc_l2 = Vector{Float64}() +Ns = 120:10:150 +for N in Ns + dt = T/N + time1 = t0 .+ (0:(N-1)) .* dt + wiggle = 0.49.*dt.*(cos.(2.0*pi./T.*time1)) + time = time1 .+ wiggle + dtime_dtau = similar(time) + axial_loading_mid = fn.(time, r, b) + circum_loading_mid = fc.(time, r, b) + axial_loading_mid_dot = similar(axial_loading_mid) + circum_loading_mid_dot = similar(circum_loading_mid) + # data = OpenFASTData{SecondOrderFiniteDiff}(time, v, azimuth, omega, pitch, axial_loading_mid, circum_loading_mid) + data = OpenFASTData{FLOWLinearInterp,SecondOrderFiniteDiff}( + time, dtime_dtau, v, azimuth, omega, pitch, + radii, radii_mid, dradii, + cs_area, cs_area_mid, + axial_loading, axial_loading_mid, axial_loading_mid_dot, + circum_loading, circum_loading_mid, circum_loading_mid_dot) + AcousticAnalogies.calculate_loading_dot!(data) + err_fn_l2 = sqrt(sum((data.axial_loading_mid_dot .- fndot.(time, r, b)).^2) / length(data.axial_loading_mid_dot)) + push!(errs_fn_l2, err_fn_l2) + err_fc_l2 = sqrt(sum((data.circum_loading_mid_dot .- fcdot.(time, r, b)).^2) / length(data.circum_loading_mid_dot)) + push!(errs_fc_l2, err_fc_l2) +end + +# err ≈ dt^p = (T/N)^p +# err2/err1 ≈ ((T/N2)^p)/((T/N1)^p) = ((T/N2)/(T/N1))^p = (N1/N2)^p +# log(err2/err1) ≈ log((N1/N2)^p) = p*log(N1/N2) +# p ≈ log(err2/err1)/log(N1/N2) ≈ (log(err2) - log(err1))/(log(N1) - log(N2)) = -(log(err2) - log(err1))/(log(N2) - log(N1)) + +# Fit a line through the errors on a log-log plot, then check that the slope is -2 (second-order). +l_fn = fit(log.(Ns), log.(errs_fn_l2), 1) +@test isapprox(-l_fn.coeffs[2], 2, atol=0.02) + +# Fit a line through the errors on a log-log plot, then check that the slope is -2 (second-order). +l_fc = fit(log.(Ns), log.(errs_fc_l2), 1) +@test isapprox(-l_fc.coeffs[2], 2, atol=0.02) diff --git a/test/openfast_time_derivatives_nonconstant_time_step_no_diff.jl b/test/openfast_time_derivatives_nonconstant_time_step_no_diff.jl new file mode 100644 index 00000000..55c0f8eb --- /dev/null +++ b/test/openfast_time_derivatives_nonconstant_time_step_no_diff.jl @@ -0,0 +1,38 @@ +using AcousticAnalogies +using Polynomials: fit +using Test + +v = azimuth = omega = pitch = nothing +cs_area = cs_area_mid = nothing +axial_loading = circum_loading = nothing +dradii = nothing +T = 2.5 +t0 = 0.3 +fn(t, r, b) = 0.2*sin(2*pi/T*t*r*b + 0.1*pi) + 0.3*cos(4*pi/T*t*r*b + 0.2*pi) +fndot(t, r, b) = 0.2*2*pi/T*r*b*cos(2*pi/T*t*r*b + 0.1*pi) - 0.3*4*pi/T*r*b*sin(4*pi/T*t*r*b + 0.2*pi) +fc(t, r, b) = 0.4*sin(2*pi/T*t*r*b + 0.3*pi) + 0.5*cos(4*pi/T*t*r*b + 0.4*pi) +fcdot(t, r, b) = 0.4*2*pi/T*r*b*cos(2*pi/T*t*r*b + 0.3*pi) - 0.5*4*pi/T*r*b*sin(4*pi/T*t*r*b + 0.4*pi) +# r = reshape(range(1.0, 2.0; length=5), 1, :) +radii = range(1.0, 2.0; length=5) +radii_mid = 0.5 .* (@view(radii[1:end-1]) .+ @view(radii[2:end])) +r = reshape(radii_mid, 1, :) +b = reshape(1:3, 1, 1, :) +N = 120 +dt = T/N +time1 = t0 .+ (0:(N-1)) .* dt +wiggle = 0.49.*dt.*(cos.(2.0*pi./T.*time1)) +time = time1 .+ wiggle +dtime_dtau = similar(time) +axial_loading_mid = fn.(time, r, b) +circum_loading_mid = fc.(time, r, b) +axial_loading_mid_dot = similar(axial_loading_mid) +circum_loading_mid_dot = similar(circum_loading_mid) +data = OpenFASTData{FLOWLinearInterp,NoTimeDerivMethod}( + time, dtime_dtau, v, azimuth, omega, pitch, + radii, radii_mid, dradii, + cs_area, cs_area_mid, + axial_loading, axial_loading_mid, axial_loading_mid_dot, + circum_loading, circum_loading_mid, circum_loading_mid_dot) +AcousticAnalogies.calculate_loading_dot!(data) +@test all(data.axial_loading_mid_dot .≈ 0) +@test all(data.circum_loading_mid_dot .≈ 0) diff --git a/test/runtests.jl b/test/runtests.jl index adea76f7..95806e08 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,17 +1,45 @@ module AcousticAnalogiesTests -include("adv_time_tests.jl") -include("combine_tests.jl") -include("f1a_tests.jl") -include("compact_f1a_constructor_tests.jl") -include("anopp2_comparison.jl") -include("forwarddiff_test.jl") -include("doppler_tests.jl") -include("boundary_layer_tests.jl") -include("bpm_shape_function_tests.jl") -include("broadband_source_element_tests.jl") -include("writevtk_tests.jl") -include("openfast_helper_tests.jl") -include("bpm_itr_tests.jl") +all_tests = isempty(ARGS) || ("all" in ARGS) + +if "adv_time" in ARGS || all_tests + include("adv_time_tests.jl") +end +if "combine" in ARGS || all_tests + include("combine_tests.jl") +end +if "f1a" in ARGS || all_tests + include("f1a_tests.jl") +end +if "f1a_constructor" in ARGS || all_tests + include("compact_f1a_constructor_tests.jl") +end +if "anopp2" in ARGS || all_tests + include("anopp2_comparison.jl") +end +if "forwarddiff" in ARGS || all_tests + include("forwarddiff_test.jl") +end +if "doppler" in ARGS || all_tests + include("doppler_tests.jl") +end +if "boundary_layers" in ARGS || all_tests + include("boundary_layer_tests.jl") +end +if "bpm_shape_functions" in ARGS || all_tests + include("bpm_shape_function_tests.jl") +end +if "broadband_source_elements" in ARGS || all_tests + include("broadband_source_element_tests.jl") +end +if "writevtk" in ARGS || all_tests + include("writevtk_tests.jl") +end +if "openfast" in ARGS || all_tests + include("openfast_helper_tests.jl") +end +if "bpm_itr" in ARGS || all_tests + include("bpm_itr_tests.jl") +end end # module