.. _userDefinedDynamics:
.. |deg| unicode:: U+000B0 .. DEGREE SIGN
User Defined Dynamics
=============================
In Flow360, users are able to specify arbitrary dynamics. This section provides an example for a Proportional-Integral (PI) controller and an example for Aero-Structural Interaction (ASI).
.. _PIcontrollerAoa:
PI Controller for Angle of Attack to Control Lift Coefficient
--------------------------------------------------------------------------
In this example, we add a controller to the :ref:`om6Wing` case study. Based on the entries of the configuration file in this case study, the resulting value of the lift coefficient, *CL*, obtained from the simulation is ~0.26. Our objective is to add a controller to this case study such that for a given target value of *CL* (e.g., 0.4), the required angle of attack, :code:`alphaAngle` , is estimated. To this end, a PI controller is defined under :code:`userDefinedDynamics` in the configuration file as:
.. literalinclude:: Files/alphaController.json
More details for some of the parameters are explained below:
1. **dynamicsName**: This is a name given by the user for the dynamics. The results of user-defined dynamics are saved to "udd_dynamicsName_v2.csv".
2. **constants**: This includes a list of all the constants related to this UDD. These constants will be used in the equations for :code:`updateLaw` and :code:`outputVars`. In the above-mentioned example, :code:`CLTarget`, :code:`Kp`, and :code:`Ki` are the target value of the lift coefficient, the proportional gain of the controller, and the integral gain of the controller, respectively.
3. **outputVars**: Each :code:`outputVars` consists of the variable name and the relations between the variable and the state variables. For this example, we set the output variable equal to the angle of attack calculated by the controller after the first 500 pseudoSteps (i.e., :code:`state[0]`).
4. **stateVarsInitialValue**: There are two state variables for this controller. The first one, :code:`state[0]`, is the angle of attack computed by the controller, and the second one, :code:`state[1]`, is the summation of error :code:`(CLTarget - CL)` over iterations. The initial values of :code:`state[0]` and :code:`state[1]` are set to :code:`alphaAngle` and :code:`0`, respectively.
The equation for the control system can be thus written as:
.. math::
:label: def_dynamics-1
\alpha^n = \alpha^{n-1} + \left[ K_p \cdot (C_{L,target} - C_L^n) + K_i \cdot e^n \right]
e^n = e^{n-1} + C_{L,target} - C_L^n
where the superscript :math:`n` denotes pseudo step, :code:`alphaAngle` (:math:`\alpha`) is :code:`state[0]` and the integral of error (:math:`e`) is :code:`state[1]`.
5. **updateLaw**: The expressions for the state variables are specified here. The first and the second expressions correspond to equations for :code:`state[0]` and :code:`state[1]`, respectively. The conditional statement forces the controller to not be run for the first 500 pseudoSteps. This ensures that the flow field is established before the controller is initialized.
6. **inputBoundaryPatches**: The related boundary conditions for the inputs of the user-defined dynamics are specified here. In this example, the noSlipWall boundary conditions are where :code:`CL` is calculated. As seen in the configuration file, the boundary name for the No-Slip boundary condition is :code:`"1"`.
The following figure shows the values of lift coefficient, :code:`CL`, and angle of attack, :code:`alphaAngle`, versus the steps, :code:`pseudoStep`. As specified in the configuration file, the initial value of angle of attack is :code:`3.06`. With the specified user-defined dynamics, :code:`alphaAngle` is adjusted to a final value of :code:`4.61` so that :code:`CL = 0.4` is achieved.
.. _fig1_alphaController:
.. figure:: Figures/alphaControllerResults.png
:scale: 35%
:align: center
*Results of alphaController for om6Wing case study.*
.. _UDDGridRotation:
Dynamic Grid Rotation using Structural Aerodynamic Load
--------------------------------------------------------------------------
In this example, we present an illustrative case where a flat plate rotates while being subjected to a rotational spring moment and damper as well as aerodynamic loads.
As you can see below, the flate plate mesh consists of two zones separated by a circular sliding interface. The inner zone will rotate while the outer zone remains stationary. This is the same mesh topology as what is used for :ref:`rotating propeller simulations`.
.. _fig_flatplate_mesh:
.. figure:: Figures/flat_plate_UDD_tutorial_mesh.png
:width: 600
:align: center
*Mesh of flat plate showing sliding interface for rotation.*
Below is a video showing the flutter motion of the plate.
.. _fig2_flutterMotion:
.. figure:: Figures/plateFlutterMotion.gif
:width: 600
:align: center
*Flutter motion of the plate under aerodynamic and structural forces.*
The :code:`userDefinedDynamics` and :code:`slidingInterfaces` section of the case JSON is presented below:
.. literalinclude:: Files/ASI.json
:linenos:
For this case the dynamics of the plate are:
.. math::
:label: def_dynamics-2
\dot{\omega} = \left[ \tau_y - K \cdot \left( \theta-\theta_0 \right) - 2\zeta\omega_N I \omega \right]/I
\omega -\omega_0= \int_{0}^{t} \dot{\omega}\left(t'\right) \, dt'
\theta -\theta_0 = \int_{0}^{t} \omega\left(t'\right) \, dt'
which corresponds to the three :code:`updateLaw` arguments in the above JSON input file. The symbols are defined in the table below.
.. table::
:widths: 25 50
+---------------------+-------------------------------------------------------------------------------+
| Symbol | Description |
+=====================+===============================================================================+
| :math:`\theta` | Rotation angle of the plate in radians. |
+---------------------+-------------------------------------------------------------------------------+
| :math:`\omega` | Rotation velocity of the plate in radians. |
+---------------------+-------------------------------------------------------------------------------+
| :math:`\dot{\omega}`| Rotation acceleration of the plate in radians. |
+---------------------+-------------------------------------------------------------------------------+
| :math:`\tau_y` | Aerodynamic moment exerted on the plate. |
+---------------------+-------------------------------------------------------------------------------+
| :math:`K` | Stiffness of the spring attached to the plate at the structural support. |
+---------------------+-------------------------------------------------------------------------------+
| :math:`\zeta` | Structural damping ratio. |
+---------------------+-------------------------------------------------------------------------------+
| :math:`\omega_N` | Structural natural angular frequency. |
+---------------------+-------------------------------------------------------------------------------+
| :math:`I` | Structural moment of inertia. |
+---------------------+-------------------------------------------------------------------------------+
In the :code:`inputVars`, the :code:`momentY` is the y-component of the rotational moment imposed by the aerodynamic loading. The moment is calculated relative to the :ref:`momentCenter ` on :code:`"plateBlock/noSlipWall"`.
The name of rotating volume zone that we want to control is :code:`plateBlock` and thus the :code:`"outputTargetName"` is set to :code:`plateBlock`. The first, second and third state variable correspond to :math:`\dot{\omega}`, :math:`\omega`, and :math:`\theta`, respectively.
.. UDE_cs_torque:
Custom set of user defined forces and moments
---------------------------------------------------
In this example we will use a simple airplane geometry with control surfaces (aileron and rudders) to showcase how a user could calculate user defined forces and moments to track the evolution of the torque on the control surfaces as the simulation converges.
To generate the geometry and associated mesh we will simply use these `input csm `_ , mesh `surface definition `_ , `volume definition `_. See :ref:`this example tutorial` or :ref:`this one` for more information on how to do that.
As you can see in :ref:`this image below`, the airplane has 2 ailerons and 2 rudders. All 4 control surfaces move in a simple rotation around an axis. For each of these we will monitor the torque around the rotation axis.
.. note::
Due to the coarseness of the mesh settings selected, the solver may give you a warning when you go to submit the case with that configuration. You can ignore it and proceed to the case submission.
.. _udd_airplane_isometric:
.. figure:: Figures/UDD_FM_monitor_airplane_isometric.png
:align: center
:scale: 60%
Isometric view of a simple airplane with 4 control surfaces
Now that we have a volume mesh ready we need to define the run parameters. In this tutorial we will focus on the definition of :code:`userDefinedDynamics`. Please download this `sample input json file `_. The setup is for our simple airplane flying at 50m/s at 10 |deg| angle of attack.
The control surfaces have their rotation centers and axes as:
.. _cs_coord_table:
.. csv-table::
:file: Tables/cs_coordinates.csv
:align: center
:header-rows: 1
:delim: @
.. _hinge_moment_ref_check:
.. note::
In this tutorial we have arbitrarily set the :code:`geometry/momentCenter` at the same X and Z location as the Aileron rotation axes ( as in, on the Ailerons' rotation axes). Since the aileron hinges axes is aligned with the Y axis, this allows us to check the hinge moment values we get for the ailerons against the Moment Y value for that same aileron automatically generated by the software. They are the same.
The interesting part is the definition of the control surface hinge moments via user defined variables and the automatic tracking of the outputs to get the evolution of those forces as the run converges to see if we have reached satisfactory convergence of the desired forces and moments.
Here, we will convert the monitored torques into SI forces (N) and moments (Nm).
The input json file has an optional :code:`userDefinedDynamics` (UDD) section which contains 4 similar items. Each control surface requires its own :code:`userDefinedDynamics` instance with the only differences being the names and geometric inputs for each control surface as per the :ref:`table above`.
If we look at the "userDefinedDynamics" section of the configuration json file we see that the UDD instance has 4 items each containing 6 sections. The 4 items are the 4 control surfaces. The 6 sections are described below.
.. _udd_json_section:
.. figure:: Figures/udd_json_screenshot.png
:align: center
UDD section of the case configuration json file.
.. _udd_constants:
constants
~~~~~~~~~~
Named constants that will be used in the update law. "newCenterX,Y,Z" and "newAxisX,Y,Z" are the coordinates of the new moment reference center and rotation axes we want to project onto. They define the aileron or rudder's hingelines.
.. literalinclude:: Files/constants.json
:linenos:
updateLaw
~~~~~~~~~~
This block is a list of expressions used to define the outputs we want to save. Because we convert the raw solver outputs to dimensional metric value we need to calculate the :math:`\rho_\infty C_\infty^2 L_\text{gridUnit}^3` :ref:` moment dimensionalization values`.
The equation on lines 2 defines that dimensionalization value. Lines 3 to 5 define the new dimensional moments as per the new moment reference center and rotation axis we defined :ref:`above`.
The last line defines the total hinge moment, as in the torque that the control surface's actuating mechanism will have to overcome to rotate the control surface.
.. literalinclude:: Files/update_law.json
:linenos:
inputBoundaryPatches
~~~~~~~~~~~~~~~~~~~~~
Patches of the mesh that we want to apply the calculations to. For example:
.. literalinclude:: Files/input_patches.json
:linenos:
visualizing the results
~~~~~~~~~~~~~~~~~~~~~~~~
All the UDD outputs can be downloaded in the :code:`results/udd_DYNAMICSNAME_v2.csv` where :code:`DYNAMICSNAME` matches the values entered in the "dynamicsName" section.
They can also be automatically visualized in the :code:`Monitors` tab for that case on our website.
Remember that the hinge torque's dimensional value is calculated in :code:`state[4]`. By clicking on each line's legend icon you can disable the values you are not interested in, leaving only the desired values' convergence history displayed.
.. _udd_r_ail_torq_hist:
.. figure:: Figures/r_ail_torq_hist.png
:align: center
Dimensional hinge torque convergence history for the right Aileron.
A good way to check that the calculations are correct is to compare the overall forces on a given component with the resulting torque values. If you :ref:`recall above` it was mentioned that we have purposefully put the simulation :code:`geometry/momentCenter` on the aileron hinge line. Both are at X=5.7542 and since the aileron hinge vector [0,1,0] is aligned with the Y axis, the left and right aileron My magnitudes are identical to the left and right aileron hinge torque. The left aileron has it hinge vector set to [0,-1,0] hence the sign will be flipped between My and hinge torque.
To check that the calculations are correct let's get the final value of the right aileron hinge torque. To do that please click on the DataView icon circled in red in the above figure. Scroll to the :code:`state[4]` column and down to the last row. Make sure to navigate to the table's final section(by clicking on the 6 below the table as circled in red below) That highlighted value of -52.65Nm is the right aileron hinge torque in dimensional Newton-meters.
.. _udd_r_ail_torq_hist_table:
.. figure:: Figures/r_ail_torq_hist_table.png
:align: center
Dimensional hinge torque table for the right Aileron.
Now by going back to the :code:`Forces` tab, down to the :code:`Forces and Moments by surface` section, we can similarly find the final CMy= :math:`-5.73 \times 10^{-4}` for the right aileron:
.. _udd_r_ail_cmy_table:
.. figure:: Figures/r_ail_cmy_table.png
:align: center
CMy table for the right Aileron.
Now if we :ref:`dimensionalize that value `, we get
.. math::
:label: results_validation
\begin{align}
-5.73 \times 10^{-4} &\cdot \frac{1}{2} \rho_\infty U_\text{ref}^2 A_\text{ref} L_\text{ref}\left[1\right] \\
&= -5.73\times 10^{-4} \cdot \frac{1}{2} \cdot 1.225\frac{kg}{m^3} \cdot (50\frac{m}{s})^2 \cdot 60m^2 \cdot 1m \\
&= -52.65N \cdot m
\end{align}