Modeling and Thermal Analysis of a Moving Spacecraft Subject to Solar Radiation Effect

: The impact of solar radiation on spacecraft can increase the cooling load, degrade the material properties of the structure and possibly lead to catastrophic failure of their missions. In this paper, we develop a computational model to investigate the effect of the exposure to solar radiation on the thermal distribution of a spacecraft with a cylindrical shape which is traveling in low earth orbit environment. This is obtained by the energy conservation between the heat conduction among the spacecraft, the heating from the solar radiation, and the radiative heat dissipation into the surroundings while accounting for the dynamics of the space vehicle (rotational motion). The model is solved numerically using the meshless collocation point method to evaluate the temperature variations under different operating conditions. The meshless method is based on approximating the unknown ﬁeld function and their space derivatives, by using a set of nodes, sprinkled over the spatial domain of the spacecraft wall and functions with compact support. Meshless schemes bypass the use of conventional mesh conﬁgurations and require only clouds of points, without any prior knowledge on their connectivity. This would relieve the computational burden associated with mesh generation. The simulation results are found in good agreement with those reported in previously-published research works. The numerical results show that spinning the spacecraft at appropriate rates ensures low and uniform temperature distribution on the spacecraft, treated as thick-walled object of different geometries. Therefore, this would extend its lifetime and protect all on-board electronic equipment needed to accomplish its mission.


Introduction
The space environment during the spacecraft travelling time and mission continues to be a major cause of anomalies in many space missions. The most common factors of such anomalies are surface charging and discharging, solar radiation heating loads, and the existence of orbital debris in space as indicated by several published studies [1][2][3][4][5][6][7][8][9][10][11][12]. In addition, the solar load, shielding material, onboard equipment location, thermal control and air conditioning systems must be considered during the design stage of spacecraft/satellites [13][14][15][16][17][18][19]. The thermal control and vehicle dynamics play a pivotal role during the design stage of the spacecraft [19]. This is mainly due to the exposure of the spacecraft and onboard equipment to high solar loads and the penetration of heat into the spacecraft structure that affect the temperature gradient of the spacecraft body and its cooling load. The design process of any space vehicles such as satellites requires detailed thermal characteristics in addition to the dynamic analysis of the vehicle motion to guarantee no thermal stresses due to high temperature fluctuations the mathemtaical model governing the heat transfer of a rotating spacecraft. Section 3 is concerned with the description of the numerical procedure based on meshless method. In Section 4, we present and discuss the results of the numerical investigation. In Section 5, we summarize the main findings of the study and conclude with some remarks for the design of spacecraft.

Heat Transfer Modeling of a Rotating Spacecraft Under Solar Radiation
The performance requirements of a spacecraft's mission could be assessed at the earliest stages of design through the development of mathematical models and deployment of computational tools. Recent advances in computer hardware and software have enabled to reduce the computational burden associated with the numerical integration of the high-fidelity equations governing the heat transfer of the aforementioned systems [2,9,26]. Software tools can be deployed to assist the design process and provide guidelines to enable enhanced operability of spacecraft.
As a first step towards developing a good understanding of the thermal behavior of a space vehicle, being in exposure to different heat sources while traveling along the atmospheric layers, we represent this system by an infinitely long and hollow cylindrical shape rotating at constant angular speed ω. The space vehicle is expected to operate in low Earth orbit environment. The inner boundary of the cylindrical vehicle is assumed fully-insulated and the outer surface is subjected to the solar radiation and the radiative heat dissipation into the surroundings. The axis of the vehicle is placed at an angle ψ to the direction of the parallel rays of the impinging radiation of intensity f s as shown in Figure 1. The objective is to examine the temperature variations within the space vehicle under different operating conditions, including the spinning speed and the altitude.

Incident solar rays
Side view Perspective 3d view The governing equation of the temperature field is given by as follows: where κ is the thermal diffusivity, T is the temperature, r is the radius coordinate, t is the time, and θ is the angular coordinate fixed in the space vehicle. Following [9], we introduce the angular coordinate as and rewrite the governing equation as Accounting for the effect of the solar ray incidence and the radiative heat dissipation, the boundary condition at the outer vehicle surface (r = r o ) can be expressed as where K is the thermal conductivity, α is the absorptivity of the surface exposed to radiation, f s is the intensity of the solar radiation, ψ is the inclination the vehicle axis to the radiation axis, σ is the Stefan-Boltzman constant, is the emissivity of the vehicle surface, and g s is given by The inner vehicle surface (r = r i ) is insulated and the corresponding boundary condition is expressed as Expanding the nonlinear term of the outer boundary condition while assuming small temperature variation with respect to a characteristic temperature T o , introducing the following dimensionless parameters [9] we obtain the following dimensionless governing equations and boundary conditions We note that the characteristic temperature T o is given by the following expression as derived in [7] T o = α f s sin(ψ) σ π 1 4 (11) The mathematical model presented above will be used to investigate the temperature distribution over a thick-walled cylindrical space vehicle rotating at different angular speeds while exposed to solar radiation.

