© A.V. Vostrikov & E.N. Prokofeva, Published by EDP Sciences 2022

1 Introduction

Work to solve the problem of electrification of spacecraft (SC) began in the early 70s of the 20th century during the development of the geostationary orbit (De Forest & Mc. Jlwain, 1971). By charging a spacecraft, we mean the accumulation of charge on the surfaces of the spacecraft (Garrett & Whittlesey, 2000). The spacecraft is charged as a result of its interaction with the surrounding plasma and solar energetic particles; therefore, charging effects are most pronounced during geomagnetic disturbances (substorms) in the Earth’s magnetosphere (Purvis et al., 1984; Lui, 2000; Lai, 2011). The penetration of the energetic electrons into the internal components or portions of the spacecraft will be stopped and accumulated, leading to high electric fields affecting a nearby victim circuit. This is what causes an electrostatic discharge (ESD). This occurs when the electric fields resulting from the charge exceed the breakdown strength of the outer surface materials of the spacecraft. ESD, a consequence of charging, adversely affects the operation of the onboard electronic equipment of the spacecraft. Electromagnetic interference generated by ESD can reach tens of volts and cause malfunctions in the operation of various onboard systems of a spacecraft, and intense discharge currents can lead to irreversible damage to equipment elements. The main receptors of impulse interference from ESD are fragments of the onboard cable network (OCN) laid along the outer surface of the spacecraft (Garrett, 1981). Table 1 shows some examples of the effects of electrostatic discharges.

According to statistics, about 30% of the loss of satellites is due to the phenomenon of electrostatic discharge. According to statistics, in a third of the cases, the loss of a spacecraft occurs due to an electrostatic discharge that has arisen on it (Fig. 1). An electrical potential on the order of tens of volts can have detrimental effects on spacecraft components.

Figure 1 Types of spacecraft damage in geostationary orbit.

Attempts to completely eliminate the possibility of ESD by selecting materials for the outer surface of the spacecraft or active protection of the spacecraft have not yet been very successful (Tyutnev et al., 2019). It is possible to reduce the frequency and power of ESD but not completely eliminate them (Saenko et al., 2015). Therefore, it is necessary to take additional measures for the trouble-free operation of the spacecraft electronics when exposed to ESD. The spacecraft electronics must be manufactured in such a way as to guarantee its flawless operation in flight conditions. EMI levels should be assessed in advance during the spacecraft design sketch phase. Then the electronics engineers will have the necessary design information. Thus, the main goal of creating new solutions is aimed at recommendations for minimizing negative phenomena and hazards caused by space weather. To do this, before the operation of the spacecraft, it is necessary to calculate the pickups in the OCN of the spacecraft. Let us present the procedure for determining interference in the OCN of the spacecraft in Figure 2 (the stage considered in the article is marked with a red marker).

Figure 2 The procedure for determining interference in the OCN SC.

First, identify the most likely locations for electrostatic discharge. The National Aeronautics and Space Administration uses a well-known computer program for this purpose. NASA Charging Analyzer Program (NASCAP).

Secondly, it is necessary to estimate the levels of impulse noise in the electronic input circuits connected with the outer surface of the spacecraft through the fragments of the cable network. To solve this problem, we use the software “Satellite-MIEM,” (developed by Novikov et al.) is used in Russia. For the software “Sattellite-MIEM” the place of the discharge is a priori information obtained earlier. The tasks of the software “Sattellite-MIEM” are to build a picture of the spreading of currents over the surface of the spacecraft and to calculate pickups in fragments of the onboard cable network of the spacecraft. This computer system uses Windows 7 or higher. With the help of the Satellite-MIEM software, we calculated the pickups in the cables of the Spektr-R space radio telescope, which was launched on July 18, 2011, and served a warranty period in space. A preliminary result of this design phase may be the current distribution resulting from an electrostatic discharge localized at a specified point.

To do this, we used the Satellite-MIEM software (Marchenkov et al., 2008), adapted for modeling and calculating noise signals generated in cable networks and electronic inputs of a spacecraft. This approach is based on the formal representation of the model in the form of a system of differential equations describing thermal, oscillatory, and structural processes in various devices. Satellite-MIEM software for a given 3D model generates a lumped element model (LEM) of the spacecraft, which was originally used to assess the emerging potentials on the spacecraft surface due to cosmic particle fluxes (De Forest & Mc. Jlwain, 1971; Rosen, 1976a). One of the first implementations can be seen in the NASCAP software (H.B. Garret, 1981). NASCAP Software Update Confirms Concerns From NASA Spacecraft Charging (V.A. Davis et al., 2004).

