Method for computing self-contamination processes of a spacecraft

Abstract

A method for computing self-contamination processes of a spacecraft by means of a data processing device comprising the following steps: receiving a first set of input parameters comprising general definitions of the spacecraft, receiving a second set of input parameters comprising control parameters for the spacecraft orbital data, physics, numeric, and a predetermined accuracy requirement of the computation, computing a self-contamination process of the spacecraft based on the received first and second sets of input data by either evaluating the analytical solution of a basic equation of emission or numerically solving the basic equation of emission for calculating a deposit of molecules outgassed from surfaces of the spacecraft with a numerical solver with the data processing device, wherein the numerical solver applies an adaptive stepsize control based on the preset accuracy requirement of the computation, and outputting the calculated deposit.

Claims

1. A method for computing self-contamination processes of a spacecraft by means of a data processing device comprising the following steps: receiving a first set of input parameters comprising general definitions of the spacecraft, receiving a second set of input parameters comprising control parameters for the spacecraft orbital data, physics, numeric, and a predetermined accuracy requirement of the computation, computing a self-contamination process of the spacecraft based on the received first and second sets of input data by either evaluating the analytical solution of an equation of emission or numerically solving the equation of emission ${\stackrel{.}{W}}_{\mathrm{dep},i}=\left(\frac{F_{\mathrm{ji}}-{?}_{\mathrm{ij}}}{{?}_{r,j}}\right)W_{\mathrm{dep},j}+F_{\mathrm{ji}}{\stackrel{.}{W}}_{\mathrm{emitted},j}$ for calculating a deposit of molecules outgassed from surfaces of the spacecraft with a numerical solver with the data processing device, wherein the numerical solver applies an adaptive stepsize control based on the predetermined accuracy requirement of the computation, wherein F.sub.ji denotes a view factor from node j, first surface mesh element, to node i, second surface mesh element, ?.sub.ij denotes the Kronecker delta from node i to node j, ?.sub.rj denotes the temperature-dependent reemission time of node j, ?.sub.ej denotes the temperature-dependent emission time of node j, and wherein calculation of a deposit of molecules outgassed from surfaces of the spacecraft with a numerical solver comprises a stepwise analytical integration of the emission term W.sub.emittable of the equation of emission according to the following equation: $?W_{\mathrm{emittable},j}=W_{\mathrm{emittable},j}\left(t_n\right)\left(1-{\left(\frac{{?}_{e,j}\left(t_n\right)}{{?}_{e,j}\left(t_n\right)+{?}_{j}h}\right)}^{\frac{1}{{?}_j}}\right),$ with ?.sub.j being the secant slope defined by ${?}_j=\frac{{?}_{e,j}\left(t_{n+1}\right)-{?}_{e,j}\left(t_n\right)}{t_{n+1}-t_n},$ ?.sub.e designating time dependent material parameters, and h designating the time interval h=t.sub.n+1?t.sub.n, and outputting the calculated deposit; wherein the numerical solver applies a theta-method for calculating an analytical solution of the equation of emission for calculating a deposit of molecules outgassed from surfaces of the spacecraft and adapts the stepsize of the theta-method by performing the following steps: a) choosing an initial stepsize h.sub.lte and setting it as actual stepsize, b) calculating the local truncation error LTE of the theta-method for one step with the actual stepsize, c) calculating the stepsize for the next step depending on the calculated local truncation error LTE and setting the calculated stepsize for the next step as actual stepsize, d) repeating steps b) and c); and wherein calculating of the stepsize h.sub.tol for the next step depending on the calculated local truncation error is performed according to the following equation: $h_{\mathrm{tol}}=\min \left(h_{\mathrm{lte}}\sqrt{\frac{{\mathrm{TOL}}_m}{{\mathrm{LTE}}_m\left({?}_{\mathrm{prop}}\right)}}\right)$ with h.sub.lte designating the chosen initial stepsize, ?.sub.prop being either ?.sub.1 or ?.sub.2 depending on which one is selected as a propagator, LTE.sub.m(?.sub.prop) designating the local truncation error calculated for a pair of nodes W.sub.m, V.sub.m for the selected theta-method ?.sub.prop, wherein TOL.sub.m designates the predetermined accuracy requirements or tolerances that are applied at time t.sub.n+1 in the calculation of the next timestep and are selected from two predetermined values TOL.sub.m,abs and TOL.sub.rel according to the following rule: if U.sub.m(t.sub.n+1,?.sub.prop)>TOL.sub.m,abs
TOL.sub.m=TOL.sub.relu.sub.m(t.sub.n+1,?.sub.prop)
else
TOL.sub.m=TOL.sub.m,abs where U.sub.m(t.sub.n+1,?.sub.prop) is either W.sub.m(t.sub.n+1,?.sub.1) or V.sub.m(t.sub.n+1,?.sub.2) depending on which one is selected as propagator, and wherein $w\left(t_{n+1},{?}_1\right)=y\left(t_n\right)+h\left[{?}_{1}f\left(t_{n+1},y\left(t_{n+1}\right)\right)+\left(1-{?}_1\right)f\left(t_n,y\left(t_n\right)\right)\right];$ $v\left(t_{n+1},{?}_2\right)=y\left(t_n\right)+h\left[{?}_{2}f\left(t_{n+1},y\left(t_{n+1}\right)\right)+\left(1-{?}_2\right)f\left(t_n,y\left(t_n\right)\right)\right];$ $\mathrm{and}$ ${\mathrm{LTE}}_m\left({?}_{\mathrm{prop}}\right)=\left(\frac{1}{2}-{?}_{\mathrm{prop}}\right)\frac{w_m\left(t_n,{?}_1\right)-v_m\left(t_n,{?}_2\right)}{{?}_1-{?}_2},$ wherein y(t.sub.n) denotes the exact solution of the LTE at a time t.sub.n, wherein y(t.sub.n+1) denotes the exact solution of the LTE at a next time t.sub.n+1 and wherein
f(t.sub.n+1,y(t.sub.n+1))={dot over (y)}(t.sub.n+1).