Theoretical Background
Meshless methods (MMs) have been recently used to solve problems arising in physics and engineering. The motivation behind meshless methods lies in relieving the burden of mesh generation. All traditional mesh-based numerical methods have difficulties on obtaining a suitable mesh. In fact, the problem of generating meshes has become more acute in recent years. Since the application of computational methods to real world problems appears to be paced by mesh generation, alleviating this bottleneck potentially impacts an enormous field of problems.
Meshless methods are based on approximating the unknown field functions and their derivatives, by using a set of nodes, "sprinkled" over the spatial domain and, functions with compact support. Meshless schemes bypass the use of a conventional mesh. The major fields of computational mechanics, such as finite element methods (FEM), finite difference methods (FDM), and finite volume methods (FVM), rely on the use of elements, interlaced grids, or finite volumes as the underlying structures upon which to discretize governing partial differential equations (PDE). On the other hand, meshless schemes only require clouds of points, without any prior knowledge on their connectivity. On this set of nodes, the spatial domain is represented and PDEs are discretized. Local clouds for each point in a domain are proximity-based subsets of the global set of points. Local clouds of points replace the more traditional forms of connectivity found in FEM, FDM, and FVM. Several formulations of the meshfree method have been formulated such as smoothed particle hydrodynamics (SPH) [31], the diffuse approximation method (DAM) [32], and the element free Galerkin method (EFG). The present study employs the discretization corrected particle strength exchange (DC PSE) [33,34]. DC PSE was formulated as an extension of the particle strength exchange (PSE) method [35]. The latter is used for the evaluation of spatial derivatives of a continuous function, discretized over scattered collocation points. A drawback from the construction procedure of the PSE operators is the introduction of an overlap condition [36], which results in a large number of particles for small kernel sizes. For the DC PSE method, the overlap condition can be relaxed by directly satisfying discrete moment conditions over collocation points. The DC PSE operators ensure that the discretization error at the collocation points does not dominate the overall order of accuracy of the approximation. Next, we introduce the DC PSE operators for strong form formulations (excluding collocation point volumes) in 2D and how to construct them. For higher dimensions and more general analysis and discussion the reader is directed to [36].
We consider the differential operator for a continuous field function where m and n are integers that determine the order of the differential operator. We define the DC PSE operator for the spatial derivative D ( m, n) f (x i ) as follows: where (x) is a local scaling parameter, η(x, (x)) is a kernel function normalized by the factor 1 (x) 2 , and N(x) is the set of points in the collocation grid within the support of the kernel function.
The objective is to construct our DC PSE operators so that as the average spacing between neighbours is decreases, h(x p ) tends to zero, around the point x p and then the operator converges to the spatial derivative D (m,n) f (x i ) with an asymptotic rate r, as expressed below where the average neighbour spacing h is defined as where N is the number of points in the support domain of x i . Additionally, we define a kernel function, η(x) and a scaling relation (x p ) which satisfy Equation (14). To achieve this, we begin by replacing the terms f (x i ) in Equation (13) with their Taylor series expansions around the point x i . This gives where are the discrete moments. The convergence behavior of Equation (16) is determined by the coefficients of the terms (x i ) j+k−m−n D (j,k) . We enforce the DC PSE operator to satisfy Equation (16) by setting the coefficients m and n of the term D (m,n) equal to one, and all other coefficients equal to zero. Finally, using the kernel function we obtain The row vector p(x) contains monomials and vector α T is the column vector of unknown coefficients. For example, if we set r = 2 and approximate the first spatial derivative in the x direction (D 1,0 ), we have l = 6 and the monomial basis is p(x) = {1, x, y, x 2 , xy, y 2 }. The discrete moment conditions can be expressed as with where k is the number of points in the support domain of the operator, l is the number of moment conditions to be satisfied and V(x) is the Vandermonde matrix, computed using the set of monomials.
The matrix E(x) is a diagonal matrix determined by the kernel's window function The matrix A, often referred to as moment matrix, gives information about the spatial distribution of the nodes x i around the centre point x. On the other hand, b determines the approximation properties of the kernels. From the mathematical point of view, the invertibility of matrix A depends entirely on that of the Vandermode matrix V.