The structural electrophysical model (SEM) of a spacecraft is an equivalent electrical circuit of the spacecraft surface made of lumped R, L, and C elements (Saenko et al., 2015). To construct an SEM, the spacecraft body with attached elements is divided into geometric bodies such as a cylinder, sphere, torus, plane, rod, etc.

Satellite-MIEM software synthesizes SEM based on a given polygonal 3D model of the spacecraft (Fig. 3). A polygonal model consisting of a collection of rectangles or triangles can be easily and quickly built in the 3dsMax package. Then the model is exported to the Satellite-MIEM software, where it is transformed into a surface mesh: a collection of connected nodes. Each connection (branch) is represented in the form of elements of an electric circuit as a whole, forming an equivalent electric circuit of the spacecraft surface.

Figure 3 Converting a collection of quadrangles to an SEM.

If there is information about possible places of ESD occurrence, it is possible to calculate the spreading pattern of transient currents along the body and attached elements of the spacecraft. In this case, ESDs are presented in the form of a pulsed current source with predetermined characteristics corresponding to the discharge parameters. These sources are connected to points between which ESD may occur. Next, the transient currents are calculated.

To calculate the magnitude of the interference signal in the OCN fragment, the software is used to lay the trace of this fragment along the outer surface of the spacecraft, and the value of the transformation ratio of the current flowing along the spacecraft surface is introduced into the pickup voltage in the OCN (Agapov et al., 2010). The calculated data on the levels of interference signals at the inputs of electronic units, obtained at the stage of the preliminary design of the spacecraft, are included in terms of reference (TOR) for manufacturing these units. In accordance with such TOR, the electronic unit must work without failures at this calculated level of the interference signal.

The cutoff frequency in the ESD spectrum is defined as the rise time of an electrostatic discharge pulse and is 200 MHz for a rise time of 5 ns (equivalent wavelength λ = 1.5 m). For acceptable computational accuracy, the grid cell size is equivalent to the length of the branch on the outer side of the CS. Surface should not exceed λ/10, i.e., 0.15 m. In this case, the use of LEM is permissible. With a typical satellite size of 10 m, the number of grid nodes will approach (1–2)10^5.

To calculate the spreading pattern of transient currents along the SEM of the spacecraft, the Satellite-MIEM software uses the well-known program for calculating electrical circuits PSpice from the OrCAD package. It should be noted that SPICE programs such as OrCad, Pspice, Micro-Cap, and LTspice compute R, L, and C circuits with this many nodes, using conventional multicore personal computers very slowly. Analysis of an EMS consisting of 150,000 nodes or more on modern PCs in such a program is impossible due to the enormous complexity of the computer spent on solving a system of ordinary differential equations (ODE) by the Runge–Kutta method (Rosen, 1976b). Studies have shown that the time spent solving a system of 2000 or more ODEs in PSpice software is unacceptable for calculations. We have carried out a test analysis of equivalent electrical network (EEN) models of spacecraft of various dimensions on a modern machine. The signal had a duration of 150 ns, and the sampling step was 1 ns. Figure 4 shows the dependence of the computer time consumption for experimental calculations of the spreading of transient currents from ESD to the EEN of the spacecraft on the number of nodes. It can be seen that with an increase in the dimension of the model by n times, the time for calculating the transient currents increases by approximately n² times.

Figure 4 The dependence of the calculation time on the number of nodes in the EEN model.

The following conclusion was made: the result of spreading currents on large EEN (the number of circuit nodes is 150,000) on modern personal computers takes a lot of computer time (about 80 h). Therefore, it becomes necessary to reduce the time of calculating the spreading of currents over the surface of the spacecraft.

2 Development of the method of highlighted areas

In this paper, we propose methods for the accelerated calculation of the spreading of currents from an electrostatic discharge over the surface of a spacecraft. A feature of the first method is the calculation of transient currents only in the local zone of the spacecraft circuit. In the course of the experiments carried out in the Satellite-MIEM program, an area of significant transient currents of 400 nodes around a given ESR site was identified. Insignificant transient currents can be considered currents that do not exceed 1–2% of the ESD value. Minor currents are considered due to the fact that the electromagnetic interference from such transient currents has too little effect on the onboard electronic equipment of the spacecraft (Garret H.B., 1981; Gaines E.E. et al., 1981).

