[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: [Orekit Developers] issue 324 is probably a serious ugly bug

Hi Luc,

On Wed, 2017-05-24 at 17:55 +0200, MAISONOBE Luc wrote:
Hi all,

I am working on fixing remaining bugs on the issue tracker prior to  
9.0 release.
One of the issues is issue 324 <https://www.orekit.org/forge/issues/324> that
states that BoxAndSolarArraySpacecraft does not consider the lift force. This
was first registered as a feature request, but I considered it a bug and
changed its nature to bug some time ago.

Now it seems to me that the problem is *much* more serious than initially
envisioned. I think BoxAndSolarArraySpacecraft not only forgot about the
across-flux component, but it also computed the along-flux component wrong!
I would like to have other developers have a look at my reasoning. It is
really elementary physics, but I would like some confirmation, as the current
code does not behave like this.

Let's consider a particle arriving at an angle on a plate and bouncing:

          \  | plate
           \ |
    n <---  /|
           / |
          /  |
         i   |

The particle starts from the lower left, goes upward and to the right
along incident vector i, bounces, and changes its path upward and to the
left along reflected vector r. The plate unit normal is n, here oriented
to the left.

The effect of the drag on the plate is the opposite of the variation of
the linear momentum of the particle, and hence can be computed as
-mV (r - i). Considering an elastic collision, r - i is collinear to n.
If we consider the incident angle theta = pi - angle(i, n), we see
that r-i = 2cos(theta) n.

The force resulting from millions of particles hitting this plate will
be proportional to the cross section of the plate considering the flux
direction, which means this cross section will be A * cos(theta) is A
is the total surface of the plate.

The fact cos(theta) appears twice implies that the contribution of one
plate to the drag is proportional to the square of cos(theta), and is
always collinear to the plate normal.

The contribution of all the facets is only a summation of individual
contributions (we do not consider cast shadows in this class, this would
imply building a full 3D model).

The previous simple elastic collision model is what I have implemented
to solve issue 324.

During validation of this new implementation, I was puzzled when
comparing with the previous implementation. After some head-scratching,
I now think this is because our former implementation was plain wrong.

Basically, what the former implementation did was compute the cos(theta)
factor (using a simple dot product), and apply it to the incoming flux
in the flux direction. Even if we consider that the 2 in 2cos(theta)
can be removed from the model and included in the overall drag coefficient
(which seems correct to me, so C will remain close to 2 which is often
a default value), there is a square missing and the direction is wrong.

Unfortunately, the wrong direction does not show up much in general
computation because when bodies are convex and symmetrical, plates on the
left hand side and right hand side as seen from the incoming flux will
create forces cancelling out across-flux and adding together along flux.
The global resultant is therefore really along flux (there is no lift
for symmetrical bodies). So we end up with the correct direction, but is
is only by chance. However, the module of the contribution is wrong for
facets inclined. For head-on flux on a cubic body cos(theta)=1 and there
is no error. For flux along the diagonal of a cubic body, the three
visible facets all have cos(theta)=1/sqrt(3), and the current code
uses a factor 0.577 when 0.333 should be used.

What do you think?

In Section 8.2 in [1] Hughes examines the different interactions the
rarefied atmosphere can have with a spacecraft panel. He agrees with
you on the cosine squared term for specular reflections. He also states
(without providing a reference or much explanation) that in most
interactions the gas particle is absorbed by the surface and then
diffusely emitted latter at a much lower velocity. Ignoring the particles leaving the surface in diffuse reflection there is only one cosine term (due to the shrinking cross sectional area) and the resultant force is always along the relative velocity vector. So I think class is correctly implementing a multi-plate model assuming that the surface absorbs all of the incoming momentum from an air particle. 

In practice detailed drag modeling is hampered by the fact that we're really bad at predicting the atmospheric density. Perhaps a more detailed model would provide a coefficient for absorption vs. reflection incoming particle, but how much of a difference it makes is going to be situation dependent given all of the other uncertainties.

[1] Hughes, Peter C. Spacecraft Attitude Dynamics. Dover Publications.

best regards,

Attachment: smime.p7s
Description: S/MIME cryptographic signature