2. The method of claim 1, wherein the calculating of the local truncation error LTE of the theta-method comprises calculating a first local truncation error LTE.sub.explicit using the explicit Euler method, calculating a second local truncation error LTE.sub.implicit using the implicit Euler method, and selecting as the local truncation error LTE the one of the calculated first and second local truncation error LTE.sub.explicit and LTE.sub.implicit, with the smaller magnitude.

3. The method of claim 1, wherein the calculating of the local truncation error LTE of the theta-method is performed according to the following equation: ${\mathrm{LTE}}_m\left({?}_{\mathrm{prop}}\right)=\left(\frac{1}{2}-{?}_{\mathrm{prop}}\right)\frac{w_m\left(t_{n+1},{?}_1\right)-v_m\left(t_{n+1},{?}_2\right)}{{?}_1-{?}_2}$ with ?.sub.1, ?.sub.2 designating different theta-methods, ?.sub.prop being either ?.sub.1 or ?.sub.2 depending on which one is selected as a propagator, W.sub.m, V.sub.m designating the deposit on a pair of nodes or mesh elements of the surface of the spacecraft, which is represented as a collection of surfaces, which are subdivided into mesh elements with outgassing properties.

4. The method of claim 1, wherein ?.sub.1, ?.sub.2 are set to the explicit and implicit Euler method with ?.sub.1=0, ?.sub.2=1 the implicit Euler method ?.sub.2 is selected as propagator ?.sub.prop.

5. The method of claim 2, wherein the theta-method comprises a higher order linear multistep method, which is used to calculate the results and/or propagate the solution.

6. The method of claim 1, wherein the calculating of the stepsize for the next step depending on the calculated local truncation error LTE comprises the following steps: if the local truncation error LTE is smaller than or equal to the predetermined accuracy requirement TOL performing the following steps: accepting the step, calculating an increased stepsize for the next step depending on the calculated local truncation error LTE and setting the increased stepsize as actual stepsize, and if the local truncation error LTE is larger than the predetermined accuracy requirement TOL performing the following steps: discarding the step, calculating a decreased stepsize for this step depending on the calculated local truncation error LTE and setting the decreased stepsize as actual stepsize.

7. The method of claim 1, wherein the equation used for the stepwise analytical integration of the emission term W.sub.emittable of the equation of emission is expanded for small values of ?.sub.j to $?W_{\mathrm{emittable},j}?W_{\mathrm{emittable},j}\left(t_n\right)\left(1-\exp \left(-\frac{?t}{{?}_{e,j}\left(t_n\right)}\right)\left(1+\frac{{\left(?t\right)}^2{?}_j}{2{?}_{e,j}^2\left(t_n\right)}\right)\right).$

Description

BRIEF DESCRIPTION OF DRAWINGS

(1) The various embodiments will hereinafter be described in conjunction with the following drawing figures, wherein like numerals denote like elements, and:

(2) FIG. 1 shows schematically the process of self-contamination by means of two surfaces i and j;

(3) FIG. 2 shows a graph with a timestep dt for demonstrating the principle the method according to the embodiment;

(4) FIG. 3 shows a table of test results compared to analytical results calculated for two nodes after a first timestep with a computer program implementing an embodiment of a method for computing self-contamination processes of a spacecraft according to an embodiment;

(5) FIGS. 4-7 show diagrams of the computation results of the method used to compute self-contamination for the full BepiColombo geometry with the complete trajectory; and

(6) FIG. 8 shows a diagram of the development of the relative error as compared to the analytical solution given a relative accuracy of 0.1% of an embodiment for a geometry whose nodes are at constant but different temperature.

(7) FIG. 9 is a flow chart that illustrates a method for computing self-contamination processes of a spacecraft.

DETAILED DESCRIPTION

(8) 0. Basics of the Computation of Self-Contamination Processes of a Spacecraft

(9) For computing self-contamination processes of a spacecraft, which particularly comprises different materials, the spacecraft is represented as a collection of surfaces, which are sub-divided into nodes (mesh elements) that have varying outgassing properties and emit and exchange contaminant mass. The contaminant molecules emitted by the surfaces and the bulk of the spacecraft move towards all spacecraft surfaces.

(10) For computing a self-contamination process of a spacecraft with a data processing device such as a computer, which executes a computer program implementing basic steps of the embodiment, the data processing device receives two sets of input parameters. A first set of input parameters contains general definitions of the spacecraft such as parameters of the geometry of the spacecraft, properties of surfaces of the spacecraft, particularly outgassing/emission and reemission properties of the surfaces or coatings of the surfaces, and temperatures. A second set of input parameters contains control parameters for the spacecraft orbital data, physics, numeric, and predetermined accuracy requirements. The spacecraft orbital data typically comprise orbit definition and interpolation data and kinematics data. The control parameters for physics comprise data for atmosphere definition and transport definition. The control parameters for numeric may comprise parameters for raytracing and view factors, parameters regarding surfaces sensitive to contamination and to outgassing, and/or volume mesh data.

(11) The data processing device processes the received input data for computing a self-contamination process of the spacecraft. The computation comprises calculating an analytical solution of a basic equation for calculating a deposit of molecules outgassed from surfaces of the spacecraft. The computation may be divided in two basic steps:

(12) computation of view factors, i.e. the proportion of the total mass emitted by a node i and reaching node j, and

(13) solution of the integral equation of the deposit.

(14) The view factors F.sub.ij between two surfaces mesh elements (nodes) i and j are calculated in a first step. The view factors can be defined at macroscopic level as the mass or the molecules emitted by node or surface mesh element i, which reach(es) node or surface mesh element j, divided by the entire mass or all molecules emitted by node or surface mesh element i. At a microscopic level, view factors F.sub.ij can be defined as the probability for a molecule emitted by node i reaching node j.

(15) As shown in FIG. 1, the transfer of mass between two nodes i and j comprises direct flux, flux scattered by other contaminant molecules (self-scattering), flux scattered by ambient gas flux (ambient gas scattering), and flux reemitted from spacecraft surfaces.

(16) The view factors F.sub.ij can be determined by a Monte-Carlo ray-tracing as briefly explained in the following: a node i randomly emits a large number of rays, which represent the outgassed molecules. The emission of rays can be performed as follows: emission points are sampled according to an iso-distribution law; emission directions are sampled according to Lambert's law. The mean free path ?.sub.a is based on Robertson formula where the cross-section is computed by the Hard Sphere model. For each ray, a collision distance d.sub.c is sampled at random according to the collision probability law P(x)=e.sup.?x/?.sup.a, where x is the distance and ?.sub.a is the mean free path. By comparing the distance, the molecule is either scattered or impinges the node j. If d.sub.c<d.sub.j, a collision is accepted, otherwise it is rejected. If a collision is rejected, the molecule reaches node j and contributes to F.sub.ij. If a collision is accepted, the molecule is reemitted from the collision point with a random direction and reaches the node k and contributes to F.sub.ik.

(17) After calculating the view factors F.sub.ij, the mass deposit can be calculated as briefly outlined in the following: the nodes i, j, k, . . . form a network where each pair of elements i, j or i, k or j, k is related by their view factors F.sub.ij, F.sub.ik or F.sub.jk. Each element has emission, deposition and reemission of deposited contaminants that depend on the temperature. Thus, this problem becomes an inhomogeneous system of differential equations of first order with time dependent coefficients, which can be solved by the inventive solution as explained in the following in detail.

(18) The computation steps according to the inventive solution and as described in the following can be implemented in a computer program such as an additional module or plug-in for the above mentioned COMOVA-software and can significantly reduce the computation time required for computing self-contamination processes of spacecraft consisting of different materials with different and time-dependent outgassing and desorbing behavior.

(19) Calculation of Deposits

(20) 0. The Basic Equation for Deposit

(21) The equation (as defined in the OUTGASSING Software User's Manual, internal document identification MOS.UM.93694.ASTR, on page 136, equation (72),

(22) $W_{d,i}^{}=W_{d,i}^n{?}_i+\frac{S_i^n}{A_{i,j}}\left({?}_i-1\right)+\underset{j}{.Math.}\frac{F_{j,i}{W}_{e,j}^n}{1+A_{i,j}{?}_{e,j}}\left({?}_i-\exp \left(-\frac{?t}{{?}_{e,j}}\right)\right);$
simplified form of the equation can be found in Van Eesbeek M., Zwaal A., Outgassing and Contamination Model Based on Residence Time, Proceedings of the Third European Symposium on Spacecraft Materials in Space Environment, Noordwijk, October 1985, (ESA SP-232)) for the deposit on a node i that is to be solved is as follows:

(23) ${\stackrel{}{W}}_{\mathrm{dep},i}=\left(\frac{F_{\mathrm{ji}}-{?}_{\mathrm{ij}}}{{?}_{r,j}}\right)W_{\mathrm{dep},j}+F_{\mathrm{ji}}{\stackrel{}{W}}_{\mathrm{emitted},j}$${\stackrel{}{W}}_{\mathrm{emitted},j}=\frac{1}{{?}_{e,j}}{W}_{\mathrm{emittable},j}=-{\stackrel{}{W}}_{\mathrm{emittable},j}$