For an accelerated calculation of the spreading of currents from an electrostatic discharge over the surface of a spacecraft, it is rational to use the method of selected areas developed by us.

The idea is to apply Kirchhoff’s law to every node of the equivalent circuit. The node gives five branches, four of which are connected to other nodes, and the fifth is connected to the ground through a capacitor. Consequently, the current drops sharply after several cell sizes from the ESD of node (node), for example, drops from 100 to 1 A or less. At this point, the cable voltage caused by electrostatic discharge can be disregarded.

The essence of the method consists in constructing a narrow area for calculating the spreading of currents. The calculation of transient currents will only occur in the user-limited EMS area without affecting the rest. To do this, it is proposed to convert the file received as input to the PSpice program. The idea of the method is as follows: based on the found node in which ESD occurred, select the area adjacent to the place of discharge of nodes (Fig. 5). The size of the selected area depends on the value of the error in calculating the spreading of transient currents over the spacecraft body, entered by the user. So, for example, specifying an error equal to 9%, an area equal to 441 nodes around the node where ESD occurred will be selected, and with an error equal to 1%, the selected area will consist of 1681 nodes.

Figure 5 Block diagram of the algorithm for selecting areas.

When the local area is selected, which means that when the dimension of the model is reduced by n times, the time for calculating the transient currents decreases by n³ times. The application of the method in selected areas is justified in the presence of a large spacecraft scheme. The algorithm of the method for calculating the interference in the selected area in the fragments of the spacecraft OCN is shown in Figure 6.

Figure 6 Block diagram of the algorithm of the method for calculating the interference in the selected area in the fragments of the OCN.

The result of using this method is a significant reduction in the time for calculating transient currents over the surface of the spacecraft in the presence of an error. The error of the currents was calculated on a plane in the region, the values of the transient currents in the branches of which are no more than 1–2% of the ESR value. It was found that with the expansion of the selected area, and hence, with an increase in the number of nodes around a given ESD, the error decreases (Fig. 7). When 1681 nodes are selected around the discharge point, the error becomes insignificant (equal to 1%) and is quite acceptable for real calculations. The main advantage of this approach is the ability to read large-scale RLC circuits on modern computers in a short time.

Figure 7 The decrease in the error with an increase in the number of nodes in the selected area.

This method is based on a heuristic approach and is used when the accuracy of calculations is 3–5%. The main idea is to apply Kirchhoff’s law to every node of an equivalent electrical circuit. The node gives five branches, four of which are connected to other nodes, and the fifth is connected to the ground through a capacitor. Consequently, the current drops sharply after a few cell sizes from the ESD node, for example, drops from 100 to 1 A or less. At this point, the cable voltage caused by electrostatic discharge can be disregarded.

The main disadvantage of this method is that third-party software (PSpice) is required for its implementation. In this work, we propose another method that does not have this drawback.

2.1 Development of a macromodel of a spacecraft based on explicit and implicit Euler methods

Let us present the problem before us in the language of mathematics. SEM of the spacecraft can be formed in an extended homogeneous coordinate basis and written in the form of a system of linear ordinary differential equations

$$C\frac{\mathrm{d}}{\mathrm{d}t}\bar{X}\left(t\right)+G\bar{X}\left(t\right)=\bar{Y}\left(t\right),\hspace{1em}\bar{X}\left(0\right)={\bar{X}}_0$$(1)

where are C, G – numeric order matrices (n × n), moreover C – capacitance and inductance matrix, G – resistance matrix, $\bar{X}\left(t\right)$ – the vector of the sought phase variables (voltages at all nodes of the circuit and currents flowing through inductive elements), $\bar{Y}\left(t\right)$ – vector of inputs.

It is necessary to carry out the solution of the system of equations at the moment of time with minimal time expenditures t*, i.e., calculate the numerical vector $\bar{X}\left(t^{\mathrm{*}}\right)$.

To calculate the transient process of the system of equations, it is necessary to solve the system (1). It is possible to use various methods for solving the problem, numerical, numerical-analytical, based on macromodeling, etc.

Explicit Euler’s Method for solving a system of equations

$$\frac{\mathrm{d}}{\mathrm{d}t}\bar{X}\left(t\right)=F\left[\bar{X}\left(t\right)\right], \hspace{1em}\bar{X}\left(0\right)={\bar{X}}_0,$$

can be obtained by expanding the solution into a truncated Taylor series in the neighborhood h_i any point t_i

