Adding powercurve actuator method

demo/documented/studies/powercurve_validation/CT_calibration.py

@@ -0,0 +1,146 @@
+import yaml
+
+import pandas as pd
+import numpy as np
+import matplotlib.pyplot as plt
+import seaborn
+
+################################################################################
+# this is a brief script for calibration of CT.
+#
+# data comes from a run at 7.0 m/s (region II) with all calibration factors set
+# to zero. a convergence study is run with single_turbine.yaml, and the results
+# are stored and saved in CTcalibration_data.csv.
+#
+# we target matching the reference CT at 7.0 m/s with the asymptotic simulation.
+# we assume geometric convergence of thrust w.r.t. dx. this allows the hand-fit
+# of the convergence data for an estimate of the asymptotic thrust, shown below.
+################################################################################
+
+# simulation properties
+rho_fluid = 1.225  # kg/m^3
+diameter_rotor = 130.0  # m
+area_rotor = np.pi/4*diameter_rotor**2  # m^2
+V_inf = 8.0  # m/s
+Tf_fluid = 0.5*rho_fluid*area_rotor*V_inf**2/1e3
+P_fluid = 0.5*rho_fluid*area_rotor*V_inf**3/1e6
+volume_domain = (780.0 - -390.0)*(390.0 - -390.0)*(520.0 - 0.04)
+
+# load, amend calibration simulation data
+fn_CT_data = "CTcalibration_data.csv"
+df_data = pd.read_csv(fn_CT_data)
+df_data["dx_elem"] = volume_domain/df_data.Nelem
+
+# fit loge convergence with linear model
+mb_Tf = np.polyfit(
+  np.log10(df_data.dx_elem.to_numpy()[:-1]),
+  np.log10(np.abs(rho_fluid*df_data.Tf.to_numpy()[:-1] - rho_fluid*df_data.Tf.to_numpy()[-1])),
+  deg=1,
+)
+mb_P = np.polyfit(
+  np.log10(df_data.dx_elem.to_numpy()[:-1]),
+  np.log10(np.abs(rho_fluid*df_data.P.to_numpy()[:-1] - rho_fluid*df_data.P.to_numpy()[-1])),
+  deg=1,
+)
+# i.e.:
+#   log10(|J-Jinf|) = m*log(dx) + b
+#   |J-Jinf| = dx**m * 10**b
+#   Jinf = J -/+ dx**m * 10**b
+
+# estimates of the true thrust based on each sample
+Tf_inf_est = rho_fluid*df_data.Tf + df_data.dx_elem**mb_Tf[0] * 10**mb_Tf[1]
+signpattern_P = np.ones_like(df_data.dx_elem) # np.array([1.0, 1.0, 1.0, 1.0])
+P_inf_est = rho_fluid*df_data.P + signpattern_P*df_data.dx_elem**mb_P[0] * 10**mb_P[1]
+
+# plot the convergence behavior for validation
+fig, axes = plt.subplots(2, sharex=True)
+
+# thrust vs. discretization length
+axes[0].semilogx(
+  df_data.dx_elem.to_numpy(),
+  rho_fluid*df_data.Tf.to_numpy(),
+)
+axes[0].semilogx(
+  df_data.dx_elem.to_numpy(),
+  Tf_inf_est.to_numpy(),
+  "--", label="corrected value",
+)
+axes[0].semilogx(
+  df_data.dx_elem.to_numpy(),
+  Tf_inf_est.to_numpy()[-1]*np.ones_like(df_data.dx_elem.to_numpy()),
+  "--", label="estimated asymptote",
+)
+axes[0].set_xlabel("discretization length, $\\Delta x$ (m)")
+axes[0].set_ylabel("thrust force, $T$ (kN)")
+axes[0].legend()
+
+# apparent error vs. discretization length
+axes[1].loglog(
+  df_data.dx_elem.to_numpy()[:-1],
+  np.abs(rho_fluid*df_data.Tf.to_numpy()[:-1] - rho_fluid*df_data.Tf.to_numpy()[-1]),
+)
+axes[1].loglog(
+  df_data.dx_elem.to_numpy()[:-1],
+  10**(mb_Tf[0]*np.log10(df_data.dx_elem.to_numpy()[:-1]) + mb_Tf[1]),
+  "--", label="line of best fit",
+)
+axes[1].set_xlabel("discretization length, $\\Delta x$ (m)")
+axes[1].set_ylabel("apparent thrust error, $e_T$ (kN)")
+axes[1].legend()
+
+# plot the convergence behavior for validation
+fig, axes = plt.subplots(2, sharex=True)
+
+# power vs. discretization length
+axes[0].semilogx(
+  df_data.dx_elem.to_numpy(),
+  rho_fluid*df_data.P.to_numpy(),
+)
+axes[0].semilogx(
+  df_data.dx_elem.to_numpy(),
+  P_inf_est.to_numpy(),
+  "--", label="corrected value",
+)
+axes[0].semilogx(
+  df_data.dx_elem.to_numpy(),
+  P_inf_est.to_numpy()[-1]*np.ones_like(df_data.dx_elem.to_numpy()),
+  "--", label="estimated asymptote",
+)
+axes[0].set_xlabel("discretization length, $\\Delta x$ (m)")
+axes[0].set_ylabel("power, $P$ (MW)")
+axes[0].legend()
+
+# apparent error vs. discretization length
+axes[1].loglog(
+  df_data.dx_elem.to_numpy()[:-1],
+  np.abs(rho_fluid*df_data.P.to_numpy()[:-1] - rho_fluid*df_data.P.to_numpy()[-1]),
+)
+axes[1].loglog(
+  df_data.dx_elem.to_numpy()[:-1],
+  10**(mb_P[0]*np.log10(df_data.dx_elem.to_numpy()[:-1]) + mb_P[1]),
+  "--", label="line of best fit",
+)
+axes[1].set_xlabel("discretization length, $\\Delta x$ (m)")
+axes[1].set_ylabel("power error, $e_P$ (MW)")
+axes[1].legend()
+
+# get the reference power curve we want to match
+with open("writeup/FLORIS_IEA-3p4-130-RWT.yaml", "r") as f_yaml:
+  data_refPC = yaml.safe_load(f_yaml)
+V_refPC = np.array(data_refPC["power_thrust_table"]["wind_speed"])
+CT_refPC = np.array(data_refPC["power_thrust_table"]["thrust_coefficient"])
+P_refPC = np.array(data_refPC["power_thrust_table"]["power"])/1e3
+
+# print out the CT values we see
+print(f"estimated asymptotic CT: {Tf_inf_est.to_numpy()[-1]/Tf_fluid}")
+print(f"reference powercurve CT: {np.interp(V_inf, V_refPC, CT_refPC)}")
+print(f"correction: {np.interp(V_inf, V_refPC, CT_refPC)/(Tf_inf_est.to_numpy()[-1]/Tf_fluid)}\n")
+
+print(f"estimated asymptotic CP: {P_inf_est.to_numpy()[-1]/P_fluid}")
+print(f"reference powercurve CP: {np.interp(V_inf, V_refPC, P_refPC)/P_fluid}")
+print(f"correction: {np.interp(V_inf, V_refPC, P_refPC)/(P_inf_est.to_numpy()[-1])}")
+
+# show and exit
+plt.show()
+
+###

