diff --git a/database/airfoils/xf-n64015a-il-100000-n5.csv b/database/airfoils/xf-n64015a-il-100000-n5.csv
new file mode 100755
index 00000000..8fd8bae1
--- /dev/null
+++ b/database/airfoils/xf-n64015a-il-100000-n5.csv
@@ -0,0 +1,132 @@
+Xfoil polar. Reynolds number fixed. Mach number fixed
+Polar key,xf-n64015a-il-100000-n5
+Airfoil,n64015a-il
+Reynolds number,100000
+Ncrit,5
+Mach,0
+Max Cl/Cd,36.8079
+Max Cl/Cd alpha,5.75
+Url,http://airfoiltools.com/polar/csv?polar=xf-n64015a-il-100000-n5
+
+Alpha,Cl,Cd,Cdp,Cm,Top_Xtr,Bot_Xtr
+-15.000,-0.7511,0.11842,0.11288,-0.0148,1.0000,0.0333
+-14.750,-0.7867,0.10562,0.09996,-0.0228,1.0000,0.0330
+-14.500,-0.8199,0.09535,0.08953,-0.0293,1.0000,0.0327
+-14.250,-0.8495,0.08708,0.08106,-0.0340,1.0000,0.0325
+-14.000,-0.8757,0.08015,0.07393,-0.0373,1.0000,0.0324
+-13.750,-0.8987,0.07420,0.06775,-0.0396,1.0000,0.0323
+-13.500,-0.9188,0.06903,0.06232,-0.0410,1.0000,0.0323
+-13.250,-0.9358,0.06449,0.05752,-0.0416,1.0000,0.0324
+-13.000,-0.9498,0.06046,0.05322,-0.0415,1.0000,0.0325
+-12.750,-0.9606,0.05686,0.04931,-0.0410,1.0000,0.0327
+-12.500,-0.9681,0.05362,0.04577,-0.0399,1.0000,0.0329
+-12.250,-0.9722,0.05071,0.04254,-0.0387,1.0000,0.0332
+-12.000,-0.9719,0.04808,0.03963,-0.0372,1.0000,0.0336
+-11.750,-0.9620,0.04587,0.03732,-0.0364,1.0000,0.0341
+-11.500,-0.9522,0.04404,0.03540,-0.0354,1.0000,0.0348
+-11.250,-0.9424,0.04240,0.03366,-0.0343,1.0000,0.0357
+-11.000,-0.9316,0.04079,0.03191,-0.0331,1.0000,0.0368
+-10.750,-0.9186,0.03916,0.03008,-0.0319,1.0000,0.0381
+-10.500,-0.9023,0.03751,0.02823,-0.0310,1.0000,0.0394
+-10.250,-0.8834,0.03595,0.02644,-0.0301,1.0000,0.0405
+-10.000,-0.8641,0.03443,0.02495,-0.0295,1.0000,0.0417
+-9.750,-0.8476,0.03317,0.02369,-0.0286,1.0000,0.0431
+-9.500,-0.8321,0.03203,0.02250,-0.0274,1.0000,0.0451
+-9.250,-0.8163,0.03096,0.02132,-0.0261,1.0000,0.0476
+-9.000,-0.8020,0.02982,0.02015,-0.0246,1.0000,0.0499
+-8.750,-0.7908,0.02879,0.01915,-0.0227,1.0000,0.0522
+-8.500,-0.7799,0.02784,0.01818,-0.0206,1.0000,0.0549
+-8.250,-0.7688,0.02702,0.01725,-0.0183,1.0000,0.0584
+-8.000,-0.7616,0.02606,0.01638,-0.0156,1.0000,0.0620
+-7.750,-0.7532,0.02530,0.01560,-0.0128,1.0000,0.0665
+-7.500,-0.7476,0.02453,0.01485,-0.0095,1.0000,0.0713
+-7.250,-0.7432,0.02385,0.01420,-0.0060,1.0000,0.0772
+-7.000,-0.7404,0.02323,0.01360,-0.0022,1.0000,0.0840
+-6.750,-0.7294,0.02252,0.01292,-0.0001,0.9968,0.0951
+-6.500,-0.7004,0.02138,0.01195,-0.0017,0.9866,0.1203
+-6.250,-0.6737,0.02012,0.01099,-0.0030,0.9760,0.1690
+-6.000,-0.6517,0.01859,0.01001,-0.0037,0.9646,0.2620
+-5.750,-0.6341,0.01723,0.00949,-0.0029,0.9524,0.3973
+-5.500,-0.6063,0.01690,0.00945,-0.0031,0.9429,0.4872
+-5.250,-0.5777,0.01681,0.00936,-0.0032,0.9324,0.5301
+-5.000,-0.5490,0.01678,0.00928,-0.0033,0.9223,0.5593
+-4.750,-0.5182,0.01679,0.00920,-0.0037,0.9137,0.5842
+-4.500,-0.4915,0.01690,0.00927,-0.0031,0.9029,0.6069
+-4.250,-0.4622,0.01708,0.00941,-0.0028,0.8947,0.6293
+-4.000,-0.4353,0.01720,0.00948,-0.0022,0.8845,0.6456
+-3.750,-0.4070,0.01723,0.00942,-0.0020,0.8756,0.6573
+-3.500,-0.3802,0.01714,0.00920,-0.0018,0.8665,0.6680
+-3.250,-0.3515,0.01714,0.00914,-0.0017,0.8575,0.6751
+-3.000,-0.3248,0.01698,0.00884,-0.0016,0.8491,0.6836
+-2.750,-0.2967,0.01694,0.00875,-0.0015,0.8402,0.6888
+-2.500,-0.2691,0.01682,0.00853,-0.0015,0.8324,0.6956
+-2.250,-0.2429,0.01672,0.00836,-0.0012,0.8234,0.7017
+-2.000,-0.2145,0.01665,0.00822,-0.0012,0.8161,0.7071
+-1.750,-0.1892,0.01655,0.00804,-0.0010,0.8077,0.7143
+-1.500,-0.1610,0.01649,0.00795,-0.0009,0.8004,0.7195
+-1.250,-0.1341,0.01645,0.00787,-0.0008,0.7927,0.7253
+-1.000,-0.1082,0.01637,0.00772,-0.0005,0.7851,0.7326
+-0.750,-0.0800,0.01636,0.00770,-0.0005,0.7786,0.7378
+-0.500,-0.0541,0.01633,0.00766,-0.0003,0.7709,0.7447
+-0.250,-0.0262,0.01629,0.00757,-0.0002,0.7652,0.7509
+0.000,0.0000,0.01632,0.00764,0.0000,0.7572,0.7572
+0.250,0.0263,0.01629,0.00757,0.0002,0.7509,0.7652
+0.500,0.0541,0.01633,0.00766,0.0003,0.7447,0.7709
+0.750,0.0801,0.01636,0.00770,0.0005,0.7378,0.7786
+1.000,0.1082,0.01637,0.00772,0.0005,0.7326,0.7851
+1.250,0.1341,0.01645,0.00787,0.0008,0.7253,0.7927
+1.500,0.1610,0.01649,0.00795,0.0009,0.7195,0.8005
+1.750,0.1892,0.01655,0.00804,0.0010,0.7143,0.8077
+2.000,0.2145,0.01665,0.00822,0.0012,0.7071,0.8161
+2.250,0.2429,0.01672,0.00836,0.0012,0.7017,0.8235
+2.500,0.2691,0.01682,0.00853,0.0015,0.6956,0.8324
+2.750,0.2967,0.01694,0.00875,0.0015,0.6888,0.8402
+3.000,0.3248,0.01698,0.00884,0.0016,0.6836,0.8491
+3.250,0.3515,0.01713,0.00914,0.0017,0.6751,0.8575
+3.500,0.3802,0.01714,0.00920,0.0018,0.6680,0.8665
+3.750,0.4070,0.01723,0.00942,0.0020,0.6573,0.8756
+4.000,0.4353,0.01720,0.00947,0.0022,0.6456,0.8845
+4.250,0.4622,0.01707,0.00941,0.0028,0.6293,0.8947
+4.500,0.4915,0.01690,0.00927,0.0031,0.6069,0.9030
+4.750,0.5182,0.01678,0.00920,0.0037,0.5842,0.9137
+5.000,0.5490,0.01678,0.00928,0.0033,0.5593,0.9224
+5.250,0.5777,0.01681,0.00936,0.0032,0.5300,0.9324
+5.500,0.6063,0.01690,0.00944,0.0031,0.4872,0.9429
+5.750,0.6342,0.01723,0.00949,0.0029,0.3973,0.9525
+6.000,0.6518,0.01859,0.01000,0.0036,0.2619,0.9646
+6.250,0.6738,0.02012,0.01099,0.0030,0.1690,0.9760
+6.500,0.7005,0.02138,0.01194,0.0017,0.1202,0.9867
+6.750,0.7295,0.02252,0.01291,0.0000,0.0951,0.9969
+7.000,0.7404,0.02322,0.01359,0.0022,0.0840,1.0000
+7.250,0.7432,0.02385,0.01420,0.0060,0.0772,1.0000
+7.500,0.7476,0.02453,0.01484,0.0095,0.0713,1.0000
+7.750,0.7533,0.02529,0.01560,0.0128,0.0665,1.0000
+8.000,0.7617,0.02605,0.01637,0.0155,0.0619,1.0000
+8.250,0.7690,0.02701,0.01725,0.0182,0.0583,1.0000
+8.500,0.7802,0.02784,0.01818,0.0205,0.0549,1.0000
+8.750,0.7911,0.02879,0.01915,0.0227,0.0522,1.0000
+9.000,0.8025,0.02982,0.02015,0.0245,0.0499,1.0000
+9.250,0.8168,0.03096,0.02132,0.0260,0.0475,1.0000
+9.500,0.8327,0.03203,0.02250,0.0273,0.0451,1.0000
+9.750,0.8482,0.03317,0.02370,0.0285,0.0431,1.0000
+10.000,0.8647,0.03443,0.02496,0.0294,0.0417,1.0000
+10.250,0.8841,0.03596,0.02645,0.0300,0.0405,1.0000
+10.500,0.9029,0.03752,0.02824,0.0309,0.0394,1.0000
+10.750,0.9192,0.03916,0.03009,0.0318,0.0381,1.0000
+11.000,0.9323,0.04080,0.03192,0.0329,0.0368,1.0000
+11.250,0.9431,0.04241,0.03367,0.0341,0.0357,1.0000
+11.500,0.9530,0.04404,0.03541,0.0352,0.0348,1.0000
+11.750,0.9629,0.04588,0.03733,0.0362,0.0341,1.0000
+12.000,0.9729,0.04811,0.03965,0.0371,0.0336,1.0000
+12.250,0.9729,0.05073,0.04258,0.0385,0.0332,1.0000
+12.500,0.9688,0.05366,0.04582,0.0398,0.0329,1.0000
+12.750,0.9614,0.05690,0.04936,0.0408,0.0327,1.0000
+13.000,0.9507,0.06050,0.05327,0.0414,0.0325,1.0000
+13.250,0.9367,0.06454,0.05758,0.0414,0.0324,1.0000
+13.500,0.9196,0.06911,0.06241,0.0408,0.0323,1.0000
+13.750,0.8995,0.07430,0.06785,0.0394,0.0323,1.0000
+14.000,0.8764,0.08028,0.07406,0.0370,0.0324,1.0000
+14.250,0.8501,0.08724,0.08122,0.0336,0.0325,1.0000
+14.500,0.8206,0.09555,0.08973,0.0289,0.0327,1.0000
+14.750,0.7872,0.10590,0.10025,0.0224,0.0330,1.0000
+15.000,0.7514,0.11884,0.11330,0.0143,0.0333,1.0000
diff --git a/docs/make.jl b/docs/make.jl
index 70a21bf5..5d50e5ed 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -50,6 +50,7 @@ makedocs(
"examples/blownwing-aero.md",
# "examples/blownwing-acoustics.md",
"examples/blownwing-asm.md",
+ "examples/prowim-aero.md",
],
"eVTOL Aircraft" => [
"examples/vahana-vehicle.md",
diff --git a/docs/src/api/flowunsteady-openvsp.md b/docs/src/api/flowunsteady-openvsp.md
index 1eb962d7..955c22b2 100644
--- a/docs/src/api/flowunsteady-openvsp.md
+++ b/docs/src/api/flowunsteady-openvsp.md
@@ -1,4 +1,4 @@
-# [Importing OpenVSP geometry](@id openvsp_import)
+# [OpenVSP geometry](@id openvsp_import)
```@docs
FLOWUnsteady.read_degengeom
diff --git a/docs/src/examples/blownwing-aero.md b/docs/src/examples/blownwing-aero.md
index a9ed7012..106b67a2 100644
--- a/docs/src/examples/blownwing-aero.md
+++ b/docs/src/examples/blownwing-aero.md
@@ -1,4 +1,4 @@
-# [Aerodynamic Solution](@id blownwingaero)
+# [Wing-on-Prop Interactions](@id blownwingaero)
```@raw html
@@ -10,6 +10,17 @@
```
+In this example we show mount propellers on a swept wing.
+The wing is modeled using the actuator line model that represents the wing
+as a lifting line.
+This wing model is accurate for capturing wing-on-prop interactions.
+For instance, the rotor will experience an unsteady blade loading (and
+increased tonal noise) caused by the turning of the flow ahead of the wing
+leading edge.
+However, this simple wing model is not adecuate for capturing
+prop-on-wing interactions (see [the next two sections](@ref asm) to
+accurately predict prop-on-wing interactions).
+
```julia
#=##############################################################################
@@ -199,8 +210,8 @@ for ri in 1:nrotors
# Account for angle of attack of wing
nrm = sqrt(x^2 + z^2)
- x = nrm*cosd(AOAwing)
- z = -nrm*sind(AOAwing)
+ x = (x==0 ? 1 : sign(x))*nrm*cosd(AOAwing)
+ z = -(z==0 ? 1 : sign(z))*nrm*sind(AOAwing)
# Translate rotor to its position along wing
O = [x, y, z] # New position
diff --git a/docs/src/examples/blownwing-asm.md b/docs/src/examples/blownwing-asm.md
index 5c05be38..c01356fe 100644
--- a/docs/src/examples/blownwing-asm.md
+++ b/docs/src/examples/blownwing-asm.md
@@ -17,13 +17,17 @@ vorticity at the three-quarter chord, as shown here:
The ALM is very accurate for isolated wings and even cases with mild wake
interactions.
However, for cases with stronger wake interactions (e.g., a wake directly
-impinging on the wing surface), we have developed an actuator surface model
-(ASM) that introduces the surface vorticity into the LES domain that better
-represents the physics.
-This is done by spreading the surface vorticity following a pressure-like
-distribution, which ends up producing a velocity field at the wing surface
-that minimizes the flow that crosses the airfoil centerline, thus better
-representing a solid surface:
+impinging on the wing surface), the ALM can lead to unphysical results as
+the flow tends to cross the airfoil centerline.
+To address this, we have developed an actuator surface model
+(ASM) to embed the wing surface in the LES domain and better
+represent the physics.
+
+The ASM spreads the surface vorticity following a pressure-like
+distribution.
+This produces a velocity field at the wing surface that minimizes the mass
+flow that crosses the airfoil centerline, thus better representing a solid
+surface:
```@raw html
@@ -32,9 +36,11 @@ representing a solid surface:
```
-For an in-depth discussion of the actuator line and surface models
-implemented in FLOWUnsteady, see Chapter 6 in
-[Alvarez' Dissertation](https://scholarsarchive.byu.edu/etd/9589).[^2]
+For an in-depth discussion of the actuator models implemented in
+FLOWUnsteady, see Chapter 6 in
+[Alvarez' Dissertation](https://scholarsarchive.byu.edu/etd/9589)[^2]
+(also published in
+[Alvarez & Ning, 2023](https://arc.aiaa.org/doi/abs/10.2514/1.C037279)[^3]).
[^2]: E. J. Alvarez (2022), "Reformulated Vortex Particle Method and
@@ -42,6 +48,10 @@ implemented in FLOWUnsteady, see Chapter 6 in
Dissertation, Brigham Young University*.
[**[VIDEO]**](https://www.nas.nasa.gov/pubs/ams/2022/08-09-22.html)
[**[PDF]**](https://scholarsarchive.byu.edu/etd/9589/)
+[^3]: E. J. Alvarez and A. Ning (2023), "Meshless Large-Eddy Simulation of
+ Propeller–Wing Interactions with Reformulated Vortex Particle Method,"
+ *Journal of Aircraft*.
+ [**[DOI]**](https://arc.aiaa.org/doi/abs/10.2514/1.C037279)[**[PDF]**](https://scholarsarchive.byu.edu/facpub/6902/)
In order to activate the actuator surface model, we define the following
parameters:
@@ -74,7 +84,7 @@ include_unsteadyforce = true # Include unsteady force
add_unsteadyforce = false # Whether to add the unsteady force to Ftot or to simply output it
include_parasiticdrag = true # Include parasitic-drag force
-add_skinfriction = true # If false, the parasitic drag is purely parasitic, meaning no skin friction
+add_skinfriction = true # If false, the parasitic drag is purely form, meaning no skin friction
calc_cd_from_cl = false # Whether to calculate cd from cl or effective AOA
wing_polar_file = "xf-rae101-il-1000000.csv" # Airfoil polar for parasitic drag
```
@@ -88,7 +98,7 @@ that uses the vortex sheet:
forces = []
# Calculate Kutta-Joukowski force
-kuttajoukowski = uns.generate_calc_aerodynamicforce_kuttajoukowski(KJforce_type,
+kuttajoukowski = uns.generate_aerodynamicforce_kuttajoukowski(KJforce_type,
sigma_vlm_surf, sigma_rotor_surf,
vlm_vortexsheet, vlm_vortexsheet_overlap,
vlm_vortexsheet_distribution,
@@ -173,10 +183,8 @@ uns.run_simulation( ...
)
```
-!!! info "ASM and High Fidelity"
- ASM uses a very high density of particles at the wing
- surface (~100k particles per wing) to accuratelly introduce the solid
- boundary into the LES.
- This increases the computational cost of the simulation considerably.
- Hence, we recommend using ASM only for high-fidelity simulations.
+!!! info "ASM Example"
+ The [next section](@ref prowimaero) shows an example on how to
+ set up and run a simulation using the actuator surface model.
+
diff --git a/docs/src/examples/prowim-aero.md b/docs/src/examples/prowim-aero.md
new file mode 100644
index 00000000..2a5833b4
--- /dev/null
+++ b/docs/src/examples/prowim-aero.md
@@ -0,0 +1,673 @@
+# [Prop-on-Wing Interactions](@id prowimaero)
+
+In this example we use the [actuator surface model](@ref asm) (ASM) to
+more accurately predict the effects of props blowing on a wing.
+This case simulates the PROWIM experiment in
+[Leo Veldhuis' dissertation](https://repository.tudelft.nl/islandora/object/uuid%3A8ffbde9c-b483-40de-90e0-97095202fbe3)
+(2005), and reproduces the validation study published in
+[Alvarez & Ning (2023)](https://arc.aiaa.org/doi/10.2514/1.C037279).
+
+In this example you can vary the fidelity of the simulation setting the
+following parameters:
+
+| Parameter | Mid-low fidelity | Mid-high fidelity | High fidelity | Description |
+| :-------: | :--------------: | :---------------: | :-----------: | :---------- |
+| `n_wing` | `50` | `50` | `100` | Number of wing elements per semispan |
+| `n_rotor` | `12` | `20` | `50` | Number of blade elements per blade |
+| `nsteps_per_rev` | `36` | `36` | `72` | Time steps per revolution |
+| `p_per_step` | `2` | `5` | `5` | Particle sheds per time step |
+| `shed_starting` | `false` | `false` | `true` | Whether to shed starting vortex |
+| `shed_unsteady` | `false` | `false` | `true` | Whether to shed vorticity from unsteady loading |
+| `treat_wake` | `true` | `true` | `false` | Treat wake to avoid instabilities |
+| `vlm_vortexsheet_overlap` | `2.125/10` | `2.125/10` | `2.125` | Particle overlap in ASM vortex sheet |
+| `vpm_integration` | `vpm.euler` | RK3``^\star`` | RK3``^\star`` | VPM time integration scheme |
+| `vpm_SFS` | None``^\dag`` | Dynamic``^\ddag`` | Dynamic``^\ddag`` | VPM LES subfilter-scale model |
+
+* ``^\star``*RK3:* `vpm_integration = vpm.rungekutta3`
+* ``^\dag``*None:* `vpm_SFS = vpm.SFS_none`
+* ``^\ddag``*Dynamic:* `vpm_SFS = vpm.SFS_Cd_twolevel_nobackscatter`
+
+(Mid-low fidelity settings may be inadequate for capturing prop-on-wing interactions, unless using `p_per_step=5`)
+
+
+
+```julia
+#=##############################################################################
+# DESCRIPTION
+ Validation of prop-on-wing interactions with twin props mounted mid span
+ blowing on a wing. This case simulates the PROWIM experiment in Leo
+ Veldhuis' dissertation (2005), “Propeller Wing Aerodynamic Interference.”
+
+ In this simulation we use the actuator surface model for the wing in order
+ to accurately capture prop-on-wing interactional effects. The rotors still
+ use the actuator line model.
+
+ The high-fidelity settings replicate the results presented in Alvarez &
+ Ning (2023), "Meshless Large-Eddy Simulation of Propeller–Wing Interactions
+ with Reformulated Vortex Particle Method," Sec. IV.B, also available in
+ Alvarez (2022), "Reformulated Vortex Particle Method and Meshless Large Eddy
+ Simulation of Multirotor Aircraft," Sec. 8.4.
+
+# ABOUT
+ * Author : Eduardo J. Alvarez (edoalvarez.com)
+ * Email : Edo.AlvarezR@gmail.com
+ * Created : January 2024
+ * Last updated : January 2024
+ * License : MIT
+=###############################################################################
+
+
+import FLOWUnsteady as uns
+import FLOWUnsteady: vlm, vpm
+
+run_name = "prowim" # Name of this simulation
+save_path = run_name*"-example2" # Where to save this simulation
+prompt = true # Whether to prompt the user
+paraview = true # Whether to visualize with Paraview
+
+add_wing = true # Whether to add wing to simulation
+add_rotors = true # Whether to add rotors to simulation
+
+# ----------------- GEOMETRY PARAMETERS ----------------------------------------
+# Wing geometry
+b = 2*0.64 # (m) span length
+ar = 5.33 # Aspect ratio b/c_tip
+tr = 1.0 # Taper ratio c_tip/c_root
+twist_root = 0.0 # (deg) twist at root
+twist_tip = 0.0 # (deg) twist at tip
+lambda = 0.0 # (deg) sweep
+gamma = 0.0 # (deg) dihedral
+thickness_w = 0.15 # Thickness t/c of wing airfoil
+
+# Rotor geometry
+rotor_file = "beaver.csv" # Rotor geometry
+data_path = uns.def_data_path # Path to rotor database
+read_polar = vlm.ap.read_polar2 # What polar reader to use
+pitch = 2.5 # (deg) collective pitch of blades
+xfoil = false # Whether to run XFOIL
+ncrit = 6 # Turbulence criterion for XFOIL
+
+# Read radius of this rotor and number of blades
+R, B = uns.read_rotor(rotor_file; data_path=data_path)[[1,3]]
+
+# Vehicle assembly
+AOAwing = 0.0 # (deg) wing angle of attack
+spanpos = [-0.46875, 0.46875] # Semi-span position of each rotor, 2*y/b
+xpos = [-0.8417, -0.8417] # x-position of rotors relative to LE, x/c
+zpos = [0.0, 0.0] # z-position of rotors relative to LE, z/c
+CWs = [false, true] # Rotation direction of each rotor: outboard up
+# CWs = [true, false] # Rotation direction of each rotor: inboard up
+nrotors = length(spanpos) # Number of rotors
+
+# Discretization
+n_wing = 50 # Number of spanwise elements per side
+r_wing = 2.0 # Geometric expansion of elements
+# n_rotor = 20 # Number of blade elements per blade
+n_rotor = 12
+r_rotor = 1/10 # Geometric expansion of elements
+
+# Check that we declared all the inputs that we need for each rotor
+@assert nrotors==length(spanpos)==length(xpos)==length(zpos)==length(CWs) ""*
+ "Invalid rotor inputs! Check that spanpos, xpos, zpos, and CWs have the same length"
+
+# ----------------- SIMULATION PARAMETERS --------------------------------------
+# Freestream
+magVinf = 49.5 # (m/s) freestream velocity
+AOA = 4.0 # (deg) vehicle angle of attack
+rho = 1.225 # (kg/m^3) air density
+mu = 1.79e-5 # (kg/ms) air dynamic viscosity
+speedofsound = 342.35 # (m/s) speed of sound
+qinf = 0.5*rho*magVinf^2 # (Pa) reference static pressure
+Vinf(X, t) = magVinf*[cosd(AOA), 0, sind(AOA)] # Freestream function
+
+# Rotor operation
+J = 0.85 # Advance ratio Vinf/(nD)
+RPM = 60*magVinf/(J*2*R) # RPM
+
+# Reference non-dimensional parameters
+Rec = rho * magVinf * (b/ar) / mu # Chord-based wing Reynolds number
+ReD = 2*pi*RPM/60*R * rho/mu * 2*R # Diameter-based rotor Reynolds number
+Mtip = 2*pi*RPM/60 * R / speedofsound # Tip Mach number
+
+println("""
+ Vinf: $(round(magVinf, digits=1)) m/s
+ RPM: $(RPM)
+ Mtip: $(round(Mtip, digits=3))
+ ReD: $(round(ReD, digits=0))
+ Rec: $(round(Rec, digits=0))
+""")
+
+
+# NOTE: Modify the variable `AOA` in order to change the angle of attack.
+# `AOAwing` will only change the angle of attack of the wing (while
+# leaving the propellers unaffected), while `AOA` changes the angle of
+# attack of the freestream (affecting both wing and props).
+
+# ----------------- SOLVER PARAMETERS ------------------------------------------
+
+# Aerodynamic solver
+VehicleType = uns.UVLMVehicle # Unsteady solver
+# VehicleType = uns.QVLMVehicle # Quasi-steady solver
+const_solution = VehicleType==uns.QVLMVehicle # Whether to assume that the
+ # solution is constant or not
+# Time parameters
+nrevs = 8 # Number of revolutions in simulation
+nsteps_per_rev = 36 # Time steps per revolution
+nsteps = const_solution ? 2 : nrevs*nsteps_per_rev # Number of time steps
+ttot = nsteps/nsteps_per_rev / (RPM/60) # (s) total simulation time
+
+# VPM particle shedding
+# p_per_step = 5 # Sheds per time step
+p_per_step = 2
+shed_starting = false # Whether to shed starting vortex (NOTE: starting vortex might make simulation unstable with AOA>8)
+shed_unsteady = false # Whether to shed vorticity from unsteady loading
+unsteady_shedcrit = 0.001 # Shed unsteady loading whenever circulation
+ # fluctuates by more than this ratio
+treat_wake = true # Treat wake to avoid instabilities
+max_particles = 1 # Maximum number of particles
+max_particles += add_rotors * (nrotors*((2*n_rotor+1)*B)*nsteps*p_per_step)
+max_particles += add_wing * (nsteps+1)*(2*n_wing*(p_per_step+1) + p_per_step)
+
+# Regularization
+sigma_vlm_surf = b/200 # VLM-on-VPM smoothing radius (σLBV of wing)
+sigma_rotor_surf= R/80 # Rotor-on-VPM smoothing radius (σ of rotor)
+lambda_vpm = 2.125 # VPM core overlap
+ # VPM smoothing radius (σ of wakes)
+sigma_vpm_overwrite = lambda_vpm * 2*pi*R/(nsteps_per_rev*p_per_step)
+
+# Rotor solver
+vlm_rlx = 0.3 # VLM relaxation <-- this also applied to rotors
+hubtiploss_correction = ( (0.75, 10, 0.5, 0.05), (1, 1, 1, 1.0) ) # Hub/tip correction
+# VPM solver
+# vpm_integration = vpm.rungekutta3 # VPM temporal integration scheme
+vpm_integration = vpm.euler
+
+vpm_viscous = vpm.Inviscid() # VPM viscous diffusion scheme
+ # Uncomment this to make it viscous
+# vpm_viscous = vpm.CoreSpreading(-1, -1, vpm.zeta_fmm; beta=100.0, itmax=20, tol=1e-1)
+
+vpm_SFS = vpm.SFS_none # VPM LES subfilter-scale model
+# vpm_SFS = vpm.DynamicSFS(vpm.Estr_fmm, vpm.pseudo3level_positive;
+ # alpha=0.999, maxC=1.0,
+ # clippings=[vpm.clipping_backscatter])
+
+# NOTE: By default we make this simulation inviscid since at such high Reynolds
+# number the viscous effects in the wake are actually negligible.
+# Notice that while viscous diffusion is negligible, turbulent diffusion
+# is important and non-negigible, so we have activated the subfilter-scale
+# (SFS) model.
+
+if VehicleType == uns.QVLMVehicle
+ # Mute warnings regarding potential colinear vortex filaments. This is
+ # needed since the quasi-steady solver will probe induced velocities at the
+ # lifting line of the blade
+ uns.vlm.VLMSolver._mute_warning(true)
+end
+
+println("""
+ Resolving wake for $(round(ttot*magVinf/b, digits=1)) span distances
+""")
+
+
+
+
+
+# ----------------- ACTUATOR SURFACE MODEL PARAMETERS (WING) -------------------
+
+# ---------- Vortex sheet parameters ---------------
+vlm_vortexsheet = true # Spread wing circulation as a vortex sheet (activates the ASM)
+vlm_vortexsheet_overlap = 2.125/10 # Particle overlap in vortex sheet
+vlm_vortexsheet_distribution = uns.g_pressure # Distribution of the vortex sheet
+
+vlm_vortexsheet_sigma_tbv = thickness_w*(b/ar) / 128 # Smoothing radius of trailing bound vorticity, σTBV for VLM-on-VPM
+
+vlm_vortexsheet_maxstaticparticle = 10^6 # Particles to preallocate for vortex sheet
+
+if add_wing && vlm_vortexsheet
+ max_particles += vlm_vortexsheet_maxstaticparticle
+end
+
+
+# ---------- Force calculation parameters ----------
+KJforce_type = "regular" # KJ force evaluated at middle of bound vortices
+# KJforce_type = "averaged" # KJ force evaluated at average vortex sheet
+# KJforce_type = "weighted" # KJ force evaluated at strength-weighted vortex sheet
+
+include_trailingboundvortex = false # Include trailing bound vortices in force calculations
+
+include_freevortices = false # Include free vortices in force calculation
+include_freevortices_TBVs = false # Include trailing bound vortex in free-vortex force
+
+include_unsteadyforce = true # Include unsteady force
+add_unsteadyforce = false # Whether to add the unsteady force to Ftot or to simply output it
+
+include_parasiticdrag = true # Include parasitic-drag force
+add_skinfriction = true # If false, the parasitic drag is purely form, meaning no skin friction
+calc_cd_from_cl = true # Whether to calculate cd from cl or effective AOA
+# calc_cd_from_cl = false
+
+# NOTE: We use a polar at a low Reynolds number (100k as opposed to 600k from
+# the experiment) as this particular polar better resembles the drag of
+# the tripped airfoil used in the experiment
+wing_polar_file = "xf-n64015a-il-100000-n5.csv" # Airfoil polar for parasitic drag (from airfoiltools.com)
+
+
+if include_freevortices && Threads.nthreads()==1
+ @warn("Free-vortex force calculation requested, but Julia was initiated"*
+ " with only one CPU thread. This will be extremely slow!"*
+ " Initate Julia with `-t num` where num is the number of cores"*
+ " availabe to speed up the computation.")
+end
+
+
+
+
+
+# ----------------- WAKE TREATMENT ---------------------------------------------
+
+wake_treatments = []
+
+# Remove particles by particle strength: remove particles neglibly weak, remove
+# particles potentially blown up
+rmv_strngth = 2.0 * magVinf*(b/ar)/2 * magVinf*ttot/nsteps/p_per_step # Reference strength (maxCL=2.0)
+minmaxGamma = rmv_strngth*[0.0001, 10.0] # Strength bounds (removes particles outside of these bounds)
+wake_treatment_strength = uns.remove_particles_strength( minmaxGamma[1]^2, minmaxGamma[2]^2; every_nsteps=1)
+
+if treat_wake
+ push!(wake_treatments, wake_treatment_strength)
+end
+
+
+
+
+
+# ----------------- 1) VEHICLE DEFINITION --------------------------------------
+
+# -------- Generate components
+println("Generating geometry...")
+
+# Generate wing
+wing = vlm.simpleWing(b, ar, tr, twist_root, lambda, gamma;
+ twist_tip=twist_tip, n=n_wing, r=r_wing);
+
+# Pitch wing to its angle of attack
+O = [0.0, 0.0, 0.0] # New position
+Oaxis = uns.gt.rotation_matrix2(0, -AOAwing, 0) # New orientation
+vlm.setcoordsystem(wing, O, Oaxis)
+
+# Generate rotors
+rotors = vlm.Rotor[]
+for ri in 1:nrotors
+
+ # Generate rotor
+ rotor = uns.generate_rotor(rotor_file;
+ pitch=pitch,
+ n=n_rotor, CW=CWs[ri], blade_r=r_rotor,
+ altReD=[RPM, J, mu/rho],
+ xfoil=xfoil,
+ ncrit=ncrit,
+ data_path=data_path,
+ read_polar=read_polar,
+ verbose=true,
+ verbose_xfoil=false,
+ plot_disc=false
+ );
+
+ # Simulate only one rotor if the wing is not in the simulation
+ if !add_wing
+ push!(rotors, rotor)
+ break
+ end
+
+ # Determine position along wing LE
+ y = spanpos[ri]*b/2
+ x = abs(y)*tand(lambda) + xpos[ri]*b/ar
+ z = abs(y)*tand(gamma) + zpos[ri]*b/ar
+
+ # Account for angle of attack of wing
+ nrm = sqrt(x^2 + z^2)
+ x = (x==0 ? 1 : sign(x))*nrm*cosd(AOAwing)
+ z = -(z==0 ? 1 : sign(z))*nrm*sind(AOAwing)
+
+ # Translate rotor to its position along wing
+ O_r = [x, y, z] # New position
+ Oaxis_r = uns.gt.rotation_matrix2(0, 0, 0) # New orientation
+ vlm.setcoordsystem(rotor, O_r, Oaxis_r; user=true)
+
+ push!(rotors, rotor)
+end
+
+
+# -------- Generate vehicle
+println("Generating vehicle...")
+
+# System of all FLOWVLM objects
+system = vlm.WingSystem()
+
+if add_wing
+ vlm.addwing(system, "Wing", wing)
+end
+
+if add_rotors
+ for (ri, rotor) in enumerate(rotors)
+ vlm.addwing(system, "Rotor$(ri)", rotor)
+ end
+end
+
+# System solved through VLM solver
+vlm_system = vlm.WingSystem()
+add_wing ? vlm.addwing(vlm_system, "Wing", wing) : nothing
+
+# Systems of rotors
+rotor_systems = add_rotors ? (rotors, ) : NTuple{0, Array{vlm.Rotor, 1}}()
+
+# System that will shed a VPM wake
+wake_system = vlm.WingSystem() # System that will shed a VPM wake
+add_wing ? vlm.addwing(wake_system, "Wing", wing) : nothing
+ # NOTE: Do NOT include rotor when using the quasi-steady solver
+if VehicleType != uns.QVLMVehicle && add_rotors
+ for (ri, rotor) in enumerate(rotors)
+ vlm.addwing(wake_system, "Rotor$(ri)", rotor)
+ end
+end
+
+# Pitch vehicle to its angle of attack (0 in this case since we have tilted the freestream instead)
+O = [0.0, 0.0, 0.0] # New position
+Oaxis = uns.gt.rotation_matrix2(0, 0, 0) # New orientation
+vlm.setcoordsystem(system, O, Oaxis)
+
+vehicle = VehicleType( system;
+ vlm_system=vlm_system,
+ rotor_systems=rotor_systems,
+ wake_system=wake_system
+ );
+
+
+
+
+
+
+# ------------- 2) MANEUVER DEFINITION -----------------------------------------
+# Non-dimensional translational velocity of vehicle over time
+Vvehicle(t) = [-1, 0, 0] # <---- Vehicle is traveling in the -x direction
+
+# Angle of the vehicle over time
+anglevehicle(t) = zeros(3)
+
+# RPM control input over time (RPM over `RPMref`)
+RPMcontrol(t) = 1.0
+
+angles = () # Angle of each tilting system (none)
+RPMs = add_rotors ? (RPMcontrol, ) : () # RPM of each rotor system
+
+maneuver = uns.KinematicManeuver(angles, RPMs, Vvehicle, anglevehicle)
+
+
+
+
+
+
+# ------------- 3) SIMULATION DEFINITION ---------------------------------------
+
+Vref = 0.0 # Reference velocity to scale maneuver by
+RPMref = RPM # Reference RPM to scale maneuver by
+Vinit = Vref*Vvehicle(0) # Initial vehicle velocity
+Winit = pi/180*(anglevehicle(1e-6) - anglevehicle(0))/(1e-6*ttot) # Initial angular velocity
+
+simulation = uns.Simulation(vehicle, maneuver, Vref, RPMref, ttot;
+ Vinit=Vinit, Winit=Winit);
+
+
+
+
+
+# ------------- *) AERODYNAMIC FORCES ------------------------------------------
+# Here we define the different components of aerodynamic of force that we desire
+# to capture with the wing using the actuator surface model
+
+forces = []
+
+# Calculate Kutta-Joukowski force
+kuttajoukowski = uns.generate_aerodynamicforce_kuttajoukowski(KJforce_type,
+ sigma_vlm_surf, sigma_rotor_surf,
+ vlm_vortexsheet, vlm_vortexsheet_overlap,
+ vlm_vortexsheet_distribution,
+ vlm_vortexsheet_sigma_tbv;
+ vehicle=vehicle)
+push!(forces, kuttajoukowski)
+
+# Free-vortex force
+if include_freevortices
+ freevortices = uns.generate_calc_aerodynamicforce_freevortices(
+ vlm_vortexsheet_maxstaticparticle,
+ sigma_vlm_surf,
+ vlm_vortexsheet,
+ vlm_vortexsheet_overlap,
+ vlm_vortexsheet_distribution,
+ vlm_vortexsheet_sigma_tbv;
+ Ffv=uns.Ffv_direct,
+ include_TBVs=include_freevortices_TBVs
+ )
+ push!(forces, freevortices)
+end
+
+# Force due to unsteady circulation
+if include_unsteadyforce
+ unsteady(args...; optargs...) = uns.calc_aerodynamicforce_unsteady(args...;
+ add_to_Ftot=add_unsteadyforce, optargs...)
+
+ push!(forces, unsteady)
+end
+
+# Parasatic-drag force (form drag and skin friction)
+if include_parasiticdrag
+ parasiticdrag = uns.generate_aerodynamicforce_parasiticdrag(
+ wing_polar_file;
+ read_path=joinpath(data_path, "airfoils"),
+ calc_cd_from_cl=calc_cd_from_cl,
+ add_skinfriction=add_skinfriction,
+ Mach=speedofsound!=nothing ? magVinf/speedofsound : nothing
+ )
+
+ push!(forces, parasiticdrag)
+end
+
+
+# Stitch all the forces into one function
+function calc_aerodynamicforce_fun(vlm_system, args...; per_unit_span=false, optargs...)
+
+ # Delete any previous force field
+ fieldname = per_unit_span ? "ftot" : "Ftot"
+ if fieldname in keys(vlm_system.sol)
+ pop!(vlm_system.sol, fieldname)
+ end
+
+ Ftot = nothing
+
+ for (fi, force) in enumerate(forces)
+ Ftot = force(vlm_system, args...; per_unit_span=per_unit_span, optargs...)
+ end
+
+ return Ftot
+end
+
+
+
+
+
+
+# ------------- 4) MONITORS DEFINITIONS ----------------------------------------
+
+# Generate wing monitor
+L_dir = [-sind(AOA), 0, cosd(AOA)] # Direction of lift
+D_dir = [ cosd(AOA), 0, sind(AOA)] # Direction of drag
+
+monitor_wing = uns.generate_monitor_wing(wing, Vinf, b, ar,
+ rho, qinf, nsteps;
+ calc_aerodynamicforce_fun=calc_aerodynamicforce_fun,
+ include_trailingboundvortex=include_trailingboundvortex,
+ L_dir=L_dir,
+ D_dir=D_dir,
+ save_path=save_path,
+ run_name=run_name*"-wing",
+ figname="wing monitor",
+ )
+
+# Generate rotors monitor
+monitor_rotors = uns.generate_monitor_rotors(rotors, J, rho, RPM, nsteps;
+ t_scale=RPM/60, # Scaling factor for time in plots
+ t_lbl="Revolutions", # Label for time axis
+ save_path=save_path,
+ run_name=run_name*"-rotors",
+ figname="rotors monitor",
+ )
+
+# Generate monitor of flow enstrophy (indicates numerical stability)
+monitor_enstrophy = uns.generate_monitor_enstrophy(;
+ save_path=save_path,
+ run_name=run_name,
+ figname="enstrophy monitor"
+ )
+
+# Generate monitor of SFS model coefficient Cd
+monitor_Cd = uns.generate_monitor_Cd(; save_path=save_path,
+ run_name=run_name,
+ figname="Cd monitor"
+ )
+
+
+# Concatenate monitors
+all_monitors = [monitor_enstrophy, monitor_Cd]
+
+add_wing ? push!(all_monitors, monitor_wing) : nothing
+add_rotors ? push!(all_monitors, monitor_rotors) : nothing
+
+monitors = uns.concatenate(all_monitors...)
+
+# Concatenate user-defined runtime function
+extra_runtime_function = uns.concatenate(monitors, wake_treatments...)
+
+
+
+# ------------- 5) RUN SIMULATION ----------------------------------------------
+println("Running simulation...")
+
+
+# Run simulation
+uns.run_simulation(simulation, nsteps;
+
+ # ----- SIMULATION OPTIONS -------------
+ Vinf=Vinf,
+ rho=rho, mu=mu, sound_spd=speedofsound,
+
+ # ----- SOLVERS OPTIONS ----------------
+ vpm_integration=vpm_integration,
+ vpm_viscous=vpm_viscous,
+ vpm_SFS=vpm_SFS,
+
+ p_per_step=p_per_step,
+ max_particles=max_particles,
+
+ sigma_vpm_overwrite=sigma_vpm_overwrite,
+ sigma_rotor_surf=sigma_rotor_surf,
+ sigma_vlm_surf=sigma_vlm_surf,
+
+ vlm_rlx=vlm_rlx,
+ vlm_vortexsheet=vlm_vortexsheet,
+ vlm_vortexsheet_overlap=vlm_vortexsheet_overlap,
+ vlm_vortexsheet_distribution=vlm_vortexsheet_distribution,
+ vlm_vortexsheet_sigma_tbv=vlm_vortexsheet_sigma_tbv,
+ max_static_particles=vlm_vortexsheet_maxstaticparticle,
+
+ hubtiploss_correction=hubtiploss_correction,
+
+ shed_starting=shed_starting,
+ shed_unsteady=shed_unsteady,
+ unsteady_shedcrit=unsteady_shedcrit,
+
+ extra_runtime_function=extra_runtime_function,
+
+ # ----- OUTPUT OPTIONS ------------------
+ save_path=save_path,
+ run_name=run_name,
+ prompt=prompt,
+ save_wopwopin=false, # <--- Generates input files for PSU-WOPWOP noise analysis if true
+
+ );
+
+
+
+
+
+
+# ----------------- 6) VISUALIZATION -------------------------------------------
+if paraview
+ println("Calling Paraview...")
+
+ # Files to open in Paraview
+ files = joinpath(save_path, run_name*"_pfield...xmf;")
+
+ if add_rotors
+ for ri in 1:nrotors
+ for bi in 1:B
+ global files *= run_name*"_Rotor$(ri)_Blade$(bi)_loft...vtk;"
+ end
+ end
+ end
+
+ if add_wing
+ files *= run_name*"_Wing_vlm...vtk;"
+ end
+
+ # Call Paraview
+ run(`paraview --data=$(files)`)
+
+end
+```
+```@raw html
+
+ Mid-low fidelity run time: 13 minutes a Dell Precision 7760 laptop.
+ Mid-high fidelity run time: 70 minutes a Dell Precision 7760 laptop.
+ High fidelity runtime: ~2 days on a 16-core AMD EPYC 7302 processor.
+
+
+```
+
+```@raw html
+
+
+ Mid-High Fidelity
+
+ 
+
+ High Fidelity
+
+ 
+
+
+```
+
+The favorable comparison with the experiment at $\alpha=0^\circ$ and
+$4^\circ$ confirms that ASM accurately predicts prop-on-wing
+interactions up to a moderate angle of attack. At $\alpha=10^\circ$ we suspect
+that the wing is mildly stalled in the experiment, leading to a larger discrepancy (further
+discussed in [Alvarez' Dissertation](https://scholarsarchive.byu.edu/etd/9589)[^1]
+and
+[Alvarez & Ning, 2023](https://arc.aiaa.org/doi/abs/10.2514/1.C037279)[^2]).
+
+!!! info "Source file"
+ Full example available under
+ [examples/prowim/](https://github.com/byuflowlab/FLOWUnsteady/blob/master/examples/prowim).
+
+
+
+[^1]: E. J. Alvarez (2022), "Reformulated Vortex Particle Method and
+ Meshless Large Eddy Simulation of Multirotor Aircraft," *Doctoral
+ Dissertation, Brigham Young University*.
+ [**[VIDEO]**](https://www.nas.nasa.gov/pubs/ams/2022/08-09-22.html)
+ [**[PDF]**](https://scholarsarchive.byu.edu/etd/9589/)
+[^2]: E. J. Alvarez and A. Ning (2023), "Meshless Large-Eddy Simulation of
+ Propeller–Wing Interactions with Reformulated Vortex Particle Method,"
+ *Journal of Aircraft*.
+ [**[DOI]**](https://arc.aiaa.org/doi/abs/10.2514/1.C037279)[**[PDF]**](https://scholarsarchive.byu.edu/facpub/6902/)
+
diff --git a/docs/src/generate_examples.jl b/docs/src/generate_examples.jl
index eb09ce04..40e9e98e 100644
--- a/docs/src/generate_examples.jl
+++ b/docs/src/generate_examples.jl
@@ -9,4 +9,5 @@ include("generate_examples_heaving.jl")
include("generate_examples_propeller.jl")
include("generate_examples_rotorhover.jl")
include("generate_examples_blownwing.jl")
+include("generate_examples_prowim.jl")
include("generate_examples_vahana.jl")
diff --git a/docs/src/generate_examples_blownwing.jl b/docs/src/generate_examples_blownwing.jl
index 5e14edc4..ea5160f7 100644
--- a/docs/src/generate_examples_blownwing.jl
+++ b/docs/src/generate_examples_blownwing.jl
@@ -8,7 +8,7 @@ remote_url = "https://edoalvar2.groups.et.byu.net/public/FLOWUnsteady/"
open(joinpath(output_path, output_name*"-aero.md"), "w") do fout
println(fout, """
- # [Aerodynamic Solution](@id blownwingaero)
+ # [Wing-on-Prop Interactions](@id blownwingaero)
```@raw html
@@ -20,6 +20,17 @@ open(joinpath(output_path, output_name*"-aero.md"), "w") do fout
```
+ In this example we show mount propellers on a swept wing.
+ The wing is modeled using the actuator line model that represents the wing
+ as a lifting line.
+ This wing model is accurate for capturing wing-on-prop interactions.
+ For instance, the rotor will experience an unsteady blade loading (and
+ increased tonal noise) caused by the turning of the flow ahead of the wing
+ leading edge.
+ However, this simple wing model is not adecuate for capturing
+ prop-on-wing interactions (see [the next two sections](@ref asm) to
+ accurately predict prop-on-wing interactions).
+
""")
println(fout, "```julia")
@@ -159,13 +170,17 @@ open(joinpath(output_path, output_name*"-asm.md"), "w") do fout
The ALM is very accurate for isolated wings and even cases with mild wake
interactions.
However, for cases with stronger wake interactions (e.g., a wake directly
- impinging on the wing surface), we have developed an actuator surface model
- (ASM) that introduces the surface vorticity into the LES domain that better
- represents the physics.
- This is done by spreading the surface vorticity following a pressure-like
- distribution, which ends up producing a velocity field at the wing surface
- that minimizes the flow that crosses the airfoil centerline, thus better
- representing a solid surface:
+ impinging on the wing surface), the ALM can lead to unphysical results as
+ the flow tends to cross the airfoil centerline.
+ To address this, we have developed an actuator surface model
+ (ASM) to embed the wing surface in the LES domain and better
+ represent the physics.
+
+ The ASM spreads the surface vorticity following a pressure-like
+ distribution.
+ This produces a velocity field at the wing surface that minimizes the mass
+ flow that crosses the airfoil centerline, thus better representing a solid
+ surface:
```@raw html
@@ -174,9 +189,11 @@ open(joinpath(output_path, output_name*"-asm.md"), "w") do fout
```
- For an in-depth discussion of the actuator line and surface models
- implemented in FLOWUnsteady, see Chapter 6 in
- [Alvarez' Dissertation](https://scholarsarchive.byu.edu/etd/9589).[^2]
+ For an in-depth discussion of the actuator models implemented in
+ FLOWUnsteady, see Chapter 6 in
+ [Alvarez' Dissertation](https://scholarsarchive.byu.edu/etd/9589)[^2]
+ (also published in
+ [Alvarez & Ning, 2023](https://arc.aiaa.org/doi/abs/10.2514/1.C037279)[^3]).
[^2]: E. J. Alvarez (2022), "Reformulated Vortex Particle Method and
@@ -184,6 +201,10 @@ open(joinpath(output_path, output_name*"-asm.md"), "w") do fout
Dissertation, Brigham Young University*.
[**[VIDEO]**](https://www.nas.nasa.gov/pubs/ams/2022/08-09-22.html)
[**[PDF]**](https://scholarsarchive.byu.edu/etd/9589/)
+ [^3]: E. J. Alvarez and A. Ning (2023), "Meshless Large-Eddy Simulation of
+ Propeller–Wing Interactions with Reformulated Vortex Particle Method,"
+ *Journal of Aircraft*.
+ [**[DOI]**](https://arc.aiaa.org/doi/abs/10.2514/1.C037279)[**[PDF]**](https://scholarsarchive.byu.edu/facpub/6902/)
In order to activate the actuator surface model, we define the following
parameters:
@@ -216,7 +237,7 @@ open(joinpath(output_path, output_name*"-asm.md"), "w") do fout
add_unsteadyforce = false # Whether to add the unsteady force to Ftot or to simply output it
include_parasiticdrag = true # Include parasitic-drag force
- add_skinfriction = true # If false, the parasitic drag is purely parasitic, meaning no skin friction
+ add_skinfriction = true # If false, the parasitic drag is purely form, meaning no skin friction
calc_cd_from_cl = false # Whether to calculate cd from cl or effective AOA
wing_polar_file = "xf-rae101-il-1000000.csv" # Airfoil polar for parasitic drag
```
@@ -230,7 +251,7 @@ open(joinpath(output_path, output_name*"-asm.md"), "w") do fout
forces = []
# Calculate Kutta-Joukowski force
- kuttajoukowski = uns.generate_calc_aerodynamicforce_kuttajoukowski(KJforce_type,
+ kuttajoukowski = uns.generate_aerodynamicforce_kuttajoukowski(KJforce_type,
sigma_vlm_surf, sigma_rotor_surf,
vlm_vortexsheet, vlm_vortexsheet_overlap,
vlm_vortexsheet_distribution,
@@ -315,12 +336,10 @@ open(joinpath(output_path, output_name*"-asm.md"), "w") do fout
)
```
- !!! info "ASM and High Fidelity"
- ASM uses a very high density of particles at the wing
- surface (~100k particles per wing) to accuratelly introduce the solid
- boundary into the LES.
- This increases the computational cost of the simulation considerably.
- Hence, we recommend using ASM only for high-fidelity simulations.
+ !!! info "ASM Example"
+ The [next section](@ref prowimaero) shows an example on how to
+ set up and run a simulation using the actuator surface model.
+
""")
end
diff --git a/docs/src/generate_examples_prowim.jl b/docs/src/generate_examples_prowim.jl
new file mode 100644
index 00000000..cacd0c80
--- /dev/null
+++ b/docs/src/generate_examples_prowim.jl
@@ -0,0 +1,122 @@
+# ------------- BLOWN WING EXAMPLE ---------------------------------------------
+output_name = "prowim"
+example_path = joinpath(uns.examples_path, "prowim")
+
+remote_url = "https://edoalvar2.groups.et.byu.net/public/FLOWUnsteady/"
+
+# -------- Aero Solution --------------------------------------------------------
+open(joinpath(output_path, output_name*"-aero.md"), "w") do fout
+
+ println(fout, """
+ # [Prop-on-Wing Interactions](@id prowimaero)
+
+ In this example we use the [actuator surface model](@ref asm) (ASM) to
+ more accurately predict the effects of props blowing on a wing.
+ This case simulates the PROWIM experiment in
+ [Leo Veldhuis' dissertation](https://repository.tudelft.nl/islandora/object/uuid%3A8ffbde9c-b483-40de-90e0-97095202fbe3)
+ (2005), and reproduces the validation study published in
+ [Alvarez & Ning (2023)](https://arc.aiaa.org/doi/10.2514/1.C037279).
+
+ In this example you can vary the fidelity of the simulation setting the
+ following parameters:
+
+ | Parameter | Mid-low fidelity | Mid-high fidelity | High fidelity | Description |
+ | :-------: | :--------------: | :---------------: | :-----------: | :---------- |
+ | `n_wing` | `50` | `50` | `100` | Number of wing elements per semispan |
+ | `n_rotor` | `12` | `20` | `50` | Number of blade elements per blade |
+ | `nsteps_per_rev` | `36` | `36` | `72` | Time steps per revolution |
+ | `p_per_step` | `2` | `5` | `5` | Particle sheds per time step |
+ | `shed_starting` | `false` | `false` | `true` | Whether to shed starting vortex |
+ | `shed_unsteady` | `false` | `false` | `true` | Whether to shed vorticity from unsteady loading |
+ | `treat_wake` | `true` | `true` | `false` | Treat wake to avoid instabilities |
+ | `vlm_vortexsheet_overlap` | `2.125/10` | `2.125/10` | `2.125` | Particle overlap in ASM vortex sheet |
+ | `vpm_integration` | `vpm.euler` | RK3``^\\star`` | RK3``^\\star`` | VPM time integration scheme |
+ | `vpm_SFS` | None``^\\dag`` | Dynamic``^\\ddag`` | Dynamic``^\\ddag`` | VPM LES subfilter-scale model |
+
+ * ``^\\star``*RK3:* `vpm_integration = vpm.rungekutta3`
+ * ``^\\dag``*None:* `vpm_SFS = vpm.SFS_none`
+ * ``^\\ddag``*Dynamic:* `vpm_SFS = vpm.SFS_Cd_twolevel_nobackscatter`
+
+ (Mid-low fidelity settings may be inadequate for capturing prop-on-wing interactions, unless using `p_per_step=5`)
+
+
+ """)
+
+ println(fout, "```julia")
+
+ open(joinpath(example_path, "prowim.jl"), "r") do fin
+
+ ignore = false
+
+ for l in eachline(fin)
+ if contains(l, "6) POSTPROCESSING")
+ break
+ end
+
+ # if l=="#=" || contains(l, "# Uncomment this")
+ if l=="#="
+ ignore=true
+ end
+
+ if !ignore && !contains(l, "save_code=")
+ println(fout, l)
+ end
+
+ if l=="=#" || contains(l, "# paraview = false")
+ ignore=false
+ end
+ end
+ end
+
+ println(fout, "```")
+
+ println(fout, """
+ ```@raw html
+
+ Mid-low fidelity run time: 13 minutes a Dell Precision 7760 laptop.
+ Mid-high fidelity run time: 70 minutes a Dell Precision 7760 laptop.
+ High fidelity runtime: ~2 days on a 16-core AMD EPYC 7302 processor.
+
+
+ ```
+
+ ```@raw html
+
+
+ Mid-High Fidelity
+
+ /prowimblown-compexp-midfi-composed.png)
+
+ High Fidelity
+
+ /prowimblown-compexp-hifi-composed.png)
+
+
+ ```
+
+ The favorable comparison with the experiment at \$\\alpha=0^\\circ\$ and
+ \$4^\\circ\$ confirms that ASM accurately predicts prop-on-wing
+ interactions up to a moderate angle of attack. At \$\\alpha=10^\\circ\$ we suspect
+ that the wing is mildly stalled in the experiment, leading to a larger discrepancy (further
+ discussed in [Alvarez' Dissertation](https://scholarsarchive.byu.edu/etd/9589)[^1]
+ and
+ [Alvarez & Ning, 2023](https://arc.aiaa.org/doi/abs/10.2514/1.C037279)[^2]).
+
+ !!! info "Source file"
+ Full example available under
+ [examples/prowim/](https://github.com/byuflowlab/FLOWUnsteady/blob/master/examples/prowim).
+
+
+
+ [^1]: E. J. Alvarez (2022), "Reformulated Vortex Particle Method and
+ Meshless Large Eddy Simulation of Multirotor Aircraft," *Doctoral
+ Dissertation, Brigham Young University*.
+ [**[VIDEO]**](https://www.nas.nasa.gov/pubs/ams/2022/08-09-22.html)
+ [**[PDF]**](https://scholarsarchive.byu.edu/etd/9589/)
+ [^2]: E. J. Alvarez and A. Ning (2023), "Meshless Large-Eddy Simulation of
+ Propeller–Wing Interactions with Reformulated Vortex Particle Method,"
+ *Journal of Aircraft*.
+ [**[DOI]**](https://arc.aiaa.org/doi/abs/10.2514/1.C037279)[**[PDF]**](https://scholarsarchive.byu.edu/facpub/6902/)
+ """)
+
+end
diff --git a/examples/blownwing/blownwing.jl b/examples/blownwing/blownwing.jl
index bdf0046a..3b0fcb45 100644
--- a/examples/blownwing/blownwing.jl
+++ b/examples/blownwing/blownwing.jl
@@ -206,8 +206,8 @@ for ri in 1:nrotors
# Account for angle of attack of wing
nrm = sqrt(x^2 + z^2)
- x = nrm*cosd(AOAwing)
- z = -nrm*sind(AOAwing)
+ x = (x==0 ? 1 : sign(x))*nrm*cosd(AOAwing)
+ z = -(z==0 ? 1 : sign(z))*nrm*sind(AOAwing)
# Translate rotor to its position along wing
O = [x, y, z] # New position
diff --git a/examples/prowim/prowim.jl b/examples/prowim/prowim.jl
new file mode 100644
index 00000000..2a8f6961
--- /dev/null
+++ b/examples/prowim/prowim.jl
@@ -0,0 +1,641 @@
+#=##############################################################################
+# DESCRIPTION
+ Validation of prop-on-wing interactions with twin props mounted mid span
+ blowing on a wing. This case simulates the PROWIM experiment in Leo
+ Veldhuis' dissertation (2005), “Propeller Wing Aerodynamic Interference.”
+
+ In this simulation we use the actuator surface model for the wing in order
+ to accurately capture prop-on-wing interactional effects. The rotors still
+ use the actuator line model.
+
+ The high-fidelity settings replicate the results presented in Alvarez &
+ Ning (2023), "Meshless Large-Eddy Simulation of Propeller–Wing Interactions
+ with Reformulated Vortex Particle Method," Sec. IV.B, also available in
+ Alvarez (2022), "Reformulated Vortex Particle Method and Meshless Large Eddy
+ Simulation of Multirotor Aircraft," Sec. 8.4.
+
+# ABOUT
+ * Author : Eduardo J. Alvarez (edoalvarez.com)
+ * Email : Edo.AlvarezR@gmail.com
+ * Created : January 2024
+ * Last updated : January 2024
+ * License : MIT
+=###############################################################################
+
+#=
+Use the following parameters to obtain the desired fidelity
+
+---- MID-LOW FIDELITY ------
+treat_wake = true # Treat wake to avoid instabilities
+n_wing = 50 # Number of wing elements per side
+n_rotor = 12 # Number of blade elements per blade
+nsteps_per_rev = 36 # Time steps per revolution
+p_per_step = 2 # Sheds per time step
+
+shed_starting = false # Whether to shed starting vortex
+shed_unsteady = false # Whether to shed vorticity from unsteady loading
+vlm_vortexsheet_overlap = 2.125/10 # Particle overlap in ASM vortex sheet
+
+vpm_integration = vpm.euler # VPM time integration scheme
+vpm_SFS = vpm.SFS_none # VPM LES subfilter-scale model
+
+---- MID FIDELITY ------
+treat_wake = true
+n_wing = 50
+n_rotor = 20
+nsteps_per_rev = 36
+p_per_step = 5
+
+shed_starting = false
+shed_unsteady = false
+vlm_vortexsheet_overlap = 2.125/10
+
+vpm_integration = vpm.rungekutta3
+vpm_SFS = vpm.DynamicSFS(vpm.Estr_fmm, vpm.pseudo3level_positive;
+ alpha=0.999, maxC=1.0,
+ clippings=[vpm.clipping_backscatter])
+
+---- HIGH FIDELITY -----
+treat_wake = false
+n_wing = 100
+n_rotor = 50
+nsteps_per_rev = 72
+p_per_step = 5
+
+shed_starting = AOA < 8 ? true : false
+shed_unsteady = true
+vlm_vortexsheet_overlap = 2.125
+
+vpm_integration = vpm.rungekutta3
+vpm_SFS = vpm.DynamicSFS(vpm.Estr_fmm, vpm.pseudo3level_positive;
+ alpha=0.999, maxC=1.0,
+ clippings=[vpm.clipping_backscatter])
+=#
+
+import FLOWUnsteady as uns
+import FLOWUnsteady: vlm, vpm
+
+run_name = "prowim" # Name of this simulation
+save_path = run_name*"-example2" # Where to save this simulation
+prompt = true # Whether to prompt the user
+paraview = true # Whether to visualize with Paraview
+
+add_wing = true # Whether to add wing to simulation
+add_rotors = true # Whether to add rotors to simulation
+
+# ----------------- GEOMETRY PARAMETERS ----------------------------------------
+# Wing geometry
+b = 2*0.64 # (m) span length
+ar = 5.33 # Aspect ratio b/c_tip
+tr = 1.0 # Taper ratio c_tip/c_root
+twist_root = 0.0 # (deg) twist at root
+twist_tip = 0.0 # (deg) twist at tip
+lambda = 0.0 # (deg) sweep
+gamma = 0.0 # (deg) dihedral
+thickness_w = 0.15 # Thickness t/c of wing airfoil
+
+# Rotor geometry
+rotor_file = "beaver.csv" # Rotor geometry
+data_path = uns.def_data_path # Path to rotor database
+read_polar = vlm.ap.read_polar2 # What polar reader to use
+pitch = 2.5 # (deg) collective pitch of blades
+xfoil = false # Whether to run XFOIL
+ncrit = 6 # Turbulence criterion for XFOIL
+
+# Read radius of this rotor and number of blades
+R, B = uns.read_rotor(rotor_file; data_path=data_path)[[1,3]]
+
+# Vehicle assembly
+AOAwing = 0.0 # (deg) wing angle of attack
+spanpos = [-0.46875, 0.46875] # Semi-span position of each rotor, 2*y/b
+xpos = [-0.8417, -0.8417] # x-position of rotors relative to LE, x/c
+zpos = [0.0, 0.0] # z-position of rotors relative to LE, z/c
+CWs = [false, true] # Rotation direction of each rotor: outboard up
+# CWs = [true, false] # Rotation direction of each rotor: inboard up
+nrotors = length(spanpos) # Number of rotors
+
+# Discretization
+n_wing = 50 # Number of spanwise elements per side
+r_wing = 2.0 # Geometric expansion of elements
+# n_rotor = 20 # Number of blade elements per blade
+n_rotor = 12
+r_rotor = 1/10 # Geometric expansion of elements
+
+# Check that we declared all the inputs that we need for each rotor
+@assert nrotors==length(spanpos)==length(xpos)==length(zpos)==length(CWs) ""*
+ "Invalid rotor inputs! Check that spanpos, xpos, zpos, and CWs have the same length"
+
+# ----------------- SIMULATION PARAMETERS --------------------------------------
+# Freestream
+magVinf = 49.5 # (m/s) freestream velocity
+AOA = 4.0 # (deg) vehicle angle of attack
+rho = 1.225 # (kg/m^3) air density
+mu = 1.79e-5 # (kg/ms) air dynamic viscosity
+speedofsound = 342.35 # (m/s) speed of sound
+qinf = 0.5*rho*magVinf^2 # (Pa) reference static pressure
+Vinf(X, t) = magVinf*[cosd(AOA), 0, sind(AOA)] # Freestream function
+
+# Rotor operation
+J = 0.85 # Advance ratio Vinf/(nD)
+RPM = 60*magVinf/(J*2*R) # RPM
+
+# Reference non-dimensional parameters
+Rec = rho * magVinf * (b/ar) / mu # Chord-based wing Reynolds number
+ReD = 2*pi*RPM/60*R * rho/mu * 2*R # Diameter-based rotor Reynolds number
+Mtip = 2*pi*RPM/60 * R / speedofsound # Tip Mach number
+
+println("""
+ Vinf: $(round(magVinf, digits=1)) m/s
+ RPM: $(RPM)
+ Mtip: $(round(Mtip, digits=3))
+ ReD: $(round(ReD, digits=0))
+ Rec: $(round(Rec, digits=0))
+""")
+
+
+# NOTE: Modify the variable `AOA` in order to change the angle of attack.
+# `AOAwing` will only change the angle of attack of the wing (while
+# leaving the propellers unaffected), while `AOA` changes the angle of
+# attack of the freestream (affecting both wing and props).
+
+# ----------------- SOLVER PARAMETERS ------------------------------------------
+
+# Aerodynamic solver
+VehicleType = uns.UVLMVehicle # Unsteady solver
+# VehicleType = uns.QVLMVehicle # Quasi-steady solver
+const_solution = VehicleType==uns.QVLMVehicle # Whether to assume that the
+ # solution is constant or not
+# Time parameters
+nrevs = 8 # Number of revolutions in simulation
+nsteps_per_rev = 36 # Time steps per revolution
+nsteps = const_solution ? 2 : nrevs*nsteps_per_rev # Number of time steps
+ttot = nsteps/nsteps_per_rev / (RPM/60) # (s) total simulation time
+
+# VPM particle shedding
+# p_per_step = 5 # Sheds per time step
+p_per_step = 2
+shed_starting = false # Whether to shed starting vortex (NOTE: starting vortex might make simulation unstable with AOA>8)
+shed_unsteady = false # Whether to shed vorticity from unsteady loading
+unsteady_shedcrit = 0.001 # Shed unsteady loading whenever circulation
+ # fluctuates by more than this ratio
+treat_wake = true # Treat wake to avoid instabilities
+max_particles = 1 # Maximum number of particles
+max_particles += add_rotors * (nrotors*((2*n_rotor+1)*B)*nsteps*p_per_step)
+max_particles += add_wing * (nsteps+1)*(2*n_wing*(p_per_step+1) + p_per_step)
+
+# Regularization
+sigma_vlm_surf = b/200 # VLM-on-VPM smoothing radius (σLBV of wing)
+sigma_rotor_surf= R/80 # Rotor-on-VPM smoothing radius (σ of rotor)
+lambda_vpm = 2.125 # VPM core overlap
+ # VPM smoothing radius (σ of wakes)
+sigma_vpm_overwrite = lambda_vpm * 2*pi*R/(nsteps_per_rev*p_per_step)
+
+# Rotor solver
+vlm_rlx = 0.3 # VLM relaxation <-- this also applied to rotors
+hubtiploss_correction = ( (0.75, 10, 0.5, 0.05), (1, 1, 1, 1.0) ) # Hub/tip correction
+# VPM solver
+# vpm_integration = vpm.rungekutta3 # VPM temporal integration scheme
+vpm_integration = vpm.euler
+
+vpm_viscous = vpm.Inviscid() # VPM viscous diffusion scheme
+ # Uncomment this to make it viscous
+# vpm_viscous = vpm.CoreSpreading(-1, -1, vpm.zeta_fmm; beta=100.0, itmax=20, tol=1e-1)
+
+vpm_SFS = vpm.SFS_none # VPM LES subfilter-scale model
+# vpm_SFS = vpm.DynamicSFS(vpm.Estr_fmm, vpm.pseudo3level_positive;
+ # alpha=0.999, maxC=1.0,
+ # clippings=[vpm.clipping_backscatter])
+
+# NOTE: By default we make this simulation inviscid since at such high Reynolds
+# number the viscous effects in the wake are actually negligible.
+# Notice that while viscous diffusion is negligible, turbulent diffusion
+# is important and non-negigible, so we have activated the subfilter-scale
+# (SFS) model.
+
+if VehicleType == uns.QVLMVehicle
+ # Mute warnings regarding potential colinear vortex filaments. This is
+ # needed since the quasi-steady solver will probe induced velocities at the
+ # lifting line of the blade
+ uns.vlm.VLMSolver._mute_warning(true)
+end
+
+println("""
+ Resolving wake for $(round(ttot*magVinf/b, digits=1)) span distances
+""")
+
+
+
+
+
+# ----------------- ACTUATOR SURFACE MODEL PARAMETERS (WING) -------------------
+
+# ---------- Vortex sheet parameters ---------------
+vlm_vortexsheet = true # Spread wing circulation as a vortex sheet (activates the ASM)
+vlm_vortexsheet_overlap = 2.125/10 # Particle overlap in vortex sheet
+vlm_vortexsheet_distribution = uns.g_pressure # Distribution of the vortex sheet
+
+vlm_vortexsheet_sigma_tbv = thickness_w*(b/ar) / 128 # Smoothing radius of trailing bound vorticity, σTBV for VLM-on-VPM
+
+vlm_vortexsheet_maxstaticparticle = 10^6 # Particles to preallocate for vortex sheet
+
+if add_wing && vlm_vortexsheet
+ max_particles += vlm_vortexsheet_maxstaticparticle
+end
+
+
+# ---------- Force calculation parameters ----------
+KJforce_type = "regular" # KJ force evaluated at middle of bound vortices
+# KJforce_type = "averaged" # KJ force evaluated at average vortex sheet
+# KJforce_type = "weighted" # KJ force evaluated at strength-weighted vortex sheet
+
+include_trailingboundvortex = false # Include trailing bound vortices in force calculations
+
+include_freevortices = false # Include free vortices in force calculation
+include_freevortices_TBVs = false # Include trailing bound vortex in free-vortex force
+
+include_unsteadyforce = true # Include unsteady force
+add_unsteadyforce = false # Whether to add the unsteady force to Ftot or to simply output it
+
+include_parasiticdrag = true # Include parasitic-drag force
+add_skinfriction = true # If false, the parasitic drag is purely form, meaning no skin friction
+calc_cd_from_cl = true # Whether to calculate cd from cl or effective AOA
+# calc_cd_from_cl = false
+
+# NOTE: We use a polar at a low Reynolds number (100k as opposed to 600k from
+# the experiment) as this particular polar better resembles the drag of
+# the tripped airfoil used in the experiment
+wing_polar_file = "xf-n64015a-il-100000-n5.csv" # Airfoil polar for parasitic drag (from airfoiltools.com)
+
+
+if include_freevortices && Threads.nthreads()==1
+ @warn("Free-vortex force calculation requested, but Julia was initiated"*
+ " with only one CPU thread. This will be extremely slow!"*
+ " Initate Julia with `-t num` where num is the number of cores"*
+ " availabe to speed up the computation.")
+end
+
+
+
+
+
+# ----------------- WAKE TREATMENT ---------------------------------------------
+
+wake_treatments = []
+
+# Remove particles by particle strength: remove particles neglibly weak, remove
+# particles potentially blown up
+rmv_strngth = 2.0 * magVinf*(b/ar)/2 * magVinf*ttot/nsteps/p_per_step # Reference strength (maxCL=2.0)
+minmaxGamma = rmv_strngth*[0.0001, 10.0] # Strength bounds (removes particles outside of these bounds)
+wake_treatment_strength = uns.remove_particles_strength( minmaxGamma[1]^2, minmaxGamma[2]^2; every_nsteps=1)
+
+if treat_wake
+ push!(wake_treatments, wake_treatment_strength)
+end
+
+
+
+
+
+# ----------------- 1) VEHICLE DEFINITION --------------------------------------
+
+# -------- Generate components
+println("Generating geometry...")
+
+# Generate wing
+wing = vlm.simpleWing(b, ar, tr, twist_root, lambda, gamma;
+ twist_tip=twist_tip, n=n_wing, r=r_wing);
+
+# Pitch wing to its angle of attack
+O = [0.0, 0.0, 0.0] # New position
+Oaxis = uns.gt.rotation_matrix2(0, -AOAwing, 0) # New orientation
+vlm.setcoordsystem(wing, O, Oaxis)
+
+# Generate rotors
+rotors = vlm.Rotor[]
+for ri in 1:nrotors
+
+ # Generate rotor
+ rotor = uns.generate_rotor(rotor_file;
+ pitch=pitch,
+ n=n_rotor, CW=CWs[ri], blade_r=r_rotor,
+ altReD=[RPM, J, mu/rho],
+ xfoil=xfoil,
+ ncrit=ncrit,
+ data_path=data_path,
+ read_polar=read_polar,
+ verbose=true,
+ verbose_xfoil=false,
+ plot_disc=false
+ );
+
+ # Simulate only one rotor if the wing is not in the simulation
+ if !add_wing
+ push!(rotors, rotor)
+ break
+ end
+
+ # Determine position along wing LE
+ y = spanpos[ri]*b/2
+ x = abs(y)*tand(lambda) + xpos[ri]*b/ar
+ z = abs(y)*tand(gamma) + zpos[ri]*b/ar
+
+ # Account for angle of attack of wing
+ nrm = sqrt(x^2 + z^2)
+ x = (x==0 ? 1 : sign(x))*nrm*cosd(AOAwing)
+ z = -(z==0 ? 1 : sign(z))*nrm*sind(AOAwing)
+
+ # Translate rotor to its position along wing
+ O_r = [x, y, z] # New position
+ Oaxis_r = uns.gt.rotation_matrix2(0, 0, 0) # New orientation
+ vlm.setcoordsystem(rotor, O_r, Oaxis_r; user=true)
+
+ push!(rotors, rotor)
+end
+
+
+# -------- Generate vehicle
+println("Generating vehicle...")
+
+# System of all FLOWVLM objects
+system = vlm.WingSystem()
+
+if add_wing
+ vlm.addwing(system, "Wing", wing)
+end
+
+if add_rotors
+ for (ri, rotor) in enumerate(rotors)
+ vlm.addwing(system, "Rotor$(ri)", rotor)
+ end
+end
+
+# System solved through VLM solver
+vlm_system = vlm.WingSystem()
+add_wing ? vlm.addwing(vlm_system, "Wing", wing) : nothing
+
+# Systems of rotors
+rotor_systems = add_rotors ? (rotors, ) : NTuple{0, Array{vlm.Rotor, 1}}()
+
+# System that will shed a VPM wake
+wake_system = vlm.WingSystem() # System that will shed a VPM wake
+add_wing ? vlm.addwing(wake_system, "Wing", wing) : nothing
+ # NOTE: Do NOT include rotor when using the quasi-steady solver
+if VehicleType != uns.QVLMVehicle && add_rotors
+ for (ri, rotor) in enumerate(rotors)
+ vlm.addwing(wake_system, "Rotor$(ri)", rotor)
+ end
+end
+
+# Pitch vehicle to its angle of attack (0 in this case since we have tilted the freestream instead)
+O = [0.0, 0.0, 0.0] # New position
+Oaxis = uns.gt.rotation_matrix2(0, 0, 0) # New orientation
+vlm.setcoordsystem(system, O, Oaxis)
+
+vehicle = VehicleType( system;
+ vlm_system=vlm_system,
+ rotor_systems=rotor_systems,
+ wake_system=wake_system
+ );
+
+
+
+
+
+
+# ------------- 2) MANEUVER DEFINITION -----------------------------------------
+# Non-dimensional translational velocity of vehicle over time
+Vvehicle(t) = [-1, 0, 0] # <---- Vehicle is traveling in the -x direction
+
+# Angle of the vehicle over time
+anglevehicle(t) = zeros(3)
+
+# RPM control input over time (RPM over `RPMref`)
+RPMcontrol(t) = 1.0
+
+angles = () # Angle of each tilting system (none)
+RPMs = add_rotors ? (RPMcontrol, ) : () # RPM of each rotor system
+
+maneuver = uns.KinematicManeuver(angles, RPMs, Vvehicle, anglevehicle)
+
+
+
+
+
+
+# ------------- 3) SIMULATION DEFINITION ---------------------------------------
+
+Vref = 0.0 # Reference velocity to scale maneuver by
+RPMref = RPM # Reference RPM to scale maneuver by
+Vinit = Vref*Vvehicle(0) # Initial vehicle velocity
+Winit = pi/180*(anglevehicle(1e-6) - anglevehicle(0))/(1e-6*ttot) # Initial angular velocity
+
+simulation = uns.Simulation(vehicle, maneuver, Vref, RPMref, ttot;
+ Vinit=Vinit, Winit=Winit);
+
+
+
+
+
+# ------------- *) AERODYNAMIC FORCES ------------------------------------------
+# Here we define the different components of aerodynamic of force that we desire
+# to capture with the wing using the actuator surface model
+
+forces = []
+
+# Calculate Kutta-Joukowski force
+kuttajoukowski = uns.generate_aerodynamicforce_kuttajoukowski(KJforce_type,
+ sigma_vlm_surf, sigma_rotor_surf,
+ vlm_vortexsheet, vlm_vortexsheet_overlap,
+ vlm_vortexsheet_distribution,
+ vlm_vortexsheet_sigma_tbv;
+ vehicle=vehicle)
+push!(forces, kuttajoukowski)
+
+# Free-vortex force
+if include_freevortices
+ freevortices = uns.generate_calc_aerodynamicforce_freevortices(
+ vlm_vortexsheet_maxstaticparticle,
+ sigma_vlm_surf,
+ vlm_vortexsheet,
+ vlm_vortexsheet_overlap,
+ vlm_vortexsheet_distribution,
+ vlm_vortexsheet_sigma_tbv;
+ Ffv=uns.Ffv_direct,
+ include_TBVs=include_freevortices_TBVs
+ )
+ push!(forces, freevortices)
+end
+
+# Force due to unsteady circulation
+if include_unsteadyforce
+ unsteady(args...; optargs...) = uns.calc_aerodynamicforce_unsteady(args...;
+ add_to_Ftot=add_unsteadyforce, optargs...)
+
+ push!(forces, unsteady)
+end
+
+# Parasatic-drag force (form drag and skin friction)
+if include_parasiticdrag
+ parasiticdrag = uns.generate_aerodynamicforce_parasiticdrag(
+ wing_polar_file;
+ read_path=joinpath(data_path, "airfoils"),
+ calc_cd_from_cl=calc_cd_from_cl,
+ add_skinfriction=add_skinfriction,
+ Mach=speedofsound!=nothing ? magVinf/speedofsound : nothing
+ )
+
+ push!(forces, parasiticdrag)
+end
+
+
+# Stitch all the forces into one function
+function calc_aerodynamicforce_fun(vlm_system, args...; per_unit_span=false, optargs...)
+
+ # Delete any previous force field
+ fieldname = per_unit_span ? "ftot" : "Ftot"
+ if fieldname in keys(vlm_system.sol)
+ pop!(vlm_system.sol, fieldname)
+ end
+
+ Ftot = nothing
+
+ for (fi, force) in enumerate(forces)
+ Ftot = force(vlm_system, args...; per_unit_span=per_unit_span, optargs...)
+ end
+
+ return Ftot
+end
+
+
+
+
+
+
+# ------------- 4) MONITORS DEFINITIONS ----------------------------------------
+
+# Generate wing monitor
+L_dir = [-sind(AOA), 0, cosd(AOA)] # Direction of lift
+D_dir = [ cosd(AOA), 0, sind(AOA)] # Direction of drag
+
+monitor_wing = uns.generate_monitor_wing(wing, Vinf, b, ar,
+ rho, qinf, nsteps;
+ calc_aerodynamicforce_fun=calc_aerodynamicforce_fun,
+ include_trailingboundvortex=include_trailingboundvortex,
+ L_dir=L_dir,
+ D_dir=D_dir,
+ save_path=save_path,
+ run_name=run_name*"-wing",
+ figname="wing monitor",
+ )
+
+# Generate rotors monitor
+monitor_rotors = uns.generate_monitor_rotors(rotors, J, rho, RPM, nsteps;
+ t_scale=RPM/60, # Scaling factor for time in plots
+ t_lbl="Revolutions", # Label for time axis
+ save_path=save_path,
+ run_name=run_name*"-rotors",
+ figname="rotors monitor",
+ )
+
+# Generate monitor of flow enstrophy (indicates numerical stability)
+monitor_enstrophy = uns.generate_monitor_enstrophy(;
+ save_path=save_path,
+ run_name=run_name,
+ figname="enstrophy monitor"
+ )
+
+# Generate monitor of SFS model coefficient Cd
+monitor_Cd = uns.generate_monitor_Cd(; save_path=save_path,
+ run_name=run_name,
+ figname="Cd monitor"
+ )
+
+
+# Concatenate monitors
+all_monitors = [monitor_enstrophy, monitor_Cd]
+
+add_wing ? push!(all_monitors, monitor_wing) : nothing
+add_rotors ? push!(all_monitors, monitor_rotors) : nothing
+
+monitors = uns.concatenate(all_monitors...)
+
+# Concatenate user-defined runtime function
+extra_runtime_function = uns.concatenate(monitors, wake_treatments...)
+
+
+
+# ------------- 5) RUN SIMULATION ----------------------------------------------
+println("Running simulation...")
+
+
+# Run simulation
+uns.run_simulation(simulation, nsteps;
+
+ # ----- SIMULATION OPTIONS -------------
+ Vinf=Vinf,
+ rho=rho, mu=mu, sound_spd=speedofsound,
+
+ # ----- SOLVERS OPTIONS ----------------
+ vpm_integration=vpm_integration,
+ vpm_viscous=vpm_viscous,
+ vpm_SFS=vpm_SFS,
+
+ p_per_step=p_per_step,
+ max_particles=max_particles,
+
+ sigma_vpm_overwrite=sigma_vpm_overwrite,
+ sigma_rotor_surf=sigma_rotor_surf,
+ sigma_vlm_surf=sigma_vlm_surf,
+
+ vlm_rlx=vlm_rlx,
+ vlm_vortexsheet=vlm_vortexsheet,
+ vlm_vortexsheet_overlap=vlm_vortexsheet_overlap,
+ vlm_vortexsheet_distribution=vlm_vortexsheet_distribution,
+ vlm_vortexsheet_sigma_tbv=vlm_vortexsheet_sigma_tbv,
+ max_static_particles=vlm_vortexsheet_maxstaticparticle,
+
+ hubtiploss_correction=hubtiploss_correction,
+
+ shed_starting=shed_starting,
+ shed_unsteady=shed_unsteady,
+ unsteady_shedcrit=unsteady_shedcrit,
+
+ extra_runtime_function=extra_runtime_function,
+
+ # ----- OUTPUT OPTIONS ------------------
+ save_path=save_path,
+ run_name=run_name,
+ prompt=prompt,
+ save_wopwopin=false, # <--- Generates input files for PSU-WOPWOP noise analysis if true
+ save_code= !contains(@__FILE__, "REPL") ? (@__FILE__) : ""
+
+ );
+
+
+
+
+
+
+# ----------------- 6) VISUALIZATION -------------------------------------------
+if paraview
+ println("Calling Paraview...")
+
+ # Files to open in Paraview
+ files = joinpath(save_path, run_name*"_pfield...xmf;")
+
+ if add_rotors
+ for ri in 1:nrotors
+ for bi in 1:B
+ global files *= run_name*"_Rotor$(ri)_Blade$(bi)_loft...vtk;"
+ end
+ end
+ end
+
+ if add_wing
+ files *= run_name*"_Wing_vlm...vtk;"
+ end
+
+ # Call Paraview
+ run(`paraview --data=$(files)`)
+
+end