$$\bar{X}\left(t_i+h_i\right)=\bar{X}\left(t_i\right)+\frac{\mathrm{d}}{\mathrm{d}t}\bar{X}\left(t_i\right)h_i+o\left(h_i^2\right).$$

From where

$$\dot{\overline{X}}\left(t_i\right)=\left(\overline{X}\left(t_{i+1}\right)-\overline{X}\left(t_i\right)\right)h^{-1}.$$(2)

The implicit Euler method is also derived from the expansion of the solution into a truncated Taylor series

$$\bar{X}\left(t_i-h_i\right)=\bar{X}\left(t_i\right)-\frac{\mathrm{d}}{\mathrm{d}t}\bar{X}\left(t_i\right)h_i+o\left(h_i^2\right).$$

From where

$$\dot{\overline{X}}\left(t_{i+1}\right)=\left(\overline{X}\left(t_{i+1}\right)-\overline{X}\left(t_i\right)\right)h^{-1}.$$(3)

Let us show the idea of constructing a macromodel using Euler’s formulas and assuming that the SEM of the spacecraft is devoid of any specificity; that is, the model will be presented in a general form. Let us represent the mathematical model of the spacecraft SEM in block form

$$\left[\begin{array}{cc}{C}_{11}& C_{12}\\ C_{21}& C_{12}\left(Q\right)\end{array}\right]\left[\begin{array}{c}{\dot{\overline{X}}}_1\\ {\dot{\overline{X}}}_2\end{array}\right]+\left[\begin{array}{cc}{G}_{11}& G_{12}\\ G_{21}& G_{12}\left(Q\right)\end{array}\right]\left[\begin{array}{c}{\overline{X}}_1\\ {\overline{X}}_2\end{array}\right]=\left[\begin{array}{c}{\overline{Y}}_1\\ {\overline{Y}}_2\end{array}\right],\hspace{1em}\overline{X}\left(0\right)={\overline{X}}_0$$(4)

where Q are variable parameters, matrices G and C are nondegenerate. After specifying the discharge points, the vector $\overline{Y}$ will have only two nonzero elements in the subvector ${\overline{Y}}_2$.

The macromodeling algorithm consists of the preliminary transformation of problems (2) and (3) into some reduced spacecraft model, the number of equations in which will be equal to m ≪ n, which reduces the complexity of the calculations.

Substituting expressions for derivatives (2) and (3) into the model (4), we obtain a system of 4 equations

$$(G_{11}-C_{11}{h}^{-1})({\overline{X}}_1{)}_i+(G_{12}-C_{12}{h}^{-1}\left)\right({\overline{X}}_2{)}_i+C_{11}{h}^{-1}({\overline{X}}_1{)}_{i+1}+C_{12}{h}^{-1}({\overline{X}}_2{)}_{i+1}={\overline{Y}}_1,$$(4.1)

$$(G_{21}-C_{21}{h}^{-1})({\overline{X}}_1{)}_i+(G_{22}\left(Q\right)-C_{22}\left(Q\right)h^{-1}\left)\right({\overline{X}}_2{)}_i+C_{21}{h}^{-1}({\overline{X}}_1{)}_{i+1}+C_{22}(Q\left)h^{-1}\right({\overline{X}}_2{)}_{i+1}={\overline{Y}}_2,$$(4.2)

$$(C_{11}{h}^{-1}+G_{11})({\overline{X}}_1{)}_{i+1}+(C_{12}{h}^{-1}+G_{12}\left)\right({\overline{X}}_2{)}_{i+1}-C_{11}{h}^{-1}({\overline{X}}_1{)}_i-C_{12}{h}^{-1}({\overline{X}}_2{)}_i={\overline{Y}}_1,$$(4.3)

$$(G_{21}+C_{21}{h}^{-1})({\overline{X}}_1{)}_{i+1}+(G_{22}\left(Q\right)+C_{22}\left(Q\right)h^{-1}\left)\right({\overline{X}}_2{)}_{i+1}-C_{21}{h}^{-1}({\overline{X}}_1{)}_i-C_{22}(Q\left)h^{-1}\right({\overline{X}}_2{)}_i={\overline{Y}}_2.$$(4.4)

Let us carry out a series of substitutions and transformations. The aim of the transformation is to reduce the system (4). Substitute (4.2) into (4.4), then substitute (4.2) into (4.1), then substitute the result of the second substitution into the result of the first substitution. The result will be a new macro model or computational scheme (5):