demo/documented/studies/powercurve_validation/CTcalibration_data.csv

@@ -0,0 +1,5 @@
+Nelem,Ndof,V,P,Tf
+63744,46048,8.0,1.404740050380867578e+00,3.039104256841761185e+02
+158074,112160,8.0,1.500207951806328710e+00,3.177386680432255730e+02
+451005,314080,8.0,1.553362153286078184e+00,3.254239347945710961e+02
+1110316,765984,8.0,1.586169749666255901e+00,3.300853155380664816e+02

demo/documented/studies/powercurve_validation/multi_turbine.yaml

@@ -0,0 +1,79 @@
+general:
+    name: 3D_Wind_Farm      # Name of the output folder
+    output: ["mesh","initial_guess","height","turbine_force","solution"]
+
+wind_farm:
+
+    ######################### Grid Wind Farm #########################
+    type: grid            #                          |
+    grid_rows: 1          # Number of rows           | -
+    grid_cols: 3          # Number of columns        | -
+    x_spacing: 390        # horizontal spacing       | m      390, 520, 650, 780, 910
+    y_spacing: 0          # vertical spacing         | m
+    ###################################################################
+
+turbines:
+
+    ######################### Turbine Properties ########################
+    type: power_disk      #                           |
+    RD: 130.0             # Turbine Diameter          | m
+    thickness: 13.0       # Effective Thickness       | m
+    CPprime0: 1.022       # Region 2 Power Coeff      | -
+    # CPprime0: 1.043       # Region 2 Power Coeff      | -
+    # CPprime0: 1.360       # Region 2 Power Coeff      | -
+    CTprime0: 1.319       # Region 2 Thrust Coeff     | -
+    # CTprime0: 1.350       # Region 2 Thrust Coeff     | -
+    # CTprime0: 1.431       # Region 2 Thrust Coeff     | -
+    # CTprime0: 1.639       # Region 2 Thrust Coeff     | -
+    Prated: 2570282.74287 # Specific Rated Power      | W
+    # Prated: 2532257.89912 # Specific Rated Power      | W
+    # Prated: 3077551.02040 # Specific Rated Power      | W
+    Trated: 368163.26530  # Specific Thrust           | N
+    force: sine           # radial force distribution | -
+    HH: 110               # Hub Height                | m
+    yaw: 0.0              # Yaw                       | rads
+    ###################################################################
+
+domain:
+
+    ########################### Box Domain ############################
+    type: box              #                       |
+    x_range: [-390, 1950]   # x-range of the domain | m          1950, 2340, 2730, 3120, 3510
+    y_range: [-390, 390]   # y-range of the domain | m
+    z_range: [0.04, 520]   # z-range of the domain | m
+    nx: 54                 # Number of x-nodes     | -          54, 63, 72, 81, 90
+    ny: 18                 # Number of y-nodes     | -
+    nz: 9                  # Number of z-nodes     | -
+    ###################################################################
+
+refine:
+    warp_type: split
+    warp_percent: 0.85     # percent of cells moved       | -
+    warp_height:  240      # move cell below this value   | m
+    refine_custom:
+        1:
+            type: box
+            x_range: [-130, 1300] # 1300, 1690, 2080, 2470, 2860
+            y_range: [-130, 130]
+            z_range: [0, 240]
+    turbine_num: 1         # number of turbine refinements| -
+    turbine_factor: 1.25   # turbine radius multiplier    | -
+
+function_space:
+    type: linear
+
+boundary_conditions:
+    vel_profile: log
+    HH_vel: 8.0
+    k: 0.4
+    inflow_angle: 270.0
+
+problem:
+    type: stabilized
+    viscosity: 5
+    lmax: 50
+
+solver:
+    type: steady
+    save_power: true
+    save_thrust: true