Generation of the Point Cloud
In the context of meshless methods, generating point clouds that represent the geometry and computing spatial derivatives is a straightforward procedure. Nodes can be selected randomly or uniformly (following a uniform (Cartesian) or irregular nodal distributions) at the interior domain and on the boundary of the geometry without no prior knowledge on their connectivity. The former is easy to construct while the second requires the use of robust triangular mesh generators.
In this study, we apply an embedded Cartesian grid method to represent the complex geometry (the two cylinders). We start the procedure by setting the boundary nodes from which we will compute the average spacing h c for the Cartesian grid (Figure 2b). The Cartesian grid step size h c , is defined as the average of the inner and outer boundary nodes spacing, as given by h ave = (h in +h out ) 2 . From the Cartesian grid, we use only the nodes located in the interior of the domain, while special care is taken for the nodes located close to the boundary nodes. These nodes, often called as degenerated nodes, are the grid nodes with distance less than 0.2h c , as shown in Figure 2c. The built-in MATLAB function inpolygon has been used to define the Cartesian grid nodes that are located inside the spatial domain. Degenerated nodes are removed and not included in the point cloud, since they give rise to ill-conditioned Vandremode matrices.

Explicit Solver
The diffusion-convection problem simulating the temperature variations of the space vehicle is described by Equations (1)- (10). To implement the meshless point collocation method, it is necessary to represent the spatial domain with a set of nodes, distributed over the interior domain and on the boundary. The spatial derivatives of the dependent variable (in our case temperature) are computed by the DC PSE method. As described above, the DC PSE method computes the derivatives of temperature field in a local sense, by using the neighboring nodes. Applying the Euler explicit time integration method, one can rewrite the governing equations in the Cartesian form as follows: where T (n+1) and T (n) are the temperature values at the current and previous time steps, respectively, and u = rω cos(φ) and v = rω sin(φ) are the velocity components. We note that the time step δt of the explicit solver should be carefully selected to guarantee stable and accurate numerical solutions.
In the present application, Neumann boundary conditions are imposed in the inner and outer surfaces of the cylinder, as given by where n = (n x , n y ) is the unit normal vector pointing outward. We compute the values on the boundary nodes with Neumann boundary conditions based on values of interior nodes that belong to the support domain, as shown in Figure 3. The temperature on the inner and outer surfaces of the cylinder are computed using the Neumann boundary conditions as follows. First, we compute the weights, w x and w y , for the first spatial derivative with respect to x− and y−, respectively, based on the support nodes. The support nodes (stencil) include one boundary point (center number i) such as the stencil shown in Figure 3 and the nodes located in the interior of the geometry. The boundary values of the temperature are computed as The computational model described above is implemented on Matlab. We discretize the spatial domain using successively finer uniform Cartesian embedded grids with spacing of h = 2.5 × 10 −3 , 1.25 × 10 −3 and 6.25 × 10 −4 m to ensure a grid independent solution, resulting in 34,580, 136,972 and 545,416 nodes, respectively. The results on the two finest grids are almost indistinguishable, indicating grid independence. Therefore, in our simulations we use the grid with a spacing of h = 1.25 × 10 −3 m. The time step used in the simulations is set to dt = 10 −3 s for all angular velocity values considered while the initial values for all flow variables at the interior points were set to zero. Computations were conducted using an Intel i7 quad core processor with 16 GB RAM.