$$\begin{array}{c}\begin{array}{cc}({\overline{X}}_2{)}_{i+1}=& \{\left(\left(G_{22}\left(Q\right)+C_{22}\left(Q\right)h^{-1}\right)+C_{22}\left(Q\right)h^{-1}{\left(C_{21}{h}^{-1}\right)}^{-1}a\right)-\left[{\left(C_{21}{h}^{-1}\right)}^{-1}\mathrm{da}+C_{21}{h}^{-1}\right]\end{array}\\ \times {{\left[c-\left(C_{11}{h}^{-1}\right)d{\left(C_{21}{h}^{-1}\right)}^{-1}\right]}^{-1}\left[{(C}_{11}{h}^{-1}{C}_{22}\left(Q\right)h^{-1}\right]{\left(C_{21}{h}^{-1}\right)}^{-1})-C_{12}{h}^{-1}]\}}^{-1}\\ \begin{array}{c}\times (\{\left[a{\left(C_{21}{h}^{-1}\right)}^{-1}b+C_{22}\left(Q\right)h^{-1}\right]+\left[{\left(C_{21}{h}^{-1}\right)}^{-1}\mathrm{da}+C_{21}{h}^{-1}\right]{\left[{c-\left(C_{11}{h}^{-1}\right)}^{-1}d+{\left(C_{21}{h}^{-1}\right)}^{-1}\right]}^{-1}\\ \times \left[-\left(G_{12}-C_{12}{h}^{-1}\right)+C_{11}{h}^{-1}b{\left(C_{21}{h}^{-1}\right)}^{-1}\right]\}({\overline{X}}_2{)}_i+{\overline{Y}}_2-a{\left(C_{21}{h}^{-1}\right)}^{-1}{\overline{Y}}_2+\left[{\left(C_{21}{h}^{-1}\right)}^{-1}\mathrm{da}+C_{21}{h}^{-1}\right]\\ \times {\left[c-\left(C_{11}{h}^{-1}\right)d{\left(C_{21}{h}^{-1}\right)}^{-1}\right]}^{-1}\left[-{\overline{Y}}_1-C_{11}{h}^{-1}{\overline{Y}}_2{\left(C_{21}{h}^{-1}\right)}^{-1}\right]).\end{array}\end{array}$$(5)

Based on the fact that the subvector ${\bar{Y}}_1$ is zero, then after multiplying all matrices in (5), the final expression will have the form

$$({\bar{X}}_2{)}_{i+1}=(A_1+A_{2}h\left)\right({\overline{X}}_2{)}_i+(A_3+A_{4}h){\overline{Y}}_2+{\overline{Y}}_2.$$(6)

where are A₁, A₂, A₃, A₄ – numeric (m × m) – matrices.

But the model (1) considered by the authors has a number of peculiarities. Let us note the specifics of the scheme model.

1. The matrix is diagonal and degenerate. The diagonal consists of two groups of coefficients, and the coefficients within each group are the same.

2. Matrix - nondegenerate, symmetric, and sparse matrix.

3. Vector $\bar{Y}\left(t\right)$ contains only one or two nonzero coefficients of the form y_it, where the subscript “i” can be any number from 1 to n.

We will assume that the original problem has the following form:

$$\frac{\mathrm{d}}{\mathrm{dt}}\bar{X}\left(t\right)+A\bar{X}\left(t\right)=\bar{Y}\left(t\right), \hspace{1em}\bar{X}\left(0\right)={\bar{X}}_0.$$(7)

Taking into account the specific form of the system of equations (7), we obtain the following formula of the Euler method

$${\overline{X}}_{i+1}\approx {\overline{X}}_i+\left(-A{\overline{X}}_i+{\overline{Y}}_i\right)h_i.$$(8)

Taking into account the specific form of the system of equations (6), we obtain the following formula of the Euler method

$$\left(E+h_{i}A\right){\overline{X}}_{i+1}\approx {\overline{X}}_i+{\overline{Y}}_i{h}_i.$$(9)

The algorithm of macromodeling developed by us consists of the preliminary transformation of problems (8) and (9) into some reduced spacecraft model, the number of equations which will be equal, which reduces the complexity of calculations. A feature of the method is the calculation of transient currents only in the local zone of the spacecraft circuit. As previously indicated, this action is correct due to too little influence on the spacecraft radio-electronic equipment of the calculated values of electromagnetic interference outside the selected area equal to 400 nodes. Therefore, for the accelerated calculation of the spreading of currents from ESD over the surface of the spacecraft, it is rational to use the method proposed in this work. To do this, we write both types of computational processes in block form:

$$\left[\begin{array}{c}{\bar{U}}_1\\ {\bar{U}}_2\end{array}\right]=\left[\begin{array}{c}{\bar{V}}_1\\ {\bar{V}}_2\end{array}\right]-\left[\begin{array}{cc}{A}_{11}& A_{12}\\ A_{21}& A_{22}\end{array}\right]\left[\begin{array}{c}{\bar{V}}_1\\ {\bar{V}}_2\end{array}\right]h+\left[\begin{array}{c}{\bar{Y}}_1\\ {\bar{Y}}_2\end{array}\right]h,$$(10)

$$\left[\begin{array}{cc}{E}_{11}+h{A}_{11}& h{A}_{12}\\ h{A}_{21}& E_{22}+h{A}_{22}\end{array}\right]\left[\begin{array}{c}{\bar{U}}_1\\ {\bar{U}}_2\end{array}\right]=\left[\begin{array}{c}{\bar{V}}_1\\ {\bar{V}}_2\end{array}\right]+\left[\begin{array}{c}{\bar{Y}}_1\\ {\bar{Y}}_2\end{array}\right]h,$$(11)

where ${\left[{\bar{U}}_1\hspace{0.5em}{\bar{U}}_2\right]}^T={\bar{X}}_{i+1}^T$, ${\left[{\bar{V}}_1\hspace{0.5em}{\bar{V}}_2\right]}^T={\bar{X}}_i^T$. The specificity is that block matrix A is sparse and has a high order. In this regard, it is proposed, before the second stage of model reduction, to transform the matrix A to a triangular form with a border using the method of determining quantities Sangiovanni-Vincentelli & Bickart (1979). After transforming the matrix A we get the blocks: A₁₁ is a non-defective lower triangular matrix of order (n*n), A₁₂ – matrix (n*m), A₂₁ – matrix (m*n), A₂₂ – matrix (m*m).

We will assume that the subvectors ${\bar{U}}_2$ and ${\bar{V}}_2$ contain on m ≪ n the desired solution coefficients. It is necessary, on the basis of expressions (10) and (11), to formulate a new computational scheme from which the subvectors ${\bar{U}}_1$ and ${\bar{V}}_1$.

For this, we write (10) and (11) in the form of the following subsystems of equations

$${\bar{U}}_1=(E_{11}-h{A}_{11}){\bar{V}}_1-h{A}_{12}{\bar{V}}_2+h{\bar{Y}}_{1,}$$(10.1)

$${\bar{U}}_2=-h{A}_{21}{\bar{V}}_1+\left(E_{22}-h{A}_{22}\right){\bar{V}}_2+h{\bar{Y}}_2,$$(10.2)

$$\left(E_{11}+h{A}_{11}\right){\bar{U}}_1+h{A}_{12}{\bar{U}}_2={\bar{V}}_1+h{\bar{Y}}_1,$$(11.1)

$$h{A}_{11}{\bar{U}}_1+\left(E_{22}+h{A}_{12}\right){\bar{U}}_2={\bar{V}}_2+h{\bar{Y}}_2.$$(11.2)

Expressions (10.1), (10.2), (11.1), and (11.2) will be considered as an analog of a system of four equations with four unknowns ${\bar{U}}_1$, ${\bar{U}}_2$, ${\bar{V}}_1$, ${\bar{V}}_2$, based, generally speaking, on computational processes of different types (explicit and implicit).

Substitute (10.1) in (11.1):

$${\bar{V}}_1\hspace{0.5em}=\hspace{0.5em}-h^{-1}(A_{11}^2{)}^{-1}\left[\left(E_{11}+h{A}_{11}\right)A_{12}{\bar{V}}_2-A_{12}{\bar{U}}_2-h{A}_{11}{\bar{Y}}_1\right].$$(12)

Now we substitute the obtained expression (12) into (10.2).

$$[E_{22}+A_{21}(A_{11}^2{)}^{-1}{A}_{12}]{\bar{U}}_2=[A_{21}\left(A_{11}^2{)}^{-1}\right(E_{11}+h{A}_{11})A_{12}+(E_{22}-h{A}_{22}\left)\right]{\bar{V}}_2+h[-A_{21}\left(A_{11}^2{)}^{-1}{A}_{11}{\bar{Y}}_1+{\bar{Y}}_2\right].$$(13)

After multiplying all matrices in (12), the final expression will have the form

