diff --git a/doc/XMLreference.rst b/doc/XMLreference.rst index aff1d69a..b605c1dc 100644 --- a/doc/XMLreference.rst +++ b/doc/XMLreference.rst @@ -2838,31 +2838,48 @@ helps clarify the role of bodies and geoms in MuJoCo. .. _body-geom-fluidshape: :at:`fluidshape`: :at-val:`[none, ellipsoid], "none"` - "ellipsoid" Activates geom-level stateless fluid interaction model based on an ellipsoidal approximation of the geom - shape. When active, the model based on :ref:`body inertia sizes ` is disabled for the parent body. + "ellipsoid" activates the geom-level fluid interaction model based on an ellipsoidal approximation of the geom + shape. When active, the model based on :ref:`body inertia sizes ` is disabled for the parent body. + See section on :ref:`ellipsoid-based` fluid interaction model for details. .. _body-geom-fluidcoef: :at:`fluidcoef`: :at-val:`real(5), "0.5 0.25 1.5 1.0 1.0"` Dimensionless coefficients of fluid interaction model, as follows. + See section on :ref:`ellipsoid-based` fluid interaction model for details. -.. table:: +.. list-table:: + :width: 60% :align: left + :widths: 1 5 2 1 + :header-rows: 1 + + * - Index + - Description + - Symbol + - Default + * - 0 + - Blunt drag coefficient + - :math:`C_{D, \text{blunt}}` + - 0.5 + * - 1 + - Slender drag coeficient + - :math:`C_{D, \text{slender}}` + - 0.25 + * - 2 + - Angular drag coeficient + - :math:`C_{D, \text{angular}}` + - 1.5 + * - 3 + - Kutta lift coeficient + - :math:`C_K` + - 1.0 + * - 4 + - Magnus lift coeficient + - :math:`C_M` + - 1.0 - +--------+-----------------------------+----------+ - | Index | Description | Default | - +========+=============================+==========+ - | 0 | Blunt drag coefficient. | 0.5 | - +--------+-----------------------------+----------+ - | 1 | Slender drag coeficient. | 0.25 | - +--------+-----------------------------+----------+ - | 2 | Angular drag coefficient. | 1.5 | - +--------+-----------------------------+----------+ - | 3 | Kutta lift coeficient. | 1.0 | - +--------+-----------------------------+----------+ - | 4 | Magnus lift coeficient. | 1.0 | - +--------+-----------------------------+----------+ .. _body-geom-user: diff --git a/doc/computation/fluid.rst b/doc/computation/fluid.rst new file mode 100644 index 00000000..c5153334 --- /dev/null +++ b/doc/computation/fluid.rst @@ -0,0 +1,436 @@ +Fluid forces +============ + +Proper simulation of fluid dynamics is beyond the scope of MuJoCo, and would be too slow for the applications we aim to +facilitate. Nevertheless we provide two phenomenological models which are sufficient for simulating behaviors +such as flying and swimming. These models are *stateless*, in the sense that no additional states are assigned to the +surrounding fluid, yet are able to capture the salient features of rigid bodies moving through a fluid medium. + +Both models are enabled by setting the :ref:`density` and :ref:`viscosity` attributes +to positive values. These parameters correspond to the density :math:`\rho` and viscosity :math:`\beta` of the medium. + +1. The :ref:`Inertia-based model`, uses only viscosity and density, inferring geometry from body + equivalent-inertia boxes. +2. The :ref:`Ellipsoid-based model ` is more elaborate, using an ellipsoid approximation of geoms. + In addition to the global viscosity and density of the medium, this model exposes 5 tunable parameters per + interacting geom. + +.. tip:: + As detailed in the :ref:`Numerical Integration` section, implicit integration significantly improves + simulation stability in the presence of velocity-dependent forces. Both of the fluid-force models described below + exhibit this property, so the ``implicit`` or ``implicitfast`` :ref:`intergrators` are + recommended when using fluid forces. The required analytic derivatives for both models are fully implemented. + +.. _flInertia: + +Inertia model +------------- + +In this model, the shape of each body, for fluid dynamics purposes, is assumed to be the *equivalent inertia box*, +which can also be visualized. Each forward-facing (relative to the linear velocity) face of the box experiences force +along its normal direction. All faces also experience torque due to the angular velocity; this torque is obtained by +integrating the force resulting from the rotation over the surface area. In this sub-section, let :math:`v` and +:math:`\omega` denote the linear and angular body velocity in the body local frame (aligned with the equivalent +inertia box), and :math:`s` the 3D vector of box sizes. When the contributions from all faces are added, the resulting +force and torque applied to the body by a fluid of density :math:`\rho`, in local body coordinates, have the +:math:`i`-th component + +.. math:: + \begin{aligned} + \text{density force} : \quad &- {1 \over 2} \rho s_j s_k |v_i| v_i \\ + \text{density torque} : \quad &- {1 \over 64} \rho s_i \left(s_j^4 + s_k^4 \right) |\omega_i| \omega_i \\ + \end{aligned} + +This model implicitly assumes high Reynolds numbers, with lift-to-drag ratio equal to the tangent of the angle of +attack. One can also specify a non-zero :ref:`wind`, which is a 3D vector subtracted from the body linear +velocity in the fluid dynamics computation. + +Each body also experiences a force and a torque proportional to the viscosity :math:`\beta` and opposite to its linear and +angular velocity. Note that viscosity can be used independent of density, to make the simulation more damped. We use the +formulas for a sphere at low Reynolds numbers, with diameter :math:`d` equal to the average of the equivalent inertia +box sizes. The resulting 3D force and torque in local body coordinates are + +.. math:: + \begin{aligned} + \text{viscosity force} : \quad &- 3 \beta \pi d v \\ + \text{viscosity torque} : \quad &- \beta \pi d^3 \omega \\ + \end{aligned} + +.. _flEllipsoid: + +Ellipsoid model +--------------- + +.. cssclass:: caption-small +.. figure:: ../images/computation/fruitfly.png + :figwidth: 50% + :align: right + + The flight-capable Drosophila Melanogaster model in this figure will be described in a + forthcoming publication. + + +In this section we describe and derive a stateless model of the forces exerted onto a moving rigid body by the +surrounding fluid, based on an ellipsoidal approximation of geom shape. This model provides finer-grained control of the +different types of fluid forces than the inertia-based model of the previous section. The motivating use-case for this +model is insect flight, see figure on the right. + + +Summary +~~~~~~~ + +The model is activated per-geom by setting the :ref:`fluidshape` attribute to ``ellipsoid``, which +also disables the inertia-based model for the parent body. The +5 numbers in the :ref:`fluidcoef` attribute correspond to the following semantics + + +.. list-table:: + :width: 60% + :align: left + :widths: 1 5 2 1 + :header-rows: 1 + + * - Index + - Description + - Symbol + - Default + * - 0 + - Blunt drag coefficient + - :math:`C_{D, \text{blunt}}` + - 0.5 + * - 1 + - Slender drag coeficient + - :math:`C_{D, \text{slender}}` + - 0.25 + * - 2 + - Angular drag coeficient + - :math:`C_{D, \text{angular}}` + - 1.5 + * - 3 + - Kutta lift coeficient + - :math:`C_K` + - 1.0 + * - 4 + - Magnus lift coeficient + - :math:`C_M` + - 1.0 + +Elements of the model are a generalization of :cite:t:`andersen2005b` to 3 dimensions. +The force :math:`\mathbf{f}_{\text{fluid}\rightarrow \text{solid}}` and torque +:math:`\mathbf{g}_{\text{fluid} \rightarrow \text{solid}}` exerted by the fluid onto the solid are +the sum of of the terms + +.. math:: + \begin{align*} + \mathbf{f}_{\text{fluid} \rightarrow \text{solid}} &= \mathbf{f}_A + \mathbf{f}_D + \mathbf{f}_M + \mathbf{f}_K \\ + \mathbf{g}_{\text{fluid} \rightarrow \text{solid}} &= \mathbf{g}_A + \mathbf{g}_D + \end{align*} + +Where subscripts :math:`A`, :math:`D`, :math:`M` and :math:`K`, denote Added mass, viscous Drag, Magnus lift and +Kutta lift, respectively. The :math:`D`, :math:`M` and :math:`K` terms are scaled by the respective +:math:`C_D`, :math:`C_M` and :math:`C_K` coefficients above, while the added mass term cannot be scaled. + +Notation +~~~~~~~~ + +We describe the motion of the object in an inviscid, incompressible quiescent fluid of density :math:`\rho`. The +arbitrarily-shaped object is described in the model as the equivalent ellipsoid of semi-axes +:math:`\mathbf{d} = \{d_x, d_y, d_z\}`. +The problem is described in a reference frame aligned with the sides of the ellipsoid and moving with it. The +body has velocity :math:`\mathbf{v} = \{v_x, v_y, v_z\}` and angular velocity +:math:`\boldsymbol{\omega} = \{\omega_x, \omega_y, \omega_z\}`. We will also use + +.. math:: + \begin{align*} + d_\text{max} &= \max(d_x, d_y, d_z) \\ + d_\text{min} &= \min(d_x, d_y, d_z) \\ + d_\text{mid} &= d_x + d_y + d_z - d_\text{max} - d_\text{min} + \end{align*} + +The Reynolds number is the ratio between inertial and viscous forces within a flow and is defined as :math:`Re=u~l/\beta`, where +:math:`\beta` is the kinematic viscosity of the fluid, :math:`u` is the characteristic speed of the flow (or, by change of frame, the +speed of the body), and :math:`l` is a characteristic size of the flow or the body. + +We will use :math:`\Gamma` to denote circulation, which is the line integral of the velocity field around a closed curve +:math:`\Gamma = \oint \mathbf{v} \cdot \textrm{d} \mathbf{l}` and, due to Stokes' Theorem, +:math:`\Gamma = \int_S \nabla \times \mathbf{v} \cdot \textrm{d}\mathbf{s}`. +In fluid dynamics notation the symbol :math:`\boldsymbol{\omega}` is often used for the +vorticity, defined as :math:`\nabla \times \mathbf{v}`, rather than the angular velocity. For a rigid-body motion, the +vorticity is twice the angular velocity. + +Finally, we use the subscripts :math:`i, j, k` to denote triplets of equations that apply symmetrically to +:math:`x, y, z`. For example :math:`a_i = b_j + b_k` is shorthand for the 3 equations + +.. math:: + \begin{align*} + a_x &= b_y + b_z \\ + a_y &= b_x + b_z \\ + a_z &= b_x + b_y + \end{align*} + +.. _flProjection: + +Ellipsoid projection +~~~~~~~~~~~~~~~~~~~~ + +We provide the following result without proof. For the derivation, contact the development team. + +.. admonition:: Lemma + :class: note + + Given an ellipsoid at the origin with semi-axes :math:`(d_x, d_y, d_z)` aligned + with the coordinate axes :math:`(x, y, z)`, and a unit vector :math:`\mathbf{u} = (u_x, u_y, u_z)`, + the area projected by the ellipsoid onto the plane normal to :math:`\mathbf{u}` is + + .. math:: + A^{\mathrm{proj}}_{\mathbf{u}} = \pi \sqrt{\frac{d_y^4 d_z^4 u_x^2 + d_x^4 d_z^4 u_y^2 + d_x^4 d_y^4 u_z^2}{d_y^2 d_z^2 u_x^2 + d_x^2 d_z^2 u_y^2 + d_x^2 d_y^2 u_z^2}} + + +Added mass +~~~~~~~~~~ + +For a body moving in a fluid, added mass or virtual mass measures the inertia of the fluid that is moved due to the +body's motion. It can be derived from potential flow theory (i.e. it is present also for inviscid flows). + +Following Chapter 5 of :cite:t:`lamb1932`, the forces :math:`\mathbf{f}_{V}` and torques :math:`\mathbf{g}_{V}` exerted +onto a moving body due to generation of motion in the fluid from rest can be written as: + +.. math:: + \begin{align*} + \mathbf{f}_{A} &= - \frac{\textrm{d}}{\textrm{d} t} \nabla_{\mathbf{v}} \mathcal{T} + \nabla_{\mathbf{v}} \mathcal{T} \times \boldsymbol{\omega} \\ + \mathbf{g}_{A} &= - \frac{\textrm{d}}{\textrm{d} t} \nabla_{\boldsymbol{\omega}} \mathcal{T} + \nabla_{\mathbf{v}} \mathcal{T} \times \mathbf{v} + \boldsymbol{\omega} \times \nabla_{\boldsymbol{\omega}} \mathcal{T} + \end{align*} + +where :math:`\mathcal{T}` is the kinetic energy of the fluid alone. These forces are often described as added or +virtual mass because they are due to the inertia of the fluid that is to moved or deflected by the accelerating body. In +fact, for a body with constant linear velocity these forces reduce to zero. We consider the body as having three planes +of symmetry because under this assumption the kinetic energy greatly simplifies and can be written as: + +.. math:: + 2 \mathcal{T} = m_{A, x} v_x^2 + m_{A, y} v_y^2 + m_{A, z} v_z^2 + + I_{A, x} \omega_x^2 + I_ {A, y} \omega_y^2 + I_{A, y} \omega_z^2 + + +For convenience we introduce the added-mass vector :math:`\mathbf{m}_A = \{m_{A, x}, m_{A, y}, m_{A, z}\}` and added-moment of +inertia vector :math:`\mathbf{I}_A = \{I_{A, x}, I_{A, y}, I_{A, z}\}`. Each of these quantities should estimate the inertia +of the moved fluid due the motion of the body in the corresponding direction and can be derived from potential flow +theory for some simple geometries. + +For a body with three planes of symmetry, we can write in compact form the forces and torques due to added inertia: + +.. math:: + \begin{align*} + \mathbf{f}_{A} &= - \mathbf{m}_A \circ \dot{\mathbf{v}} + \left(\mathbf{m}_A \circ \mathbf{v} \right) \times \boldsymbol{\omega} \\ + \mathbf{g}_{A} &= - \mathbf{I}_A \circ \dot{\boldsymbol{\omega}} + \left(\mathbf{m}_A \circ \mathbf{v} \right) \times \mathbf{v} + \left(\mathbf{I}_A \circ \boldsymbol{\omega} \right) \times \boldsymbol{\omega} + \end{align*} + +Here :math:`\circ` denotes an element-wise product, :math:`\dot{\mathbf{v}}` is the linear acceleration and +:math:`\dot{\boldsymbol{\omega}}` is the angular acceleration. :math:`\mathbf{m}_A \circ \mathbf{v}` and +:math:`\mathbf{I}_A \circ \boldsymbol{\omega}` are the virtual linear and angular momentum respectively. + +For an ellipsoid of semi-axis :math:`\mathbf{d} = \{d_x, d_y, d_z\}` and volume :math:`V = 4 \pi d_x d_y d_z / 3`, the +virtual inertia coefficients were derived by :cite:t:`tuckerman1925`. Let: + +.. math:: + \kappa_i = \int_0^\infty \frac{d_i d_j d_k}{\sqrt{(d_i^2 + \lambda)^3 (d_j^2 + \lambda) (d_k^2 + \lambda)}} \textrm{d} \lambda + + +It should be noted that these coefficients are non-dimensional (i.e. if all semi-axes are multiplied by the same scalar +the coefficients remain the same). The virtual masses of the ellipsoid are: + +.. math:: + m_{A, i} = \rho V \frac{\kappa_i}{2 - \kappa_i} + +And the virtual moments of inertia are: + +.. math:: + I_{A, i} = \frac{\rho V}{5} \frac{(d_j^2 - d_k^2)^2 (\kappa_k-\kappa_j)}{2(d_j^2 - d_k^2) + (d_j^2 + d_k^2) (\kappa_j-\kappa_k)} + +Viscous drag +~~~~~~~~~~~~ + +The drag force acts to oppose the motion of the body relative to the surrounding flow. We found that viscous forces +serve also to reduce the stiffness of the equations of motion extended with the fluid dynamic terms. For this reason, we +opted to err on the conservative side and chose approximations of the viscous terms that may overestimate dissipation. + +Despite being ultimately caused by viscous dissipation, for high Reynolds numbers the drag is independent of the +viscosity and scales with the second power of the velocity. It can be written as: + +.. math:: + \begin{align*} + \mathbf{f}_\text{D} = - C_D~\rho~ A_D ~ \|\mathbf{v}\|~ \mathbf{v}\\ + \mathbf{g}_\text{D} = - C_D \rho~ I_D ~ \|\boldsymbol{\omega}\| ~ \boldsymbol{\omega} + \end{align*} + +Where :math:`C_D` is a drag coefficient, and :math:`A_D` is a reference surface area (e.g. a measure of the projected +area on the plane normal to the flow), and :math:`I_D` a reference moment of inertia. + +Even for simple shapes, the terms :math:`C_D`, :math:`A_D` and :math:`I_D` need to be tuned to the problem-specific +physics and dynamical scales :cite:p:`duan2015`. For example, the drag coefficient :math:`C_D` generally decreases with +increasing Reynolds numbers, and a single reference area :math:`A_D` may not be sufficient to account for the skin +drag for highly irregular or slender bodies. For example, experimental fits are derived from problems ranging from +falling playing cards :cite:p:`wang2004,andersen2005a,andersen2005b` to particle transport :cite:p:`loth2008, +bagheri2016`. + +We derive a formula for :math:`\mathbf{f}_\text{D}` based on two surfaces :math:`A^\text{proj}_\mathbf{v}` and +:math:`A_\text{max}`. The first, :math:`A^\text{proj}_\mathbf{v}`, is the cylindrical projection of the body onto a +plane normal to the velocity :math:`\mathbf{v}`. The second is the maximum projected surface +:math:`A_\text{max} = 4 \pi d_{max} d_{min}`. + +.. math:: + \mathbf{f}_\text{D} = - \rho~ \big[ C_{D, \text{blunt}} ~ A^\text{proj}_\mathbf{v} ~ + + C_{D, \text{slender}}\left(A_\text{max} - A^\text{proj}_\mathbf{v} \right) \big] ~ \|\mathbf{v}\|~ \mathbf{v} + +The formula and derivation for :math:`A^\text{proj}_\mathbf{v}` is given in the :ref:`lemma` above. + +We propose an analogous model for the angular drag. For each Cartesian axis we consider the moment of inertia of the +maximum swept ellipsoid obtained by the rotation of the body around the axis. The resulting diagonal entries of the +moment of inertia are: + +.. math:: + \mathbf{I}_{D,ii} = \frac{8\pi}{15} ~d_i ~\max(d_j, ~d_k)^4 . + +Given this reference moment of inertia, the angular drag torque is computed as: + +.. math:: + \mathbf{g}_\text{D} = - \rho ~ \boldsymbol{\omega} ~ \Big( \big[ C_{D, \text{angular}} ~ \mathbf{I}_D ~ + + C_{D, \text{slender}} \left(\mathbf{I}_\text{max} - \mathbf{I}_D \right) \big] \cdot \boldsymbol{\omega} \Big) + + +Here :math:`\mathbf{I}_\text{max}` is a vector with each entry equal to the maximal component of :math:`\mathbf{I}_D`. + +The viscosity :math:`\beta` +For Reynolds numbers around or below :math:`O(10)`, the drag is best approximated as linear in the flow velocity +(e.g. Stokes' law). For example, for a sphere the drag force :cite:p:`stokes1850` and torque :cite:p:`lamb1932` are: + +.. math:: + \begin{align*} + \mathbf{f}_\text{S} &= - 6 \pi r_D \rho ~ \beta \mathbf{v}\\ + \mathbf{g}_\text{S} &= - 8 \pi r_D^3 \rho ~ \beta \boldsymbol{\omega} + \end{align*} + +Here, :math:`r_D` is the radius of the sphere and :math:`\beta` is the kinematic viscosity of the medium (e.g. +:math:`1.48~\times 10^{-5}~m^2/s` for ambient-temperature air and :math:`0.89 \times 10^{-4}~m^2/s` for water). Here, +for simplicity, we estimate the radius of the equivalent sphere as :math:`r_D = (d_x + d_y + d_z)/3`. To make a +quantitative example, Stokes' law become accurate for room-temperature air if +:math:`u\cdot l \lesssim 2 \times 10^{-4}~m^2/s`, where :math:`u` is the speed and :math:`l` a characteristic length of +the body. + +Viscous lift +~~~~~~~~~~~~ + +The Kutta-Joukowski theorem calculates the lift :math:`L` of a two-dimensional body translating in a uniform flow with +speed `u` as :math:`L = \rho u \Gamma`. Here, :math:`\Gamma` is the circulation around the body. In the next +subsections we define two sources of circulation and the resulting lift forces. + +Magnus force +^^^^^^^^^^^^ + +.. cssclass:: caption-small +.. figure:: ../images/computation/magnus.png + :figwidth: 45% + :align: right + + Smoke flow visualization of the flow past a rotating cylinder (WikiMedia Commons, CC BY-SA 4.0). Due to viscosity, + the rotating cylinder deflects the incoming flow upward and receives a downwards force (red arrow). + +The Magnus effect describes the motion of a rotating object moving through a fluid. Through viscous effects, a spinning +object induces rotation in the surrounding fluid. This rotation deflects the trajectory of the fluid past the object +(i.e. it causes linear acceleration), and the object receives an equal an opposite reaction. For a cylinder, the Magnus +force per unit length of the cylinder can be computed as :math:`F_\text{M} / L = \rho v \Gamma`, where :math:`\Gamma` +is the circulation of the flow caused by the rotation and :math:`v` the velocity of the object. We estimate this force +for an arbitrary body as: + +.. math:: + \mathbf{f}_{\text{M}} = C_M ~\rho~ V~ \boldsymbol{\omega}\times\mathbf{v} , + +where :math:`V` is the volume of the body and :math:`C_M` is a coefficient for the force, typically set to 1. + +It's worth making an example. To reduce the number of variables, suppose a body rotating in only one direction, e.g. +:math:`\boldsymbol{\omega} = \{0, 0, \omega_z\}`, translating along the other two, e.g. :math:`\mathbf{v} = \{v_x, v_y, 0\}`. The +sum of the force due to added mass and the force due to the Magnus effect along, for example, :math:`x` is: + +.. math:: + \frac{f}{\pi \rho d_z} = v_y \omega_z \left(2 d_x \min\{d_x, d_z\} - (d_x + d_z)^2\right) + +Note that the two terms have opposite signs. + +Kutta condition +^^^^^^^^^^^^^^^ + +A stagnation point is a location in the flow field where the velocity is zero. For a body moving in a flow (in 2D, in +the frame moving with the body) there are two stagnation points: in the front, where the stream-lines separate to either +sides of the body, and in the rear, where they reconnect. A moving body with a sharp trailing (rear) edge will generate +in the surrounding flow a circulation of sufficient strength to hold the rear stagnation point at the trailing edge. +This is the Kutta condition, a fluid dynamic phenomenon that can be observed for solid bodies with sharp corners, such +as slender bodies or the trailing edges of airfoils. + +.. cssclass:: caption-small +.. figure:: ../images/computation/kutta_cond_plate.svg + :figwidth: 95% + :align: left + + Sketch of the Kutta condition. Blue lines are streamlines and the two magenta points are the stagnation points. The + dividing streamline, which connects the two stagnation points, is marked in green. The dividing streamline and the + body inscribe an area where the flow is said to be "separated" and recirculates within. This circulation produces an + upward force acting on the plate. + +For a two-dimensional flow sketched in the figure above, the circulation due to the Kutta condition can be estimated as: +:math:`\Gamma_\text{K} = C_K ~ d_x ~ \| \mathbf{v}\| ~ \sin(2\alpha)`, +where :math:`C_K` is a lift coefficient, and :math:`\alpha` is the angle between the velocity vector and its projection +onto the surface. The lift force per unit length can be computed with the Kutta–Joukowski theorem as +:math:`\mathbf{f}_K / L = \rho \Gamma_\text{K} \times \mathbf{v}`. + +In order to extend the lift force equation to three-dimensional motions, we consider the normal +:math:`\mathbf{n}_{s, \mathbf{v}} = \{\frac{d_y d_z}{d_x}v_x, \frac{d_z d_x}{d_y}v_y, \frac{d_x d_x}{d_z}v_z\}` +to the cross-section of the body which generates the body's projection :math:`A^\text{proj}_\mathbf{v}` onto a plane +normal to the velocity given in the :ref:`lemma` above and the corresponding unit vector +:math:`\hat{\mathbf{n}}_{s, \mathbf{v}}`. +We use this direction to decompose :math:`\mathbf{v} = \mathbf{v}_\parallel ~+~ \mathbf{v}_\perp` with +:math:`\mathbf{v}_\perp = \left(\mathbf{v} \cdot \hat{\mathbf{n}}_{s, \mathbf{v}}\right) \hat{\mathbf{n}}_{s, \mathbf{v}}`. +We write the lift force as: + +.. math:: + \begin{align*} + \mathbf{f}_\text{K} &= \frac{C_K~\rho~ A^\text{proj}_\mathbf{v}}{\|\mathbf{v}\|} + \left( \mathbf{v} \times \mathbf{v}_\parallel\right)\times \mathbf{v} \\ + &= C_K~\rho~ A^\text{proj}_\mathbf{v} \left(\hat{\mathbf{v}} \cdot \hat{\mathbf{n}}_{s, \mathbf{v}}\right) + \left( \hat{\mathbf{n}}_{s, \mathbf{v}} \times \mathbf{v} \right)\times \mathbf{v} + \end{align*} + +Here, :math:`\hat{\mathbf{v}}` is the unit-normal along :math:`\mathbf{v}`. Note that the direction of :math:`\hat{\mathbf{n}}_{s, +\mathbf{v}}` differs from :math:`\hat{\mathbf{v}}` only on the planes where the semi-axes of the body are unequal. So for +example, for spherical bodies :math:`\hat{\mathbf{n}}_{s, \mathbf{v}} \equiv \hat{\mathbf{v}}` and by construction +:math:`\mathbf{f}_\text{K} = 0`. + +Let's unpack the relation with an example. Suppose a body with :math:`d_x = d_y` and :math:`d_z \ll d_x`. Note that the vector +:math:`\hat{\mathbf{n}}_{s, \mathbf{v}} \times \hat{\mathbf{v}}` gives the direction of the circulation induced by the +deflection of the flow by the solid body. Along :math:`z`, the circulation will be proportional to :math:`\frac{d_y d_z}{d_x}v_x v_y +- \frac{d_z d_x}{d_y}v_x v_y = 0` (due to :math:`d_x = d_y`). Therefore, on the plane where the solid is blunt, the motion +produces no circulation. + +Now, for simplicity, let :math:`v_x = 0`. In this case also the circulation along :math:`y`, proportional +to :math:`\frac{d_y d_z}{d_x}v_x v_z - \frac{d_y d_x}{d_y}v_x v_z`, is zero. The only non-zero component of the circulation +will be along :math:`x` and be proportional to :math:`\left(\frac{d_x d_z}{d_y} - \frac{d_x d_y}{d_z}\right) v_y v_z \approx +\frac{d_x^2}{d_z} v_y v_z`. + +We would have :math:`\mathbf{v}_\parallel = \{v_x, 0, v_z\}` and +:math:`\Gamma \propto \{d_z v_y v_z, ~ 0,~ - d_x v_x v_y \} / \|\mathbf{v}\|`. +The motion produces no circulation on the plane where the solid is blunt, and on the other two planes +the circulation is +:math:`\Gamma \propto r_\Gamma ~ \|\mathbf{v}\|~ \sin(2 \alpha) ~ = ~2 r_\Gamma ~\|\mathbf{v}\| ~\sin(\alpha)~\cos(\alpha)` +with :math:`\alpha` the angle between the velocity and its projection on the body on the plane (e.g. on the plane +orthogonal to :math:`x` we have :math:`\sin(\alpha) = v_y/\|\mathbf{v}\|` and +:math:`\cos(\alpha) = v_z/\|\mathbf{v}\|`), and :math:`r_\Gamma`, the lift surface on the plane (e.g. :math:`d_z` for +the plane orthogonal to :math:`x`). Furthermore, the direction of the circulation is given by the cross product (because +the solid boundary "rotates" the incoming flow velocity towards its projection on the body). + +Acknowledgements +~~~~~~~~~~~~~~~~ + +The design and implementation of the model in this section are the work of Guido Novati. + +References +~~~~~~~~~~ + +.. bibliography:: diff --git a/doc/computation.rst b/doc/computation/index.rst similarity index 95% rename from doc/computation.rst rename to doc/computation/index.rst index 51c0cc19..7753555f 100644 --- a/doc/computation.rst +++ b/doc/computation/index.rst @@ -1,6 +1,12 @@ Computation =========== +.. toctree:: + :hidden: + + fluid + + Introduction ------------ @@ -389,59 +395,16 @@ by MuJoCo are also passive in the sense of physics, i.e., they do not increase e callback :ref:`mjcb_passive` and add forces to ``mjData.qfrc_passive`` that may increase energy. This will not interfere with MuJoCo's operation as long as such user forces depend only on position and velocity. -MuJoCo can compute three types of passive forces: spring-dampers in joints and tendons, gravity compensation forces, and -fluid dynamics. +MuJoCo can compute three types of passive forces: -When Euler or the implicit integator are used, joint damping is integrated implicitly which significantly increases -stability. Thus, even though damping can be modeled as an actuator property, it is better to model it as a joint -property. Note also the XML :ref:`joint ` attribute springdamper which automates the creation of mass- -spring-dampers with desired time constants and damping ratios; in that case the compiler computes the stiffness and -damping coefficients of the joint by taking the joint inertia into account. - -Gravity compensation is a force applied to a body's center of mass opposing gravity, see :ref:`body gravcomp` for -details. - -Proper simulation of fluid dynamics is beyond the scope of MuJoCo, and would be too slow for the applications we aim to -facilitate. Nevertheless we provide a phenomenological model which is sufficient for simulating behaviors such as flying -and swimming. It is enabled by setting ``mjModel.opt.viscosity`` and ``mjModel.opt.density`` to positive values (they -are zero by default). These parameters specify the viscosity :math:`\beta` and density :math:`\rho` of the medium and -apply to all bodies. The shape of each body for fluid dynamics purposes is assumed to be the equivalent inertia box, -which can also be visualized. Each forward-facing (relative to the linear velocity) face of the box experiences force -along its normal direction. All faces also experience torque due to the angular velocity; this torque is obtained by -integrating the force resulting from the rotation over the surface area. In this sub-section only, let :math:`v` and -:math:`\omega` denote the linear and angular body velocity in its local frame (which is aligned with the equivalent -inertia box), and :math:`s` the 3D vector of box sizes. When the contributions from all faces are added, the resulting -force and torque in local body coordinates have :math:`i`-th component - -.. math:: - \begin{aligned} - \text{density force} : \quad &- {1 \over 2} \rho s_j s_k |v_i| v_i \\ - \text{density torque} : \quad &- {1 \over 64} \rho s_i \left(s_j^4 + s_k^4 \right) |\omega_i| \omega_i \\ - \end{aligned} - -This model implicitly assumes high Reynolds numbers, with lift-to-drag ratio equal to the tangent of the angle of -attack. One can also specify a non-zero ``mjModel.opt.wind``, which is a 3D vector subtracted from the body linear -velocity in the fluid dynamics computation. - -Each body also experiences a force and a torque proportional to viscosity and opposite to its linear and angular -velocity. Note that viscosity can be used independent of density, to make the simulation more damped. We use the -formulas for a sphere at low Reynolds numbers, with diameter :math:`d` equal to the average of the equivalent inertia -box sizes. The resulting 3D force and torque in local body coordinates are - -.. math:: - \begin{aligned} - \text{viscosity force} : \quad &- 3 \beta \pi d v \\ - \text{viscosity torque} : \quad &- \beta \pi d^3 \omega \\ - \end{aligned} - -Existing phenomenological models of fluid dynamics tend to be valid for one regime (e.g. flat vs. spherical objects) and -it is difficult to construct a model which can be simulated efficiently and is broadly valid. For example, drag forces -are known to transition from being linear in velocity at low Reynolds numbers, to being quadratic in velocity at high -Reynolds numbers. However it is not clear how this transition should occur for lift forces. This motivated the above -separation, using density to compute all quadratic forces and viscosity to compute all linear forces. Users who need a -more detailed simulation of fluid dynamics should leave the density and viscosity parameters set to zero so as to -disable the built-in mechanism. A custom model of fluid dynamics (or any other force field that depends only on position -and velocity) can be implemented in the callback :ref:`mjcb_passive`. +- Spring-dampers in joints and tendons. See the following attribues for details. + |br| **Joints:** + :ref:`stiffness`, :ref:`springref`, + :ref:`damping`, :ref:`springdamper`. + |br| **Tendons:** :ref:`stiffness`, + :ref:`springlength`, :ref:`damping`. +- Gravity compensation forces. See the body :ref:`gravcomp` attribute for details. +- Fluid forces exerted by the surrounding medium. See the :doc:`Fluid forces ` chapter for details. .. _geIntegration: @@ -786,7 +749,6 @@ dimensionality of the constraint residual in each case. computed by the compiler in the initial model configuration ``mjModel.qpos0``. ``distance`` : 1 - .. attention:: Distance equality constraints were removed in MuJoCo version 2.2.2. If you are using an earlier version, please switch to the corresponding version of the documentation. @@ -968,7 +930,7 @@ each component of :math:`f`. The basis vectors are spatial vectors: 3D force fol vectors into the columns of a matrix :math:`E`, the force/torque generated by the constraint force in the contact frame is :math:`E f`. The matrix of basis vectors is constructed as follows. -.. image:: images/computation/contact_frame.svg +.. image:: ../images/computation/contact_frame.svg :width: 700px :align: center @@ -1330,7 +1292,7 @@ representations of the constraint Jacobian and related matrices. the computational complexity of one matrix-vector multiplication (although the constants are larger). It has first-order convergence but nevertheless makes rapid progress in a few iterations. - .. image:: images/computation/gPGS.svg + .. image:: ../images/computation/gPGS.svg :width: 500px :align: center @@ -1435,7 +1397,7 @@ coefficients varying from left to right. Mathematically, the penalty in the pyra the penalty in the elliptic case contains pieces that are quadratics minus square roots of quadratics - allowing circular contours around the tip of the cone. -.. image:: images/computation/softcontact.png +.. image:: ../images/computation/softcontact.png :width: 600px :align: center @@ -1466,18 +1428,25 @@ others can be pruned quickly without a detailed check. MuJoCo has flexible mecha checked in detail. The decision process involves two stages: generation and filtering. Generation - First we generate a list of candidate geom pairs in two ways: + First we generate a list of candidate geom pairs in one of two ways: "pair" or "dynamic". The user can also specify + "all" which merges both sources (and is the default). This is done via the setting ``mjModel.opt.collision``. "Pair" + refers to an explicit list of geom pairs defined with the :ref:`pair ` element in MJCF. It gives the + user full control, however it is a static mechanism (independent of the spatial arrangement of the geoms at runtime) + and can be tedious for large models. It is normally used to supplement the output of the "dynamic" mechanism. Dynamic + generation works with bodies rather than geoms; when a body pair is included this means that all geoms attached to + one body can collide with all geoms attached to the other body. - - Using predefined geom pairs, specified with the :ref:`contact-pair ` element. This mechanism is - independent of the spatial arrangement of the geoms at runtime. - - Generating candidate geom pairs by examining the runtime locations of bodies. First, a sweep-and-prune broadphase - collision stage prunes bodies that cannot collide. Then, a mid-phase collision stage further prunes geom pairs from - remaining body pairs using a static Bounding Volume Hierarchy (a binary tree) of axis-aligned bounding boxes. - Each body is equipped with an AABB tree of its geoms, aligned with the body inertial or geom frames for all inner - or leaf nodes, respectively. - - The user can also explicitly exclude certain body pairs using the :ref:`exclude ` element. + The body pairs are generated via broad-phase collision detection based on a modified sweep-and-prune algorithm. The + modification is that the axis for sorting is chosen as the principal eigenvector of the covariance matrix of all geom + centers - which maximizes the spread. Then, for each body pair, a mid-phase collision detection using a static + bounding volume hierarchy (a BVH binary tree) of axis-aligned bounding boxes (AABB) is performed. Each body is + equipped with an AABB tree of its geoms, aligned with the body inertial or geom frames for all inner or leaf nodes, + respectively. - These two lists are then merged. + Finally, the user can explicitly exclude certain body pairs using the :ref:`exclude ` element + in MJCF. Exclusion is applied when "dynamic" or "all" are selected, but not when "pair" is selected. At the end of + this step we have a list of geoms pairs that is typically much smaller than :math:`n (n-1)/2`, but can still be + pruned further before detailed collision checking. Filtering Next we apply four filters to the list generated in the previous step. Filters 1 and 2 are applied to all geom pairs. diff --git a/doc/conf.py b/doc/conf.py index 0f93efde..0c200d21 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -42,6 +42,7 @@ master_doc = 'index' # extensions coming with Sphinx (named 'sphinx.ext.*') or your custom # ones. extensions = [ + 'sphinxcontrib.bibtex', 'sphinxcontrib.katex', 'sphinxcontrib.youtube', 'sphinx_copybutton', @@ -51,6 +52,9 @@ extensions = [ 'mujoco_include', ] +# Bibtex references for sphinxcontrib.bibtex +bibtex_bibfiles = ['references.bib'] + # Add any paths that contain templates here, relative to this directory. templates_path = ['templates'] diff --git a/doc/css/theme_overrides.css b/doc/css/theme_overrides.css index 8ce81468..44a37261 100644 --- a/doc/css/theme_overrides.css +++ b/doc/css/theme_overrides.css @@ -22,6 +22,10 @@ body { font-size: 93%; } +.caption-small { + font-size: 85%; +} + .sidebar-brand-text { display: none; } diff --git a/doc/images/computation/fruitfly.png b/doc/images/computation/fruitfly.png new file mode 100644 index 00000000..cdb73907 Binary files /dev/null and b/doc/images/computation/fruitfly.png differ diff --git a/doc/images/computation/kutta_cond_plate.svg b/doc/images/computation/kutta_cond_plate.svg new file mode 100644 index 00000000..d130f0bc --- /dev/null +++ b/doc/images/computation/kutta_cond_plate.svg @@ -0,0 +1,56 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +v + + + + +y +x + + + + + +α + + diff --git a/doc/images/computation/magnus.png b/doc/images/computation/magnus.png new file mode 100644 index 00000000..800d47fd Binary files /dev/null and b/doc/images/computation/magnus.png differ diff --git a/doc/index.rst b/doc/index.rst index 192b8f99..e4225fd3 100644 --- a/doc/index.rst +++ b/doc/index.rst @@ -8,7 +8,7 @@ :hidden: overview - computation + computation/index.rst modeling XMLreference programming/index.rst diff --git a/doc/overview.rst b/doc/overview.rst index 0a5bf7e4..a2f980a7 100644 --- a/doc/overview.rst +++ b/doc/overview.rst @@ -48,11 +48,12 @@ Soft, convex and analytically-invertible contact dynamics In the modern approach to contact dynamics, the forces or impulses caused by frictional contacts are usually defined as the solution to a linear or non-linear complementarity problem (LCP or NCP), both of which are NP-hard. MuJoCo is based on a different formulation of the physics of contact which reduces to a convex optimization problem, as - explained in detail in the :doc:`computation` chapter. Our model allows soft contacts and other constraints, and has - a uniquely-defined inverse facilitating data analysis and control applications. There is a choice of optimization - algorithms, including a generalization to the projected Gauss-Seidel method that can handle elliptic friction cones. - The solver provides unified treatment of frictional contacts including torsional and rolling friction, frictionless - contacts, joint and tendon limits, dry friction in joints and tendons, as well as a variety of equality constraints. + explained in detail in the :doc:`computation/index` chapter. Our model allows soft contacts and other constraints, + and has a uniquely-defined inverse facilitating data analysis and control applications. There is a choice of + optimization algorithms, including a generalization to the projected Gauss-Seidel method that can handle elliptic + friction cones. The solver provides unified treatment of frictional contacts including torsional and rolling + friction, frictionless contacts, joint and tendon limits, dry friction in joints and tendons, as well as a variety of + equality constraints. Tendon geometry MuJoCo can model the 3D geometry of tendons - which are minimum-path-length strings obeying wrapping and via-point @@ -746,13 +747,13 @@ multi-threading, but we do so in a way that is different from how object-oriente Softness and slip ~~~~~~~~~~~~~~~~~ -As we will explain at length in the :doc:`computation` chapter, MuJoCo is based on a mathematical model of the physics -of contact and other constraints. This model is inherently soft, in the sense that pushing harder against a constraint -will always result in larger acceleration, and so the inverse dynamics can be uniquely defined. This is desirable -because it yields a convex optimization problem and enables analyses that rely on inverse dynamics, and furthermore, -most contacts that we need to model in practice have some softness. However once we allow soft constraints, we are -effectively creating a new type of dynamics -- namely deformation dynamics -- and now we must specify how these dynamics -behave. This calls for elaborate parameterization of contacts and other constraints, involving the attributes +As we will explain at length in the :doc:`computation/index` chapter, MuJoCo is based on a mathematical model of the +physics of contact and other constraints. This model is inherently soft, in the sense that pushing harder against a +constraint will always result in larger acceleration, and so the inverse dynamics can be uniquely defined. This is +desirable because it yields a convex optimization problem and enables analyses that rely on inverse dynamics, and +furthermore, most contacts that we need to model in practice have some softness. However once we allow soft constraints, +we are effectively creating a new type of dynamics -- namely deformation dynamics -- and now we must specify how these +dynamics behave. This calls for elaborate parameterization of contacts and other constraints, involving the attributes :at:`solref` and :at:`solimp` that can be set per constraints and will be described later. An often confusing aspect of this soft model is that gradual contact slip cannot be avoided. Similarly, frictional diff --git a/doc/programming/simulation.rst b/doc/programming/simulation.rst index 6ff67094..2782923d 100644 --- a/doc/programming/simulation.rst +++ b/doc/programming/simulation.rst @@ -888,8 +888,8 @@ function is essentially free in terms of CPU cost; so do not hesitate to use thi Contacts ~~~~~~~~ -Collision detection and solving for contact forces were explained in detail in the :doc:`../computation` chapter. Here -we further clarify contact processing from a programming perspective. +Collision detection and solving for contact forces were explained in detail in the :doc:`../computation/index` chapter. +Here we further clarify contact processing from a programming perspective. The collision detection stage finds contacts between geoms, and records them in the array ``mjData.contact`` of :ref:`mjContact` data structures. They are sorted such that multiple contacts between the same pair of bodies are diff --git a/doc/references.bib b/doc/references.bib new file mode 100644 index 00000000..e827bad7 --- /dev/null +++ b/doc/references.bib @@ -0,0 +1,90 @@ +@book{lamb1932, + title={Hydrodynamics. Sixth edition.}, + author={Lamb, Horace}, + year={1932}, + publisher={Cambridge University Press} +} + +@book{kennard1967, + title={Irrotational Flow of Frictionless Fluids: Mostly of Invariable Density}, + author={Kennard, Earle Hesse}, + volume={2299}, + year={1967}, + publisher={US Government Printing Office} +} + +@book{tuckerman1925, + title={Inertia factors of ellipsoids for use in airship design}, + author={Tuckerman, LB}, + year={1925}, + publisher={US Government Printing Office} +} + +@article{wang2004, + title={Unsteady forces and flows in low Reynolds number hovering flight: two-dimensional computations vs robotic wing experiments}, + author={Wang, Z Jane and Birch, James M and Dickinson, Michael H}, + journal={Journal of Experimental Biology}, + volume={207}, + number={3}, + pages={449--460}, + year={2004}, + publisher={Company of Biologists} +} +@article{andersen2005a, + title={Unsteady aerodynamics of fluttering and tumbling plates}, + author={Andersen, A and Pesavento, U and Wang, Z Jane}, + journal={Journal of Fluid Mechanics}, + volume={541}, + pages={65--90}, + year={2005}, + publisher={Cambridge University Press} +} +@article{andersen2005b, + title={Analysis of transitions between fluttering, tumbling and steady descent of falling cards}, + author={Andersen, Anders and Pesavento, Umberto and Wang, Z Jane}, + journal={Journal of Fluid Mechanics}, + volume={541}, + pages={91--104}, + year={2005}, + publisher={Cambridge University Press} +} + +@article{stokes1850, + title={On the effect of internal friction of fluids on the motion of pendulums}, + author={Stokes, GG}, + journal={Trans. Camb. phi1. S0c}, + volume={9}, + number={8}, + pages={106}, + year={1850} +} + +@article{loth2008, + title={Drag of non-spherical solid particles of regular and irregular shape}, + author={Loth, E}, + journal={Powder Technology}, + volume={182}, + number={3}, + pages={342--353}, + year={2008}, + publisher={Elsevier} +} +@article{bagheri2016, + title={On the drag of freely falling non-spherical particles}, + author={Bagheri, Gholamhossein and Bonadonna, Costanza}, + journal={Powder Technology}, + volume={301}, + pages={526--544}, + year={2016}, + publisher={Elsevier} +} +@article{duan2015, + title={Sphere drag and heat transfer}, + author={Duan, Zhipeng and He, Boshu and Duan, Yuanyuan}, + journal={Scientific reports}, + volume={5}, + number={1}, + pages={1--7}, + year={2015}, + publisher={Nature Publishing Group} +} diff --git a/doc/requirements.txt b/doc/requirements.txt index 046b0642..ae4efc06 100644 --- a/doc/requirements.txt +++ b/doc/requirements.txt @@ -1,5 +1,6 @@ Sphinx==4.5.0 furo==2022.9.29 +sphinxcontrib-bibtex==2.6.1 sphinxcontrib-katex==0.9.4 sphinxcontrib-youtube==1.2.0 sphinx-copybutton==0.5.2