demo/documented/studies/powercurve_validation/multi_turbine_fmt.yaml

@@ -0,0 +1,79 @@
+general:
+    name: CASENAME_TO_REPLACE # Name of the output folder
+    output: ["mesh","initial_guess","height","turbine_force","solution"]
+
+wind_farm:
+
+    ######################### Grid Wind Farm #########################
+    type: grid            #                          |
+    grid_rows: 1          # Number of rows           | -
+    grid_cols: 3          # Number of columns        | -
+    x_spacing: 910        # horizontal spacing       | m      390, 520, 650, 780, 910
+    y_spacing: 0          # vertical spacing         | m
+    ###################################################################
+
+turbines:
+
+    ######################### Turbine Properties ########################
+    type: power_disk      #                           |
+    RD: 130.0             # Turbine Diameter          | m
+    thickness: 13.0       # Effective Thickness       | m
+    CPprime0: 1.022       # Region 2 Power Coeff      | -
+    # CPprime0: 1.043       # Region 2 Power Coeff      | -
+    # CPprime0: 1.360       # Region 2 Power Coeff      | -
+    CTprime0: 1.319       # Region 2 Thrust Coeff     | -
+    # CTprime0: 1.350       # Region 2 Thrust Coeff     | -
+    # CTprime0: 1.431       # Region 2 Thrust Coeff     | -
+    # CTprime0: 1.639       # Region 2 Thrust Coeff     | -
+    Prated: 2570282.74287 # Specific Rated Power      | W
+    # Prated: 2532257.89912 # Specific Rated Power      | W
+    # Prated: 3077551.02040 # Specific Rated Power      | W
+    Trated: 368163.26530  # Specific Thrust           | N
+    force: sine           # radial force distribution | -
+    HH: 110.0             # Hub Height                | m
+    yaw: 0.0              # Yaw                       | rads
+    ###################################################################
+
+domain:
+
+    ########################### Box Domain ############################
+    type: box              #                       |
+    x_range: [-390, 3510]  # x-range of the domain | m          1950, 2340, 2730, 3120, 3510
+    y_range: [-390, 390]   # y-range of the domain | m
+    z_range: [0.04, 520]   # z-range of the domain | m
+    nx: 90                 # Number of x-nodes     | -          54, 63, 72, 81, 90
+    ny: 18                 # Number of y-nodes     | -
+    nz: 9                  # Number of z-nodes     | -
+    ###################################################################
+
+refine:
+    warp_type: split
+    warp_percent: 0.85     # percent of cells moved       | -
+    warp_height:  240      # move cell below this value   | m
+    refine_custom:
+        1:
+            type: box
+            x_range: [-130, 2860] # 1300, 1690, 2080, 2470, 2860
+            y_range: [-130, 130]
+            z_range: [0, 240]
+    turbine_num: 1         # number of turbine refinements| -
+    turbine_factor: 1.25   # turbine radius multiplier    | -
+
+function_space:
+    type: linear
+
+boundary_conditions:
+    vel_profile: log
+    HH_vel: VINF_TO_REPLACE
+    k: 0.4
+    inflow_angle: 270.0
+
+problem:
+    type: stabilized
+    viscosity: 5
+    lmax: 50
+
+solver:
+    type: steady
+    save_power: true
+    save_thrust: true