$$B_1({\bar{X}}_2{)}_{i+1}=[B_{2}h+B_3\left]\right({\bar{X}}_2{)}_i+h\left[B_4\right({\bar{Y}}_1{)}_i+\left({\bar{Y}}_2{)}_i\right],$$(14)

where B₁, B₂, B₃ – numeric (m × m) – matrices, B₄ – numeric (m × n) – matrix. In practical problems, the vector $\bar{Y}\left(t\right)$ contains only a few nonzero coefficients. Therefore, in practice, expression (14) will have a simpler form

$$B_1({\bar{X}}_2{)}_{i+1}=[B_{2}h+B_3\left]\right({\bar{X}}_2{)}_i+h\left({\bar{Y}}_2{)}_i\right].$$(15)

By construction, the matrix B₁ nondegenerate and independent of the parameter h. Therefore, the inverse matrix can be calculated from it, after multiplying by which expressions (15), we obtain a computational scheme in the form

$$({\bar{X}}_2{)}_{i+1}=[C_{1}h+C_2\left]\right({\bar{X}}_2{)}_i+h{C}_3\left({\bar{Y}}_2{)}_i\right].$$(16)

As a result, a macromodel was built, and a new computational scheme was obtained, in which only m ≪ n vector coefficients $\bar{X}\left(t\right)$.

The complexity T of the process of constructing a computational scheme (13):

$$T=4n^2+2m^2+2m{n}^2+m^{2}n+n \mathrm{RMO}.$$

The complexity of calculations according to the scheme (16) with a changed step will be 3m² real multiplicative operations (RMO), and with a constant step – m² RMO. If n = 1,50,000 and for dense matrices, the computation speed will increase 7,031,250,000 times as compared to the implicit Euler method and 46,875 times as compared to the explicit Euler’s method.

The application of this approach to the macromodeling of the EES of the spacecraft is justified in the presence of a spacecraft of a large dimension, with a significant time interval t^*, that is when it is required to calculate the numerical vector $\bar{X}\left(t^{*}\right)$, having taken a large number of steps.

According to the proposed method, test calculations were carried out and compared with calculations performed by traditional methods. A comparison of the calculation results showed their complete adequacy.

2.2 Comparison of time costs

The practice has shown that calculating large mathematical models according to the computational scheme is 2–3 orders of magnitude faster than the methods used in specialized software. The time for analyzing an EPS, consisting of 1000 nodes, by different methods with a step h = 0.001 s and the number of steps 1000 is given in Table 2, and the time for constructing a reduced computational scheme is given in Table 2. The experiments were carried out on a PC with a dual-core processor (clock frequency 1.8 GHz per core) and 2 GB of RAM. In this case, the results of the calculation by the reduced computational scheme differ from the results of the implicit method ≪ 1%. Note that the time of constructing the macromodel is comparable to the time of a single analysis of the circuit by the implicit Euler method (Table 3). Thus, the use of the proposed reduced scheme is justified in case of its repeated use.

2.3 Practical example

The design scheme was tested on large problems, but it does not seem rational to illustrate these examples due to a large amount of data. To confirm the results of the operability of the given computational circuit, we will conduct an experiment without any impact on the investigated electrical circuit in the laboratory (room conditions). To do this, on the model, we will build a test equivalent electrical circuit consisting of resistors with a nominal value of 100 Ω, capacitors of 100 pF, and intactivities of 1 μH each, as a result, we get a circuit consisting of 23 nodes. We will connect a signal generator to the 15th node of the assembled circuits, to which we will supply a pulsed periodic signal (meander) as a disturbing effect; we will also connect an oscilloscope to fix the potentials in the circuit nodes. The meander has the following characteristics (Fig. 8): Von [V] = 5, T_rise [s] = 0.3e–6, T_fall [s] = 0.3e–6, T_on [s] = 0.3e–6, Tperiod [s] = 1e–6. In Figure 8, the oscilloscope captures signals at nodes 15 and 23 of the circuit. Figure 9 shows the applied signal (meander), while the abscissa shows the time with a grid step of 0.2e–6 s, and the y-axis shows the voltage with a grid step of 0.5 V.

Figure 8 Signal captured at node 15.

Figure 9 EES containing 23 nodes.