Results and Discussion
The developed model, given by Equations (8)- (10), is obtained by energy conservation between the heat conduction among the space vehicle, the heating from the solar radiation, the radiative heat dissipation into the surroundings while accounting for the dynamics of the space vehicle (rotational motion). This model is solved numerically using the meshless collocation point method (MCP) to evaluate the temperature distribution under different operating conditions. The simulation results are expected to provide baseline and guidance for the design of thermal control systems for space vehicles.
To verify the numerical predictions of the meshless collocation point method, we consider the simulated case of rotating cylinder exposed to solar radiation reported in [9]. The numerical values of the geometry and thermal parameters are presented in Table 1. The vehicle is assumed to be a black body made of natural rubber with a low conductivity. We simulate the thermal response of the space vehicle while varying the spinning speed ω from 0 to 1000 rad/h. The current simulation results are compared against those obtained from the finite difference method (FDM) by Gadalla and Wehba [9]. We plot in Figure 4 the variations of the steady-state temperature at the outer surface for different spinning speeds as obtained from the two numerical approaches: MCP (current study) and FDM (as used in [9]). The two sets of data show good agreement. This demonstrates the capability of MCP method to perform thermal analysis of moving space vehicles subjected to heat sources and dissipation. For a fixed space vehicle (ω = 0 rad/h), we observe a symmetric distribution of the temperature about the location of the radiant heat source. Once the rotational motion is initiated, this symmetry is broken. As the spinning speed increases, the maximum and minimum temperatures shift to lower and higher values, respectively. Clearly, the rotational motion of the space vehicle at higher angular speeds enables further reduction in the temperature gradients. At ω = 1000 rad/h, the temperature distribution at the outer surface is found almost insensitive to the circumferential position φ. Figure 5 depicts the variations of the temperature with the angular position φ at the internal circumferential surface located at ξ = 0.8. The current simulation results compare well with those obtained using FDM [9]. The simulations show that the rotational motion tends to reduce more the temperature gradients in the internal circumferential surfaces in comparison to those achieved at the outer surface, directly exposed to the solar radiation. At ω = 2 rad/h, a slight variation in the temperature distribution is observed.    Figure 6 shows the maximum and minimum temperatures at the outer surface (ξ = 1) for different values of the spinning speed ω. Rotating the space vehicle at higher speeds results in shrinking the bandwidth between the temperature extrema. At very high rotational speeds (ω ≥ 1000 rad/h), the maximum and minimum temperatures converge to the reference temperature T o , given by Equation (11). To gain an insight into the impact of the rotational motion on the temperature distribution over the full space vehicle, we plot in Figure 7 the temperature contours for varying spinning speeds. Again, a symmetric temperature field is obtained for a fixed space vehicle. The rotational motion breaks the symmetry and induces reduction in the temperature range. For ω higher than 10 rad/h, the temperature is observed to be varying only within a thin layer near the outer surface, otherwise constant over the internal domain. We note that these results are consistent with those obtained in [9]. The simulation results reveal the importance of incorporating the rotational motion of the spacecraft to achieve lower temperature with a uniform distribution over the entire structure. This would enable a safer operability of the spacecraft. The Earth's atmosphere comprises four main layers. Moving upward from ground level, these layers are named the troposphere, stratosphere, mesosphere, and thermosphere, each with its own specific traits. The exposure to different heat and pressure loads within these layers present a major cause of anomalies in many space missions. The present numerical analysis shows that the design of spacecraft can be then significantly enhanced in terms of thermal resistance to the exposure to heat sources by incorporating the rotational motion.
We discretize the spatial domain shown in Figure 8 using successively finer uniform Cartesian embedded grids with a spacing of h = 1.875 × 10 −3 , 1.25 × 10 −3 and 9.3755 × 10 −4 m to verify the convergence behavior of the numerical solution and obtain a grid independent temperature solutions. We note the aforementioned spacings lead to 114,490, 256,875 and 456,057 nodes, respectively. The results on the two finest grids are almost indistinguishable, indicating grid independence. As such, in the subsequent simulations, we use a grid with spacing of h = 1.25 × 10 −3 m. The time step is set to dt = 10 −6 s, while the initial values for all flow variables at the interior points were set to zero. Computations were conducted using an Intel i7 quad core processor with 16 GB RAM. We show in Figure 9 the contour plots of the temperature over the thick-walled vehicle for varying angular speeds. Keeping the vehicle at static position leads to a large temperature gradients. Activating the rotational motion results in a significant reduction in temperature variations and enables more uniform temperature distribution on the body of the space vehicle.

Conclusions
In this work, we developed a computational model based on the meshless collocation point method to simulate the temperature distribution over a rotating spacecraft under solar radiation effect. The spacecraft is treated as a thick-walled object of different geometries. The numerical results compared well with those reported in the literature. The simulations revealed that spacecraft may undergo large temperature gradients when traveling over the Earth's atmospheric layers or in low Earth orbit. We showed the usefulness of spinning the spacecraft to avoid large temperature variations that may degrade its performance. Rotating the spacecraft at an angular speed of or higher than 10 rad/h enables a quasi-uniform temperature distribution over its internal domain. The numerical study is expected to provide a baseline for thermal control of space vehicles.