demo/documented/studies/powercurve_validation/single_turbine.yaml

@@ -0,0 +1,79 @@
+general:
+    name: 3D_Wind_Farm      # Name of the output folder
+    output: ["mesh","initial_guess","height","turbine_force","solution"]
+
+wind_farm:
+
+    ######################### Grid Wind Farm #########################
+    type: grid            #                          |
+    grid_rows: 1          # Number of rows           | -
+    grid_cols: 1          # Number of columns        | -
+    x_spacing: 0          # horizontal spacing       | m
+    y_spacing: 0          # vertical spacing         | m
+    ###################################################################
+
+turbines:
+
+    ######################### Turbine Properties ########################
+    type: power_disk      #                           |
+    RD: 130.0             # Turbine Diameter          | m
+    thickness: 13.0       # Effective Thickness       | m
+    CPprime0: 1.022       # Region 2 Power Coeff      | -
+    # CPprime0: 1.043       # Region 2 Power Coeff      | -
+    # CPprime0: 1.360       # Region 2 Power Coeff      | -
+    CTprime0: 1.319       # Region 2 Thrust Coeff     | -
+    # CTprime0: 1.350       # Region 2 Thrust Coeff     | -
+    # CTprime0: 1.431       # Region 2 Thrust Coeff     | -
+    # CTprime0: 1.639       # Region 2 Thrust Coeff     | -
+    Prated: 2570282.74287 # Specific Rated Power      | W
+    # Prated: 2532257.89912 # Specific Rated Power      | W
+    # Prated: 3077551.02040 # Specific Rated Power      | W
+    Trated: 368163.26530  # Specific Thrust           | N
+    force: sine           # radial force distribution | -
+    HH: 110               # Hub Height                | m
+    yaw: 0.0              # Yaw                       | rads
+    ###################################################################
+
+domain:
+
+    ########################### Box Domain ############################
+    type: box              #                       |
+    x_range: [-390, 780]    # x-range of the domain | m
+    y_range: [-390, 390]   # y-range of the domain | m
+    z_range: [0.04, 520]   # z-range of the domain | m
+    nx: 54 # 19 # 27 # 38 # 54  # Number of x-nodes     | -
+    ny: 36 # 13 # 18 # 25 # 36  # Number of y-nodes     | -
+    nz: 18 #  6 #  9 # 13 # 18  # Number of z-nodes     | -
+    ###################################################################
+
+refine:
+    warp_type: split
+    warp_percent: 0.85     # percent of cells moved       | -
+    warp_height:  240      # move cell below this value   | m
+    refine_custom:
+        1:
+            type: box
+            x_range: [-130, 130]
+            y_range: [-130, 130]
+            z_range: [0, 240]
+    turbine_num: 1         # number of turbine refinements| -
+    turbine_factor: 1.25   # turbine radius multiplier    | -
+
+function_space:
+    type: linear
+
+boundary_conditions:
+    vel_profile: log
+    HH_vel: 15.0 # 8.0
+    k: 0.4
+    inflow_angle: 270.0
+
+problem:
+    type: stabilized
+    viscosity: 5
+    lmax: 50
+
+solver:
+    type: steady
+    save_power: true
+    save_thrust: true