In the future, we will obtain a similar result with the same electrical circuit in the LTSpice program, the adequacy of which in the field of computer modeling of digital and analog and digital electrical circuits has been confirmed by world experience (Fig. 9). Figure 9 shows the applied signal (meander), while the abscissa shows the time with a grid step of 0.2e–6 s, and the y-axis shows the voltage with a grid step of 0.5 V. Figure 9 below shows the diagram assembled in the LTSpice software, corresponding to the layout. Applying the tester to each of the 23 nodes, we get 23 marked signals. The signals obtained in the LTSpice software completely coincide with the experiment.

Next, we will carry out a similar calculation with the same electrical circuit according to the proposed reduction method. Let us evaluate the results obtained after analyzing the scheme with various computational schemes: the explicit Euler method, the implicit Euler method, and the proposed reduced method. Let the perturbation be associated with the 15th node, and the desired local area for calculation according to the proposed calculation scheme is located at nodes: 15, 16, and 17. Let us estimate the values of the potentials at the nodes at one point in time. Moment 400 ns by various calculation methods. (Table 4).

The result of calculating the circuit using the LTSpice software coincides with the result of calculating the implicit Euler method. The analysis results obtained by the implicit Euler method and the reduced computational scheme differ by 0.0006%. At the stage of the preliminary design of the pattern of current spreading through the body of the spacecraft, such an error is more than acceptable (Table 4).

3 Conclusion

1. The research questions or problem defined. We considered the process of charging spacecraft and the electrostatic discharge associated with it. Then we examined the consequences of an electrostatic discharge on the onboard radio-electronic equipment of a spacecraft and presented statistics according to which up to a third of spacecraft fail due to charging. Therefore, before launching a spacecraft into outer space, it is necessary to carry out calculations of pickups delivered to the input of electronic units. Next, we discussed the existing software for calculating the current spreading pattern. In this case, the calculation takes too much time on a conventional personal computer. The urgent problem is to develop a new accelerated method for constructing a pattern of current spreading over the spacecraft body.

2. Data and assumptions. To solve the problem posed, it is proposed to reduce the system of equations describing the process of current spreading based on the well-known explicit Euler method and implicit Euler method. The computer time of system transformation should be no more than the time of analysis by traditional methods. The calculation according to the proposed scheme should be faster than traditional methods.

3. Rationale of the proposed method. To justify the proposed method, we write the electrical circuit model in the extended homogeneous coordinate basis (1). Next, we write the computational process in the form of explicit and implicit Euler methods (4). After the reduction of computational processes, a macromodel or a new computational scheme (6) was obtained. Taking into account the peculiarities, the new computational scheme takes the form (13) or (16).

4. Results. To confirm the performance of the new computational scheme (16), we conducted a laboratory experiment. To do this, we assembled a working circuit from resistors, inductors, capacitors, and a source of disturbing influences. The scheme consisted of 23 nodes. A meander was applied as a disturbing influence. The result of the signal passing through the circuit is shown in Figure 8. Further, to check the adequacy of the result, the same circuit was built in the LTSpice simulator. We obtained the values of the currents in the nodes of the circuit at a time instant of 400 ns, then the same result with an error of 0.0006% was obtained using a new computational scheme.

5. Assessment based on author’s validation. The calculation time for the new computational scheme is shown in Table 2. For the 23-node scheme, the implicit Euler method took 11.4 s; the explicit Euler method took 1057.01 s, and the new computational scheme 0.063 s. The time for calculating unknown quantities in a user-specified local area has been reduced by several orders of magnitude compared to the calculation of unknown full models. It took 851.4 s to build a new computational circuit. The calculation error according to the proposed computational scheme was 0.0006%; such an error is more than acceptable. Thus, the use of a new computational scheme is justified if it is used repeatedly.

Acknowledgments

This article is an output of a research project implemented as part of the Basic Research Program at the National Research University Higher School of Economics (HSE University). The editor thanks two anonymous reviewers for their assistance in evaluating this paper.

References

All Tables

All Figures

	Figure 1 Types of spacecraft damage in geostationary orbit.
In the text

	Figure 2 The procedure for determining interference in the OCN SC.
In the text

	Figure 3 Converting a collection of quadrangles to an SEM.
In the text

	Figure 4 The dependence of the calculation time on the number of nodes in the EEN model.
In the text

	Figure 5 Block diagram of the algorithm for selecting areas.
In the text

	Figure 6 Block diagram of the algorithm of the method for calculating the interference in the selected area in the fragments of the OCN.
In the text

	Figure 7 The decrease in the error with an increase in the number of nodes in the selected area.
In the text

	Figure 8 Signal captured at node 15.
In the text

	Figure 9 EES containing 23 nodes.
In the text