(24) wherein F.sub.ji denotes a view factor from node j (first surface mesh element) to node i (second surface mesh element), ?.sub.r,j denotes the temperature-dependent reemission time of node j, ?.sub.e,j denotes the temperature-dependent emission time of node j. In the following, the derivation of a numerical solver for SYSTEMA/OUTGASSING that applies adaptive stepsize control based on accuracy requirements to minimize calculation time according to the embodiment is described. A general formulation of the numerical scheme is presented using the theta (-method. It allows solving the deposit equation for time dependent material parameters. The explicit and implicit first order Euler schemes were identified as the most promising combination. The algorithms according to an embodiment were implemented in MATLAB and C++ and applied to the various models.

(25) 2 Calculation of the Local Truncation Error (LTE) of the Theta (?)-Method

(26) The principle of the theta-method is depicted by FIG. 2:

(27) The timestep dt is split into two parts. During the first part dt(1??) the solution is propagated along the slope at time t. At time t+dt(1??) the solution is propagated along the slope at time t+dt.

(28) The explicit as well as the implicit Euler-scheme can be considered as special cases of the ?-method.

(29) Thus, the LTE for this method can be derived first and the LTEs for the two Euler methods can be deduced from it.

(30) Given the Initial Value Problem
{dot over (y)}(t)=f(t, y(t)), t ?[0, T.sub.final], y(0)=y.sub.0

(31) and assuming that at time t.sub.n=nh the solution y(t.sub.n) is exact, the local truncation error LTE is defined by
y(t.sub.n+1)=y(t.sub.n)+h[?f(t.sub.n+1, y(t.sub.n+1))+(1??)f(t.sub.n, y(t.sub.n))]+LTE (5)

(32) Observing that f(t.sub.n+1, y(t.sub.n+1))={dot over (y)}(t.sub.n+1) and expanding forward

(33) $\begin{array}{cc}y\left(t_{n+1}\right)=y\left(t_n\right)=h\stackrel{}{y}\left(t_n\right)+\frac{1}{2}{h}^2\stackrel{}{y}\left(t_n\right)+\frac{1}{6}{h}^3\stackrel{}{y}\left(t_n\right)+O\left(h^4\right)& \left(6\right)\end{array}$

(34) and using

(35) $\begin{array}{cc}\frac{d}{dt}y\left(t_{n+1}\right)=\stackrel{}{y}\left(t_n\right)+h\stackrel{}{y}\left(t_n\right)+\frac{1}{2}{h}^2\stackrel{}{y}\left(t_n\right)+O\left(h^4\right)& \left(7\right)\end{array}$

(36) Solving eq. (5) for LTE and inserting eqs. (6) and (7) yields

(37) $\begin{array}{cc}\mathrm{LTR}=\stackrel{}{y}\left(t_n\right)h^2\left(\frac{1}{2}-?\right)+\stackrel{}{y}\left(t_n\right)h^3\left(\frac{1}{6}-\frac{?}{2}\right)+O\left(h^4\right)& \left(8\right)\end{array}$

(38) In case

(39) $\frac{1}{2}=?$
the trapezoidal rule is get, which is consistent of order 2 with LTE?h.sup.3, in all other cases the ?-method is of order 1.

(40) Now, the LTEs for the Euler methods can be calculated:

(41) Explicit Euler method (?=0):

(42) $\begin{array}{cc}\mathrm{LTE}=\frac{1}{2}\stackrel{}{y}\left(t_n\right)h^2+O\left(h^3\right)& \left(9\right)\end{array}$

(43) Implicit Euler method (?=1):

(44) $\begin{array}{cc}\mathrm{LTE}=-\frac{1}{2}\stackrel{}{y}\left(t_n\right)h^2+O\left(h^3\right)& \left(10\right)\end{array}$

(45) 0. Calculation of the stepsize using the local truncation error. The LTE can be written in the form of a remainder term in the Taylor series. Two series with two different values

(46) $\frac{1}{2}?{?}_1,\frac{1}{2}?{?}_2$
can be defined:

(47) Assuming that y(t.sub.n) was the exact solution of eq. (5), the exact solution at time t.sub.n+1 can be expanded as a Taylor series with remainder:

(48) 0$\begin{array}{cc}y\left(t_{n+1},{?}_1\right)=y\left(t_n\right)+h\left[{?}_{1}f\left(t_{n+1},y\left(t_{n+1}\right)\right)+\left(1-{?}_1\right)f\left(t_n,y\left(t_n\right)\right)\right]+\stackrel{}{y}\left({?}_1\right)h^2\left(\frac{1}{2}-{?}_1\right)\mathrm{and}y\left(t_{n+1},{?}_2\right)=y\left(t_n\right)+h\left[{?}_{2}f\left(t_{n+1},y\left(t_{n+1}\right)\right)+\left(1-{?}_2\right)f\left(t_n,y\left(t_n\right)\right)\right]+\stackrel{}{y}\left({?}_2\right)h^2\left(\frac{1}{2}-{?}_2\right)\mathrm{with}{?}_i?\left[t_n,t_{n+1}\right].& \left(11\right)\end{array}$

(49) As the remainder term is unknown we can only approximate the solution at time t.sub.n+1:
w(t.sub.n+1, ?.sub.1)=y(t.sub.n)+h[?.sub.1f(t.sub.n+1, y(t.sub.n+1))+(1??.sub.1)f(t.sub.n, y(t.sub.n))](11)
v(t.sub.n+1, ?.sub.2)=y(t.sub.n)+h[?.sub.2f(t.sub.n+1, y(t.sub.n+1))+(1??.sub.2)f(t.sub.n, y(t.sub.n))](12)

(50) Subtracting eqs. (11) from (12):

(51) $\begin{array}{cc}w\left(t_{n+1},{?}_1\right)-y\left(t_{n+1},{?}_1\right)=-\stackrel{}{y}\left(?\right)h^2\left(\frac{1}{2}-{?}_1\right)v\left(t_{n+1},{?}_2\right)-y\left(t_{n+1},{?}_2\right)=-\stackrel{}{y}\left(?\right)h^2\left(\frac{1}{2}-{?}_2\right)& \left(13\right)\end{array}$

(52) Here it was used that for sufficiently smooth functions ?.sub.1??.sub.2??.

(53) As it has been assumed that y(t.sub.n+1, ?.sub.i) was the exact solution at time t.sub.n+1, eqs. (13) can be subtracted from each other to calculate an estimate for the second derivative:

(54) $\begin{array}{cc}\stackrel{}{y}\left(?\right)h^2?\frac{w\left(t_{n+1},{?}_1\right)-v\left(t_{n+1},{?}_2\right)}{{?}_1-{?}_2}& \left(14\right)\end{array}$

(55) Comparing eqs. (8) with eqs. (14) it can be seen (only the smallest power of h is relevant):

(56) $\begin{array}{cc}\mathrm{LTE}\left({?}_i\right)=\stackrel{}{y}\left(?\right)h^2\left(\frac{1}{2}-{?}_i\right)=\left(\frac{1}{2}-{?}_i\right)\frac{w\left(t_{n+1},{?}_1\right)-v\left(t_{n+1},{?}_2\right)}{{?}_1-{?}_2}& \left(15\right)\end{array}$

(57) The truncation error of any of the two series expansions at a certain step can thus be expressed by the difference of the two calculated function values at that step.

(58) Using the explicit and the implicit Euler method (?.sub.1=0, ?.sub.2=1) yields truncation errors of

(59) $\begin{array}{cc}\begin{array}{c}{\mathrm{LTE}}_{\mathrm{explicit}}=\left(\frac{1}{2}-{?}_1\right)\frac{w\left(t_{n+1},{?}_1\right)-v\left(t_{n+1},{?}_2\right)}{{?}_1-{?}_2}\\ =-\frac{1}{2}(w\left(t_{n+1},{?}_1\right)-v\left(t_{n+1},{?}_2\right)\end{array}\mathrm{and}\begin{array}{c}{\mathrm{LTE}}_{\mathrm{implicit}}=\left(\frac{1}{2}-{?}_2\right)\frac{w\left(t_{n+1},{?}_1\right)-v\left(t_{n+1},{?}_2\right)}{{?}_1-{?}_2}\\ =\frac{1}{2}(w\left(t_{n+1},{?}_1\right)-v\left(t_{n+1},{?}_2\right)\end{array}& \left(16\right)\end{array}$

(60) for each of the expansions.

(61) In the above example both LTE's have the same magnitude. Usually one would chose the LTE which has the smallest magnitude to determine the series expansion along which the solution is propagated. In the sequel this expansion shall be called propagator and the other one estimator. In addition to the LTE the stability has to be taken into account. In the above example, although the magnitudes of both LTEs are equal, one would choose the implicit expansion to propagate the solution because the explicit Euler scheme becomes instable if the timesteps become too large.

(62) Provided that

(63) $\frac{1}{2}?{?}_i$
the LTE can made arbitrarily small by reducing the stepsize h. For practical purposes it is sufficient that the LTE does not exceed a given tolerance TOL. The local truncation error TOL that is compatible with the tolerance requirement can be realized by a certain timestep h.sub.tol so that

(64) $\begin{array}{cc}\mathrm{TOL}=\stackrel{}{y}\left({?}_{\mathrm{tol}}\right)h_{\mathrm{tol}}^2\left(\frac{1}{2}-{?}_i\right)\mathrm{while}\mathrm{LTE}\left({?}_i\right)=\stackrel{}{y}\left({?}_{\mathrm{lte}}\right)h_{\mathrm{lte}}^2\left(\frac{1}{2}-{?}_i\right)& \left(17\right)\end{array}$

(65) Now requesting LTE(?.sub.i )=TOL and assuming ?(?.sub.tol)??(?.sub.lte) leads to

(66) $\begin{array}{cc}{h}_{\mathrm{tol}}=h_{\mathrm{lte}}\sqrt{\frac{\mathrm{TOL}}{\mathrm{LTE}\left({?}_i\right)}}& \left(18\right)\end{array}$

(67) with the LTE given by eq. (16).

(68) This equation defines the strategy for choosing the stepsize:

(69) Chose a stepsize h.sub.lte and calculate the pertaining LTE

(70) If LTE<=TOL accept the step, calculate the next stepsize by eq. (18)

(71) If LTE>TOL repeat the step, calculate the stepsize for the retry by eq. (18)

(72) This strategy may not always be successful. If the stiffness increases as compared to the last step the LTE may exceed the tolerance and the calculation will be repeated with a smaller timestep. If the LTE is small enough, the step is accepted and the next step is calculated using the timestep of the last successful step. If not, the timestep is decreased once more.

(73) If on the other hand the LTE is smaller than the tolerance, the next timestep is increased according to eq. (18).

(74) The method outlined is capable of detecting stiffness changes and can adapt itself. It is important to note, that the calculation of the next timestep depends on the LTE of the current timestep and thus is a mere estimate.

(75) 0. Stepsize Strategy and to Tolerance Criteria for a System of Differential Equations

(76) If a system of equations has to be solved, a LTE can be calculated for each nodepair , according to eq. (16).

(77) In order to get the process started a stepsize h.sub.lte is chosen that gives rise to a LTE according to eq. (16).

(78) Then at time t.sub.n a vector of LTEs is calculated from all propagator-estimator pairs according to

(79) $\begin{array}{cc}{\mathrm{LTE}}_m\left({?}_{\mathrm{prop}}\right)=\left(\frac{1}{2}-{?}_{\mathrm{prop}}\right)\frac{w_m\left(t_n,{?}_1\right)-v_m\left(t_n,{?}_2\right)}{{?}_1-{?}_2}& \left(19\right)\end{array}$

(80) where ?.sub.prop is either ?.sub.1 or ?.sub.2 depending on which one was chosen as the propagator.

(81) The next stepsize is then taken to be the minimum of all that were calculated:

(82) $\begin{array}{cc}{h}_{\mathrm{tol}}=\min \left(h_{\mathrm{lte}}\sqrt{\frac{\mathrm{TOL}}{{\mathrm{LTE}}_m\left({?}_{\mathrm{prop}}\right)}}\right)& \left(20\right)\end{array}$

(83) In the above formula TOL does not necessarily have to be a scalar, it can be specified for each node so that

(84) 0$\begin{array}{cc}{h}_{\mathrm{tol}}=\min \left(h_{\mathrm{lte}}\sqrt{\frac{{\mathrm{TOL}}_m}{{\mathrm{LTE}}_m\left({?}_{\mathrm{prop}}\right)}}\right)& \left(21\right)\end{array}$

(85) The tolerances TOL.sub.m that are applied at time t.sub.n in the calculation of the next timestep are calculated from two user supplied values TOL.sub.m,abs and TOL.sub.rel according to the following rule:

(86) If u.sub.m(t.sub.n, ?.sub.prop)>TOL.sub.m,abs
TOL.sub.m=TOL.sub.relu.sub.m(t.sub.n, ?.sub.prop) (22)

(87) else TOL.sub.m=TOL.sub.m,abs

(88) where u.sub.m(t.sub.n, ?.sub.prop) is either w.sub.m(t.sub.n, ?.sub.1) or v.sub.m(t.sub.n, ?.sub.2) depending on which one was selected as propagator. Eq. (22) implies that the tolerance requirements relax when the solution is close to zero.

(89) 0. Discretizing the Basic Equation for Deposit

(90) i. Direct Integration of the Emission Term

(91) In the case of N outgassing nodes the equation to be solved has the form:

(92) $\begin{array}{cc}{\stackrel{}{W}}_{\mathrm{dep},i}=\left(\frac{F_{\mathrm{ji}}-{?}_{\mathrm{ij}}}{{?}_{r,j}}\right)W_{\mathrm{dep},j}+F_{\mathrm{ji}}{\stackrel{}{W}}_{\mathrm{emitted},j}\mathrm{with}& \left(23\right)\\ {\stackrel{}{W}}_{\mathrm{emitted},j}=\frac{1}{{?}_{e,j}}{W}_{\mathrm{emittable},j}=-{\stackrel{}{W}}_{\mathrm{emittable},j}& \left(24\right)\end{array}$

(93) Eq. (23) can be written more conveniently using matrix notation:
.sub.dep,i=G.sub.ilR.sub.lmW.sub.dep,m+F.sub.jiE.sub.jlW.sub.emittable,l (25)

(94) Where G.sub.il=F.sub.li??.sub.il is the geometriematrix and R.sub.lm, E.sub.jl are diagonalmatrices of the form

(95) $R_{\mathrm{lm}}=\left(\begin{array}{cccc}{?}_{r,1}^{-1}& 0& & 0\\ 0& \mathrm{.Math.}& & \\ & & \mathrm{.Math.}& 0\\ 0& & 0& {?}_{r,N}^{-1}\end{array}\right)$$\mathrm{and}$$E_{\mathrm{jl}}=\left(\begin{array}{cccc}{?}_{e,1}^{-1}& 0& & 0\\ 0& \mathrm{.Math.}& & \\ & & \mathrm{.Math.}& 0\\ 0& & 0& {?}_{e,N}^{-1}\end{array}\right)$

(96) Matrices R.sub.lm and E.sub.jl are time dependent because they contain the time dependent material parameters ?.sub.r, ?.sub.e. Time dependency enters through the temperature of the node which depends on the trajectory and can be highly nonlinear. In such a case an analytical expression for W.sub.emittable cannot be found and the equations for W.sub.emittable and W.sub.dep have to be solved simultaneously:

(97) Introducing a New Variable

(98) W=(W.sub.dep,W.sub.emittable).sup.T eqs. (23) and (24) can now be written as a linear, homogeneous system of differential equations:

(99) $\begin{array}{cc}\stackrel{}{W}=\left(\begin{array}{cc}\mathrm{GR}& \mathrm{FE}\\ 0& -E\end{array}\right)\left(\begin{array}{c}{W}_{\mathrm{dep}}\\ W_{\mathrm{emittable}}\end{array}\right)=\mathrm{MW}& \left(26\right)\end{array}$

(100) Applying the theta-method to this equation results in a recursion formula for W(t.sub.n+1) at time t.sub.n+1=(n+1)h
W.sub.l(t.sub.n+1, ?.sub.i)=W.sub.l(t.sub.n)+h[?.sub.iM.sub.lm(t.sub.n+1)W.sub.m(t.sub.n+1)+(1??.sub.i)M.sub.lm(t.sub.n)W.sub.m(t.sub.n)](27)

(101) The relative and absolute tolerances needed to calculate the tolerance requirement of every node according to eq. (22), however, would now be applied to all nodes. This is a drawback because even in the case of constant temperature equation (24) describing the emission is always only calculated with the precision goal required for the deposit while it could have been calculated to arbitrary precision using the analytical solution. Thus a lot of precision is lost while at the same time the computation time is increased disproportionately because the matrix size has quadrupled.

(102) i. Stepwise Analytical Integration of the Emission Term

(103) 1. General Solution

(104) If, on the other hand, a solution for W.sub.emittable were available, the computation time could be decreased significantly because the theta-method usually requires the solution of an equation system due to the implicit character of the method. Halving the number of equations can thus reduce the computation time by up to a factor of 8.

(105) It shall now be derived an analytical solution for W.sub.emittable based on the linear approximation of ?.sub.e in the time interval h=t.sub.n+1?t.sub.n where h are the timesteps that are calculated according to eq. (21).

(106) Eq. (23) can be written in the form
{dot over (W)}.sub.dep,i=A.sub.ij(t)W.sub.dep,j?F.sub.ji{dot over (W)}.sub.emittable,j (28)

(107) The equation can be integrated in the interval t ? [t.sub.N, t.sub.N+1]:

(108) $\begin{array}{cc}{?}_{t_n}^{t_{n+1}}d{W}_{\mathrm{dep},i}={?}_{t_n}^{t_{n+1}}{A}_{\mathrm{ij}}\left(t\right)W_{\mathrm{dep},j}dt-F_{\mathrm{ij}}{?}_{t_n}^{t_{n+1}}{\stackrel{}{W}}_{\mathrm{emittable},j}dt& \left(29\right)\end{array}$

(109) Applying the theta-method to the first term on the right side:

(110) $\begin{array}{cc}{W}_{\mathrm{dep},i}\left(t_{n+1}\right)-W_{\mathrm{dep},i}\left(t_n\right)=h\left[?A_{\mathrm{ij}}\left(t_{n+1}\right)W_{\mathrm{dep},j}\left(t_{n+1}\right)+\left(1-?\right)A_{\mathrm{ij}}\left(t_n\right)W_{\mathrm{dep},j}\left(t_n\right)\right]-F_{ji}{?}_{t_n}^{t_{n+1}}d{W}_{\mathrm{emittable},j}& \left(30\right)\\ \mathrm{with}{F}_{\mathrm{ji}}{?}_{t_n}^{t_{n+1}}d{W}_{\mathrm{emittable},j}=F_{ji}\left(W_{\mathrm{emittable},j}\left(t_{n+1}\right)-W_{\mathrm{emittable},j}\left(t_n\right)\right)=F_{\mathrm{ji}}{\mathrm{?W}}_{\mathrm{emittable},j}& \left(31\right)\end{array}$

(111) Here ?W.sub.emittable,j is the quantity of material emitted from the bulk during the time interval t.sub.n+1?t.sub.n which is then distributed in the structure through the visibility matrix.

(112) In order to calculate ?W.sub.emittable,j the linear transformation is introduced
t=t?t.sub.n, t ? [0,t.sub.n+1?t.sub.n]

(113) so that the emission equation takes the form

(114) $\begin{array}{cc}{\stackrel{.}{W}}_{\mathrm{emittable},j}\left(t^{}\right)=-\frac{W_{\mathrm{emittable},j}\left(t^{}\right)}{{?}_{e,j}\left(t^{}=0\right)+{?}_j{t}^{}}\mathrm{with}{W}_{\mathrm{emittable},j}\left(t^{}=0\right)=W_{\mathrm{emittable},j}\left(t_n\right)& \left(32\right)\end{array}$

(115) where ?.sub.j is the secant slope defined by

(116) ${?}_j=\frac{{?}_{e,j}\left(t_{n+1}\right)-{?}_{e,j}\left(t_n\right)}{t_{n+1}-t_n}$
and is constant in the interval. (33)

(117) Equation (32) can be integrated in closed form to yield:

(118) $\begin{array}{cc}{W}_{\mathrm{emittable},j}\left(t^{}\right)=W_{\mathrm{emittable},j}\left(t_n\right){\left(\frac{{?}_{e,j}\left(t_n\right)}{{?}_{e,j}\left(t_n\right)+{?}_j{t}^{}}\right)}^{\frac{1}{{?}_j}}\mathrm{As}\begin{array}{c}{W}_{\mathrm{emittable},j}\left(t=t_{n+1}\right)=W_{\mathrm{emittable},j}\left(t^{}=h\right)\\ =W_{\mathrm{emittable},j}\left(t_n\right){\left(\frac{{?}_{e,j}\left(t_n\right)}{{?}_{e,j}\left(t_n\right)+{?}_{j}h}\right)}^{\frac{1}{{?}_j}}\end{array}& \left(34\right)\end{array}$

(119) it is now possible to calculate the desired quantity ?W.sub.emittable,j:

(120) $\begin{array}{cc}?W_{\mathrm{emittable},j}=W_{\mathrm{emittable},j}\left(t_n\right)\left(1-{\left(\frac{{?}_{e,j}\left(t_n\right)}{{?}_{e,j}\left(t_n\right)+{?}_{j}h}\right)}^{\frac{1}{{?}_j}}\right)& \left(35\right)\end{array}$

(121) Starting at t=0 with the known quantity W.sub.emittable,j(t=0) the amount of emitted material during each timestep h can be recursively calculated and inserted into eq. (30) to calculate W.sub.dep,i(t.sub.n+1).

(122) 1. Approximation for Small Slopes

(123) In the form given by eq. (34) vanishing slopes lead to a singularity so that the case of constant temperature cannot be calculated. This can be overcome by expanding eq. (34) for small values of ?.sub.j:

(124) 0$\begin{array}{cc}?W_{\mathrm{emittable},j}?W_{\mathrm{emittable},j}\left(t_n\right)\left(1-\exp \left(-\frac{?t}{{?}_{e,j}\left(t_n\right)}\right)\left(1+\frac{{\left(?t\right)}^2{?}_j}{2{?}_{e,j}^2\left(t_n\right)}\right)\right)& \left(36\right)\end{array}$

(125) For vanishing slope the known analytical solution is thus recovered. In the numerical implementation a threshold value for the slope can be defined below which eq. (36) is used instead of eq. (34).

(126) 0. Selecting ?.sub.1 and ?.sub.2

(127) Numerical experiments showed that setting ?.sub.1=0, ?.sub.2=1 (explicit and implicit Euler method) with the implicit method as the propagator was the best tradeoff between speed and accuracy: Choosing ?.sub.1>0, ?.sub.2=1 yielded a slightly better accuracy at the expense of a doubled computation time because both equations are now implicit. Setting ?.sub.1?0, ?.sub.2<1 on the other hand produced negative deposits when the values approached zero because the propagator now is only A-stable instead of L-stable.

(128) The timesteps that the algorithm computes do not meet the stability criteria for an explicit method. In order to prevent the estimator from diverging, its new value was computed using the propagator value of the last timestep. The new values for propagator and estimator were computed according to the expression

(129) $\begin{array}{cc}{W}_{\mathrm{dep},i}\left(t_{n+1}\right)-W_{\mathrm{dep},i}\left(t_n\right)=h?A_{\mathrm{ij}}\left(t_{n+1}\right)W_{\mathrm{dep},j}\left(t_{n+1}\right)-F_{ji}{?}_{t_n}^{t_{n+1}}d{W}_{\mathrm{emittable},j}& \left(37\right)\end{array}$

EXAMPLES

(130) 2.1 Constant Temperature

(131) In a first test results calculated with the numerical solver were compared to those calculated with the analytical solver (Analytical solver for SYSTEMA/OUTGASSING ACE82.TN.ASTR.00004). The test case was an electronic box at a uniform constant temperature of 377? K.

(132) The table in FIG. 3 outlines the results for nodes 1 and 2 after the first timestep, in the middle and at the end of the iterations:

(133) Although in the numerical solution a relative error of only 1e-2 was required, up to 5 digits were computed correctly.

(134) The CPU times needed were 246 sec (analytical) vs. 7 sec (numerical).

(135) 2.2 Variable Temperature

(136) As an example for variable temperature the full BepiColombo geometry with the complete trajectory was chosen. The model was comprised of 9909 nodes. The aim was the calculation of deposit on the hottest node of the structure.

(137) The parameters were:

(138) Relative error: 1e-2

(139) Absolute error: 2.56e-9 kg/m.sup.2 (corresponds to 1e-11 kg on node 6132)

(140) Max. stepsize: 250 h The results are depicted in FIGS. 4-7

(141) The deposits faithfully follow the form of the temperature curve becoming lower when temperatures rise and vice versa (FIGS. 4 and 5). The timesteps become smaller when the temperatures rise and are limited to 250 h (FIG. 6). The relative error was computed by dividing the LTE by the deposit. The requested accuracy is maintained in a stable way except for the times when the timesteps are limited to 250 h: Here the algorithm wanted to use a larger timestep and answers to the limitation with a better accuracy (FIG. 7).

(142) The results were reproduced with a C++-implementation using the standard LAPACK library with the results being identical to at least 12 digits.

(143) Calculation Times:

(144) C++ (1CPU+standard LAPACK): 536055 sec

(145) MATLAB (4CPUs+multithreading+INTEL-MKL lib): 47290 sec

(146) 0. Improving Accuracy

(147) It is in principle not possible to calculate the absolute error of a numerical method because an analytical solution is not available, else one would have used it in the first place. For the problem at hand an analytical solution can be formulated (Analytical solver for SYSTEMA/OUTGASSING ACE82.TN.ASTR.00004) albeit it is valid only for constant temperatures and needs numerical linear algebra to calculate the diagonalization of the problemmatrix. The deviation of the numerically derived solution from the true solution can therefore be assessed and different numerical methods can be compared. This approach was in fact used to determine the choice of ?.sub.1=0, ?.sub.2=1 as the final setting in the implementation of the method. It had been proofed that through the use of the implicit Euler method as the propagator of the solution a stable behaviour could be achieved. Nevertheless the order of convergency is only one. To improve the order of convergency while at the same time retaining A-stability one's choice is limited: The only linear multistep method is the trapezium method which is of order two and exhibits A-stability. It was therefore decided to implement the trapezium method in addition to the two theta-methods and use its results to propagate the solution while the theta-methods were used to calculate the timesteps. Integration took place according to the following formula:

(148) $\begin{array}{cc}\left({?}_{\mathrm{ij}}-\frac{h}{2}{A}_{\mathrm{ij}}\left(t_{n+1}\right)\right)W_{\mathrm{dep},j}\left(t_{n+1}\right)=\left({?}_{\mathrm{ij}}+\frac{h}{2}{A}_{\mathrm{ij}}\left(t_n\right)\right)W_{\mathrm{dep},j}\left(t_n\right)-F_{\mathrm{ji}}{?}_{t_n}^{t_{n+1}}d{W}_{\mathrm{emittable},j}& \left(38\right)\end{array}$

(149) The integral on the right hand side was also calculated using the trapezium rule. To speed up the calculation the equations system was solved using the PARDISO solver which is a part of the INTEL MKL library. This solver is optimized for sparse matrices and speeds up the calculation by a factor of up to 20.

(150) 0. Results

(151) The geometry of an electronic box was used as an example. It comprises 738 nodes of which some were kept at constant temperatures between 0 and 300? K while the rest was kept at 325? K.

(152) FIG. 8 depicts the results for the nodes at lower temperatures for a required relative accuracy of 10.sup.?3.

(153) The relative error with respect to the correct analytical value that was calculated according to

(154) $1-\frac{W_{\mathrm{dep},\mathrm{numerical}}}{W_{\mathrm{dep},\mathrm{analytical}}}$
stayed below the threshold that was specified for the LTE. The relative error can thus be used to control the absolute error.

(155) The calculation took 87 secs needing 666 iterations.

(156) The present embodiment enables the computation of self-contamination with a computation time, which can be significantly reduced compared to known methods by applying an adaptive stepsize control of a numerical solver used for solving the basic equation for deposit. Moreover, the computed solution is exact within given accuracy requirements. The method can be implemented as part of a software tool suite for spacecraft computations and simulations.

(157) While at least one exemplary embodiment has been presented in the foregoing detailed description, it should be appreciated that a vast number of variations exist. It should also be appreciated that the exemplary embodiment or exemplary embodiments are only examples, and are not intended to limit the scope, applicability, or configuration of the embodiment in any way. Rather, the foregoing detailed description will provide those skilled in the art with a convenient road map for implementing an exemplary embodiment, it being understood that various changes may be made in the function and arrangement of elements described in an exemplary embodiment without departing from the scope of the embodiment as set forth in the appended claims and their legal equivalents.