Nonlinear analysis of thermal behavior for a small satellite in Low Earth Orbit using many-Node model

66 SCIENCE & TECHNOLOGY DEVELOPMENT, Vol 20, No.K2- 2017  Abstract— In this paper, nonlinear thermal responses of a small satellite in Low Earth Orbit (LEO) are analyzed using many-node model. The main elements of primary structure of the satellite include six rectangular cover plates and a solar array linking with satellite's body. These elements can be modeled as different lumped thermal nodes. We use an eight-node model for estimating temperatures at nodal elements i.e. six nodes for cover plates, and two nodes for front and rear surfaces of the solar array. The nodes absorb three major heat energy sources from the space environment consisting of solar irradiation, Earth’s albedo and infrared radiation. The established system of thermal balance equations for nodes is nonlinear and is solved by a numerical algorithm. For simulation purpose, it is assumed that the satellite always remains Earth-pointing attitude during motion. Temperature evolutions of nodes in time are explored in details. The obtained results show that the predictive temperature values of nodes are within the allowable temperature limit range of the satellite. Index Terms—small satellite, Earth- pointing, thermal response, temperature limits. 1 INTRODUCTION hemal analysis is one of the most important tasks in processes of designing, manufacturing and launching a satellite. It guarantees that all kind Manuscript Received on July 13th, 2016. Manuscript Revised December 06th, 2016. This research is funded by Vietnam National Foundation for Science and Technology Development (NAFOSTED) under grant number: "107.04-2015.36". Pham Ngoc Chung, Faculty of Basic Sciences, University of Mining and Geoology, Duc Thang, Bac Tu Liem Dist., Hanoi, Vietnam (e-mail: chunghumg86@gmail.com). Nguyen Dong Anh, Institute of Mechanics, Vietnam Academy of Science and Technology, 264 Doi Can Str., Ba Dinh Dist., Hanoi, Vietnam (e-mail: ndanh@imech.vast.vn). Nguyen Nhu Hieu, Institute of Mechanics, Vietnam Academy of Science and Technology, 264 Doi Can Str., Ba Dinh Dist., Hanoi, Vietnam (e-mail: nhuhieu1412@gmail.com). of equipment of satellite will work within allowable temperature limits [1,2,3]. The prediction of temperature fluctuations under the effect of space environment aims to design thermal for satellites in the early stage of space missions. One can use a single-node, two-node or many-node model for estimating temperatures of satellite. For simple thermal models such as single-node or two-node model, analytical methods can be used, for example, the Fourier analysis method [4], techniques of linearization method [5,6]. In the work by Grande et al. [7], they utilize a technique for linearizing nonlinear terms relating to the thermal radiation of a two-node model. Their obtained linearized system takes the form that is similar to conventional mechanical system subjected to periodic excitations and easy to solve analytically. By employing perturbation and numerical methods, Gaite et al. [8] showed that the temperature response of the satellite model approaches an attracting limit cycle. In 2012, Gaite et al. [9] studied a simple single-node model for a small satellite orbiting around a solar system planet. More recently, Anh et al. [10] have extended an equivalent linearization technique based on the dual replacement concept for finding approximate thermal responses of a single-node satellite model in the Low Earth Orbit. Some other analytical techniques for analyzing satellite thermal can be found in work by Gaite [11]. In fact, due to the complexity of the geometrical model of satellites, thermal equations of satellites are discretized into many-node in which each node is characterized by a temperature at any time. If using analytical methods, solving such a many-node system is not easy. The current study is devoted to the use of a numerical method to analyze thermal behavior for a small rectangular parallelepiped satellite in Low Earth Orbit using eight-node model. The obtained result shows temperature evolutions of nodes in Nonlinear analysis of thermal behavior for a small satellite in Low Earth Orbit using many-node model Pham Ngoc Chung, Nguyen Dong Anh and Nguyen Nhu Hieu T TẠP CHÍ PHÁT TRIỂN KH&CN, TẬP 20, SỐ K2-2017 67  time  and  temperature  values  are  within  allowable  limit.  We  also  examine  the  effect  of  the  solar  absorptivity  and  the  emissivity  on  the  satellite’s  thermal responses.  It  is  noticed  that  modeling  of  thermal  loadings  acting  on  satellite  is  not  available  for  almost  thermal analyses using many-node models. This  is  because  of  the  complexity of  modeling of  thermal  inputs  and  of  the  dependence  on  orbital  configuration,  motion  characteristics  and  mission  of  satellites.  Thermal  loadings  are  usually  simulated  numerically  and  integrated  in  commercial softwares [1]. In this research, we have  constructed  an  appropriate  modeling  for  thermal  loadings  in  the  case  of  satellite's  attitude  being  Earth-pointing. Thermal characteristics of nodes are  analyzed  based  on  these  input  thermal  loadings.  Our  modeling  is  an  extension  of  a  previous  paper  on  thermal  radiation  analysis  for  solar  arrays  of  a  small satellite in Low Earth Orbit [15].  2 THERMAL BALANCE EQUATIONS FOR EIGHT- NODE MODEL OF SMALL SATELLITE   2.1 A satellite model and its orbit  In Fig. 1  illustrates a small satellite moving in a  Low  Earth  Orbit  (LEO)  at  altitude  of  680  km.  Satellite's orbit  is Sun-synchronous and orbit plane  is  parallel  with  solar  rays.  The  orbital  period  is  5902.25 seconds and the eclipse duration is 2121.2  seconds.  For  simulation,  we  suppose  that  the  satellite  always  remains  Earth-pointing  attitude  during motion.  Simulation time starts at beginning  of eclipse.  The  small  satellite  is  modeled  as  illustrated  in  Fig.  2.  The  satellite  includes  a  body  of  size  B B BL W H    and  a  solar  array  of  size  A AL W   .  The  distance  from  the  solar  array  to  satellite body   is  AB . Assume that the solar array is perpendicular  to  a  side  of  satellite  body.  In  fact,  the  solar  array  may  be  placed  at  different  positions  on  the  body  depending on  the configuration and mission of  the  satellite. Because the satellite thermal calculation is  quite complex, the above model is a simplified one  and  will  be  a  basis  for  the  more  complex  satellite  model.  The  body  of  the  satellite  is  made  from  composite  materials  which  have  specific  material  and  geometric  parameters.  Here,  we  suppose  that  material  is homogeneous, which can be considered  as  a  result  after  material  homogenization.  The  absorptivity and emissivity coefficients of the body  material  are  B   and  B ,  respectively.  The  solar  array  is  composed  of  many  different  materials.  It  includes  two  surfaces:  front  surface  (surface  8)  contains  solar  cells  absorbed  energy  directly  from  Sun's  rays;  the  absorptivity  of  the  front  surface  is  F  whereas emissivity coefficient is denoted to be  F ;  and  rear  surface  (surface  7)  is  coated  by  a  material layer with absorptivity  R , and emissivity  R . The cover plates 1, 2, 3, 4, 5, 6 are numbered  as shown in Fig.2. Plates 1 and 3 are opposite each  other,  in which plate 1  is closer  to  the  solar array.  Plates  2  and  4  are  parallel  each  other  and  perpendicular  to plates 1 and 3. Plates 5 and 6 are  upper and base plates of the satellite, respectively.  Figure 1.  Earth-pointing attitude of the satellite  A local coordinate is associated with satellite and  its origin  is at  the  intersection of  three planes 2, 3  and  5  (see  Fig.  2).  The  axis  x   is  along  the  intersection between  two planes 3  and 5,    the axis  y   is  along  the  intersection  between  two  planes  2  and 5,  the axis  z  is along the intersection between  two  planes  2  and  3.  Numbers  1  to  8  indicate  that  the  satellite  structure  is  separated  into  eight-node  with thermal characteristics assigned to each node.   Figure  2.  A  model  of  a  small  satellite  with  size  B B BL W H  and nodes numbered from 1 to 8   68          SCIENCE & TECHNOLOGY DEVELOPMENT, Vol 20, No.K2- 2017 Figure 3. A orbital model for thermal calculations  In  reality,  each  of  plates  numbered  by  indices  from  1  to  6  has  two  sides:  inside  and  outside  surfaces.  This  way  of  numbering  is  because  satellite  plates  have  been  assumed  to  be  homogeneous  material  and  therefore  the  material  properties of the inside and outside are the same. It  means  that we can consider  the  inside and outside  surfaces  as  a  representative  surface  with  an  assigned  index.  However,  the  satellite  array  is  modeled and numbered with  two separated  indices  (7 and 8). This  is because  the physical property of  the  front  surface  that  contains  solar  cells  is  very  different from the rear surface. 2.2 Thermal sources to satellite  As  mentioned  above,  there  are  three  main  heat  sources from space environment consisting of solar  irradiation,  Earth’s  albedo  and  infrared  radiation  that  affect  to  the  thermal  behavior  of  the  satellite.  From the Earth-pointing attitude of the satellite (see  Fig.  1  or  Fig.  3),  we  can  obtain  thermal  fluxes  acting  to  nodes.  The  order  of  nodes  in  thermal  calculation  is  shown  in Tab. 1.  Parameter  data  for  our calculations are given by Tab. 2. For detail, the  values  and  physical  characteristics  of  these  parameters can be seen in [1- 3, 10, 12].  During  motion,  only  six  surfaces  receive  the  thermal  loadings  from  the  space  environment  are  X+,  X-,  Z+,  Z-,  front  and  rear  surfaces;  also  for  other  two  sides  Y+  and  Y-,  the  applied  thermal  loadings are considered to equal zero.   TABLE 1. THE ORDER OF NODES IN THE THERMAL  CALCULATION  TABLE 2  MATERIAL PARAMETERS FOR THERMAL CALCULATION  System parameters Values Length of the body  BL  (m)   0.5  Width of the body  BL  (m)  0.5  Height of the body  BH  (m)  0.5  Mass density of body  plates B  (kg/m 3)  158.9  Specific heat capacity of body plates  BpC  (J/kgK)  883.70  Thickness of body plates  B  (m)  0.02  Material conductivity of the body  B  (W/mK)  5.39  Emissivity of the body material  B    0.82  Absorbsivity of the body material  B    0.65  Length of  the solar array  AL  (m)  0.7  Width of the solar array  AW  (m)   0.5  Material conductivity of the solar array A   (W/mK)  2.79  Thickness of the solar array A  (m)   0.03  Mass density of the solar array  A  (kg/m 3)  111.7  Specific heat capacity of the solar array  ApC   (J/kgK)  844.40  Emissivity of front surface F    0.82  Absorbsivity of  front surface F    0.69  Emissivity of rear surface  R    0.872  Absorbsivity of  rear surface  R    0.265  The distance from array to body AB  (m)   0.02  Solar constant  sG  (Wm -2)  1440  Earth albedo coefficient a    0.65  Earth black-body equivalent temperature  eT  (K)   259  Attitude of the orbit h (km)  680  Radius of the Earth  eR  (km)   6400  Orbital period  orbP  (s)  5902.25  Eclipse duration  ecP  (s)  2121.2  2.2.1  Surface  X+  For the surface X+, only solar thermal flux  ,s Xq    is present. It is determined as follows  TẠP CHÍ PHÁT TRIỂN KH&CN, TẬP 20, SỐ K2-2017 69     s,X 0 0 if 0, cos if , 2 0 if , 2 ec ec orb B s ec ec ec orb orb t P P P q G t P P P t P             = - -             (1)  where  sG  is the mean solar radiation,  ecP  and  orbP   are  the  eclipse  and  orbital  periods,  respectively.  The basis for formulating this expression can view  in [1-3]. In (1), we have denoted  t = ,  2 orbP   = ,   0 arccos e e R R h    =     ,  02ec  = -       (2)                                   The  solar  flux on  the  surface  X    is a periodic  function  of  time  t  as  presented  in  Fig.  4  in  an  orbital period.  Figure  4.  Solar  flux  on  surface  X+  with  parameters  given   in Tab. 2  Figure  5.  Solar  flux  on  surface  X-  with  parameters  given   in Tab. 2  2.2.2  Surface  X-  For the surface X-,  it absorbs only solar thermal  flux  ,s Xq - and is calculated as follows       0 s,X 0 0 0 0 if 0, 2 1 cos 2 if , 2 2 2 1 cos if 2 , 2 ec orb ec orb B s orb B s ec orb orb P P t P P q G t P G t P P              -            = - -  -          -    -       (3)  The graph of   ,s Xq -  is plotted in Fig. 5. It is also  a periodic function of time t.  Similarly,  the  expressions  for  thermal  loadings  acting  on  the  surface  Z+  and  Z-  can  be  easily  obtained.  Here  we  present  two  other  loadings  for  the front and rear surfaces of solar array.  2.2.3  Front surface  For  the  front  surface,  only  solar  thermal  flux  ,s FSq  is present. It is determined as follows  0 0 0 , 0 1 0 if 0, 1 2 1 1 cos if 1 , 2 2 2 2 1 0 if 2 , 2 orb ec s FS F s orb orb orb orb t P q G t P P t P P                  -               = - -  - -                    -              (4)  The graph of   ,s FSq  is illustrated in Fig. 6.  Figure 6. Solar flux on the front surface  2.2.4  Rear surface  There  are  three  thermal  fluxes  that  affect  to  the  rear surface: solar, albedo and  infrared  fluxes. The  mathematical  representation  for  them  can  be  estimated as follows  70          SCIENCE & TECHNOLOGY DEVELOPMENT, Vol 20, No.K2- 2017   0 s,RS 0 0 0 0 if 0, 1 cos if , 1 2 2 1 1 0 if 1 , 2 2 2 1 cos 2 if 2 , 2 2 ec ec R s ec orb orb orb ec R s orb orb t P G t P P q t P P G t P P                        -  -          =        - -                 -   -            (5) 0 0 0 , , 0 1 0 if 0, 1 2 1 1 cos if 1 , 2 2 2 2 1 0 if 2 , 2 orb ec a RS R s RS e orb orb orb orb t P q G aF t P P t P P                  -               = -   - -                    -        (6) 4 ,RS RS,IR R e eq F T =                                             (7)  Here,  a   is  the albedo  factor of  the Earth;  eT   is  the Earth's black-body equivalent temperature. The  graphs  of  ,RSsq ,  ,RSaq   and  ,RSIRq   are  delineated  in  Fig. 7. In calculations of infrared loadings, we pay  attention  to  view  factor  between  the  satellite  and  the  Earth.  There  are  two  surfaces  (surface  Z-  and  rear  surface)  which  have  view  factors  differ  from  zero,  i.e.  ,e 0ZF -    and  , 0RS eF  .  Other  surface’s  view  factor  are  considered  equal  zero.  The  view  factor depends on the altitude of the satellite orbit.  Because the surface areas of the satellite are rather  small,  they  can  be  considered  as    differential  surfaces. So that, we  can compute view factor from  a  differential  surface  to  the  Earth’s  sphere.  On  calculating  view  factor,  readers  can  be  seen  in  detail in the book by Howell et al. [12].  Figure 7. Thermal fluxes on the rear surface  TABLE 3  VALUES OF  iC AND   ,dis iQ  FOR THERMAL CALCULATION  Node iC , J/K  ,dis iQ , W  1   702.1  10  2   702.1  10  3   702.1  10  4   702.1  10  5   702.1  10  6   702.1  10  7 1131.8  15  8 1131.8  15  2.2 Thermal balance equations  As  has  been  stated  in  the  previous  section,  our  satellite  can  be  thermally  modeled  with  eight   nodes.  Let  iC     be  the  thermal  capacities  of  the  nodes,  and  iT   be  their  temperatures  ( 1,...,8i = ).  The geometric model corresponding to this thermal  mathematical  model  is  shown  in  Fig.2.  The  nodes  are  thermally  coupled  by  both  conduction  and  radiation,  and  also  radiation  interaction  to  space  environment.  Let  ijK   be  the  conductive  coupling  coefficient  and  ijR     the  radiative  coupling  coefficient.  The  energy  balance  equations  for  the  nodes are       4 4 4 4 , , 1 1 n n i i ij j i ij j i i i ext i dis i j j C T K T T R T T R T T Q Q  = = = -  -  -    (8)  where  ,ext iQ     represents  the  external  thermal  load  on  the  node  i,  and  ,dis iQ   represents  the  heat  dissipation of the node i.   Taking  into  account  the  input  parameter  information in Tab. 2, and nodes are assumed to be  undergone a constant heat dissipation level (in W),  it  is  possible  to  calculate  the  capacities  iC ,  i=1,,8,  with  results  given  in  Tab.  3.  We  obtain  the following conduction (in W/K) and radiation (in  W/K4) matrices:  8 8ij K   = = K     0 0.1078 0 0.1078 0.1078 0.1078 0 0 0.1078 0 0.1078 0 0.1078 0.1078 0 0 0 0.1078 0 0.1078 0.1078 0.1078 0 0 0.1078 0 0.1078 0 0.1078 0.1078 0 0 0.1078 0.1078 0.1078 0.1078 0 0 0 0 0.1078 0.1078 0.1078 0.1078 0 0 0 0 0 0 0 0 0 0 0 32.55 0 0 0 0 0 0 32.55 0                          (9)  TẠP CHÍ PHÁT TRIỂN KH&CN, TẬP 20, SỐ K2-2017 71  8 8 8 10ijR -   = =  R    0 0.1971 0.1968 0.1971 0.1971 0.1971 0.1531 0 0.1971 0 0.1971 0.1968 0.1971 0.1971 0 0 0.1968 0.1971 0 0.1971 0.1971 0.1971 0 0 0.1971 0.1968 0.1971 0 0.1971 0.1971 0 0 0.1971 0.1971 0.1971 0.1971 0 0.1968 0 0 0.1971 0.1971 0.1971 0.1971 0.1968 0 0 0 0.1531 0 0 0 0 0 0.1531 0 0 0 0 0 0 0 0 0                    -       (10)  9 8 1 10iR -     = =  R     0.1162 0.1162 0.1162 0.1162 0.1162 0.1162 0.1730 0.1627 T (11)  The system (8) can be rewritten in the following  matrix form for eight-node model     4 4 4 ext dCT = KT + RT - R T + R T + Q + Q     (12)  where   1 2 8... T T T T=T is  a  generalized  vector; C is the extended thermal capacity and  K ,  R ,  R  are  the extended conduction and radiation  matrices  obtained  from  the  original  matrices  K ,  R ,  R  by  rearranging elements  ,i jT T  of Eq.  (8)  to  get  the  form  (12).     textQ   is  a  vector  of  external  thermal  loadings;  dQ   is  a  vector  of  the  dissipation  of  nodes;  4 4 4 4 1 2 8... T T T T =  T denotes  a  vector  of  radiation terms.   Assume  that  the  thermal  capacity  matrix  is  not  singular  (i.e.   det 0C ).  Pre-multiplying  both  side of Eq. (12)  by  1-C , we obtain   1 t-    =   4 4 4 ext dT C KT + RT - R T + R T + Q + Q              (13)  In  the  next  section,  we  will  solve  Eq.  (13)  using  Runge-Kutta  algorithm  to  get  numerical  solutions  of thermal responses of nodes.  3 RUNGE-KUTTA METHOD  The  Runge-Kutta  (RK)  method  is  one  of  the  most well-known methods for finding approximate  solutions  of  ordinary  differential  equations  in  problems  of  numerical  analysis.  It  was  developed  around  1900  by  two  German  mathematicians  C.  Runge and M. W. Kutta [13,14]. In this section, we  use  fourth-order  RK  method  to  find  temperature  responses  of  nodes  of  the  small  satellite  under  consideration. Our attention  is  to solve Eq.  (13)  to  explore  thermal  characteristics  of  nodes  when  thermal  nodes  are  subjected  to  radiation  loadings  from space environment.   We  consider  a  generalized  ordinary  differential  equations system in the following form   ,t=T F T                                                       (14)  where   1 2 ... T nT T T=T   is  a  generalized  vector,   1 2 ... T nF F F=F   is a nonlinear n-vector  function  which contains linear and nonlinear terms   1 2, , ,...,i i nF F t T T T= ,  1,i n= ,                      (15)  and  1 2, , ... , nT T T  are functions of time   0 , tFt t .  To  calculate  numerical  solutions,  we  divide  the  time  interval   0 , tFt   into  n  equal  segments  by  (n+1)  points  it : 0it t ih=  ;  n Ft t= ;  0Ft th n - = .  The  set  of  points  it     creates  a  "differential  net",  each  point  is  called  a  grid  node,  h   is  called  the  mesh step. We can estimate the approximate value  of   1 1i iT T t =   from   i iT T t=  as follows   1 1 2 3 4k 2 k 2 k k , 6 i i h T T =                          (16)  where      1 2 1 3 2 4 3 k , ; k , k ; 2 2 k , k ; 2 2 k , k . i i i i i i i i t T h h t T h h t T t h T h =   =        =      =   F F F F                                  (17)  The Runge-Kutta method will be applied to system  (13) where the function  F  is given by     1,t t-    =   4 4 4 ext dF T C KT + RT - R T + R T + Q + Q  (18)  Numerical results for Eq. (13) will be presented in  next section.  72          SCIENCE & TECHNOLOGY DEVELOPMENT, Vol 20, No.K2- 2017 4 NUMERICAL RESULTS AND DISCUSSION        Temperature  evolutions  in  time  of  eight  nodes  of satellite are shown in Fig. 8. It  is observed that,  when  the  satellite  is  in  the  illuminated  region  of  orbit,  temperatures  of  solar  array  are  larger  than  that  of  other  nodes.  The  maximum  temperature  is  predicted  to  be  82.4722  0C  for  node  8  (the  front  surface)  and  80.2318  0C  for  node  7  (the  rear  surface) of solar array (see Tab. 4). There is a slight  difference  between  temperatures  of  nodes  7  and  8  because  the  conduction  thickness  between  these  two  nodes  is  very  small  in  comparison  with  the  length of solar array. In the eclipse region of orbit,  the minimum temperatures of nodes 1, 2, 3, 4 and 6  are  nearly  the  same,  about  -60  0C.  In  this  region,  the  received  thermal  of  the  nodes 1,  2, 3,  4  and 6  from  environment  is  very  low.  The  change  of  temperatures  between  nodes  is  due  to  the  thermal  interactions  via  the  conduction  and  radiation.  For  three nodes 5, 7 and 8,  the minimum temperatures  are  higher  than  others.  By  taking  the  estimated  mean of maximum and minimum temperatures, we  can see that the estimated mean of node 5 is highest  because  it  always  remains  a  thermal  flux  acted by  Earth's infrared radiation (see Tab. 4).   Figure  8.  Temperature  evolutions  in  time  of  eight  nodes  of  satellite  TABLE 4  MINIMUM AND MAXIMUM TEMPERATURES OF NODES FOR THE  PARAMETERS’ VALUES GIVEN IN TABLE 2  Node Min ( C ) Max ( C ) Estimated mean ( C ) 1 -57.1998   2.2558  -27.4720  2 -59.0681  64.4913     2.7116  3 -59.0294  -3.5636  -31.2965  4 -58.3320  45.2248    -6.5536  5 -22.4952  67.4134    22.4591  6 -61.5191  62.9609     0.7209  7 -22.6112  80.2318   11.5160  8 -24.1244  82.4722   12.6362  It  is observed that,  the predicted maximum and  minimum temperatures of satellite's body and solar  array  belong  to  the  range  of  temperature  requirements given in Tab 5.  TABLE 5. THERMAL REQUIREMENTS [3]  Temperature Min (oC) Max (oC) Solar arrays  -100  +120  Inactive structure  -100  +100  Fig.  9  depicts  the  change  of  temperature  evolutions  of  node  8  with  various  values  of  absorbtivity  F .  The  value  of  F   shows  the  absorbed  part  of  incoming  radiation  to  the  total  incoming radiation of the front surface. In the case  of  large  value  of  F ,  maximum  temperature  of  node  8  grows  rapidly.  This  shows  that  a  strong  effect  of  the  absorbitivity  on  the  temperature  of  node  8  when  the  satellite  is  in  the  illuminated  region  of  orbit.  In  Fig.  10,  the  difference  between  the  maximum  temperatures  of  front  surface  with  various  values  of  F   when  compared  to  the  case  0.92F =   is  computed.  This  difference  value  is  largest  if  absorbtivity  F   takes  small  value,  for  example  0.1F = .  The  value  of  the  temperature  difference  is  reduced  if  the  solar  absorbtivity  F   increases.   Figure.  9.  Temperature  evolutions  of  front  surface  (node 8)  with various absorbtivity  αF  TẠP CHÍ PHÁT TRIỂN KH&CN, TẬP 20, SỐ K2-2017 73  Figure 10. The errors between the maximum temperatures of  front surface (node 8) with various values of αF when compared  to the case αF = 0.92 Fig.  11  presents  the  temperature  evolution  of  node  8  with  different  values  of  emissitivity  coefficient  F .  The  emissitivity  is  taken  from  0.5  to  0.82.  Three  values  are  selected,  0.5F = ,  0.7F =   and  0.82F = .  It  is  seen  that  as  the  emissitivity is increasing, the temperature of node 8  decreases.   The thermal interaction between two nodes 7 and 8  is  shown  in  Fig.  12.  The  dependence  of  node  8's  temperature  on  the  node  7's  temperature  is  nearly  linear  because  the  difference  of  temperatures  between them is quite small.    Figure 11. Temperature evolutions of  front  surface  (node 8)  with various emissitivity εF   Figure  12.  Temperature  of  front  surface  T8    versus  temperature of  rear surface T7 Figure  13.  Temperature  of  rear  surface  (node  7)  T7  versus  temperature  of    surface  Y+    (node  1)  T1,    with  various  absorbtivity αR In  Fig.  13,  we  portray  the  thermal  interaction  between nodes 1 (of satellite's body) and 7 (of solar  array)  with various values of node 7's absorbtivity  R . For each value of  R ,  in  the steady-state,  the  temperatures  of  node  1  and  7  can  reach  a  limit  cycle.  The  shape  of  the  limit  cycle  in  this  case  is  not a circle or ellipse because the obtained thermal  responses  of  node  1  and  7  are  not  harmonic  but  almost periodic.   5 CONCLUSION        The  results  of  analysis  for  thermal  characteristics  of  a  satellite  structure  have  turned  out  to  be  extremely  useful  in  the  framework  of  a  satellite  mission.  In  this  research,  a  simplified  model of satellite are carried out using  the method  of  lumped  parameters  for  nodes.  The  numerical  results of nodal temperatures are implemented with  the  use  of  the  Runge-Kutta  method.  Several  main  74          SCIENCE & TECHNOLOGY DEVELOPMENT, Vol 20, No.K2- 2017 results  are  obtained  and  can  be  summated  as  follows:  -  A  modeling  of  thermal  loadings  from  space  environment  is  established  in  the  framework  of  Low Earth Orbit.  - A simplified model with eight nodes representing  for  the  body  plates  and  solar  array  is  constructed  based  on  the  geometrical  dimensions  and  material  properties of satellite.  -  The  thermal  balance  equations  for  nodes  are  derived  from  the  characteristics of  conduction  and  radiation  interactions  between  nodes  and  external  thermal loadings.  -  The  temperature  evolutions  in  time  of  nodes  are  obtained  using  the  Runge-Kutta  algorithm  with  representations  of  the  extended  thermal  capacity,  conduction  and  radiation  matrices  obtained  from  the  rearrangement  of  thermal  nodes  in  thermal  balance equations.  -  The  effects  of  material  properties  such  as  absorbtivity  and  emissivity  on  the  thermal  responses of nodes are explored.  -  The  maximum  and  minimum  temperature  information  of  nodes  shows  that  the  predicted  temperatures of the satellite obtained from numeral  analyses are within the allowable temperature limit  range of satellite.  REFERENCES  [1] D.  G.  Gilmor,  "Spacecraft  Thermal  Control  Handbook",  The Aerospace Corporation, 2002.  [2] P.  Fortescue,  G.  Swinerd,  J.  Stark,  "Spacecraft  System  Engineering", John Wiley & Son Ltd., 2003.  [3] J.Meseguer,  I.  Pérez-Grande  and  A.  Sanz-Andrés,  “Spacecraft thermal control”, Woodhead Publishing, 2012  [4] K.  Oshima,  Y.  Oshima,  "Analytical  approach  to  the  thermal  design  of  spacecraft",  Institute  of  Space  and  Aeronautical Science of Tokyo, Report No. 419, 1968.  [5] C. Arduini, G.  Laneve,  S.  Folco,  "Linearized  techniques  for  solving  the  inverse  problem  in  the  satellite  thermal  control", Acta Astronautica, vol. 43, pp. 473-4789, 1998.  [6] J. Gaite, A. S. Andres,  I. P. Grande,  "Nonlinear analysis  of a simple model of temperature evolution in a satellite",  Nonlinear Dynamics, vol. 58, pp. 405-415, 2009.  [7] I.  P.  Grande,  A.  S.  Andres,  C.  Guerra,  G.  Alnonso,  "Analytical study of the thermal behaviour and stability of  a  small  satellite", Applied Thermal Engineering, vol. 29,  pp. 2567-2573, 2009.  [8] J. Gaite, A. S. Andres,  I. P. Grande,  "Nonlinear analysis  of a simple model of temperature evolution in a satellite",  Nonlinear Dynamics, vol. 58, pp. 405-415, 2009.  [9] J.  Gaite,  G.  Fernández-Rico,  "Linear  approach  to  the  orbiting  spacecraft  thermal  problem",  Journal  of  Thermophysics  and  Heat Transfer,  vol.  26, pp.  511-522,  2012.  [10] N. D. Anh, N. N. Hieu, P. N. Chung, N. T. Anh, "Thermal  radiation  analysis  for  small  satellites  with  single-node  model  using  techniques  of  equivalent  linearization",  Applied Thermal Engineering, vol. 94, pp. 607-614, 2016.  [11] J.  Gaite,  "Nonlinear  analysis  of  spacecraft  thermal  models",  Nonlinear  Dynamics,  vol.  65,  pp.  283-300,  2011.  [12] J. R.Howell, R. Siegel, M P. Menguc, "Thermal Radiation  Heat Transfer", Taylor & Francis Group, 2010.  [13] J.  C.  Butcher,  "Numerical  Methods  for  Ordinary  Differential Equations", New York:  John Wiley & Sons,  2008, ISBN 978-0-470-72335-7.  [14] A.  Iserles,  "A First Course  in  the Numerical Analysis of  Differential  Equations",  Cambridge  University  Press,  1996, ISBN 978-0-521-55655-2.  [15] P. N. Chung, N. N. Hieu, N. D. Anh, "Thermal radiation  analysis for solar arrays of a small satellite in Low Earth  Orbit", The 4th  International Conference on Engineering  Mechanics  and  Automation  (ICEMA  4)  Hanoi,  August  25÷26, 2016, pp. 146-153.  Pham Ngoc Chung was born in  Ninh  Binh  province  in  Viet  Nam. He graduated from Faculty  of  Mathematics,  Mechanics  and  Informatics  in  VNU  University  of  Science.  He  got  the  B.S.  and  M.S  degree  in  Mechanics  in  2010  and  2014,  respectively.  His  current  job  is  a  lecturer  in University of Mining and Geology. His  research  interest  consists  algebraic  systems,  numerical  simulation,  nonlinear  dynamical  systems, thermal analysis and control for satellites.  He  has  published  about  ten  scientific  papers  in  National Conferences and International Journals.   Nguyen Dong Anh  was  born  in  Hanoi,  Vietnam,  in  1954.  He  received  the  D.Sc  degree  in  Vibration  in  1986.  He  was  promoted  to  Professor  in  1996.  Currently,  he  is  working  at  the  Institute  of  Mechanics,  Vietnam  Academy  of  Science and Technology as Chairman of the Board  of  Scientists.  His  research  interest  includes  vibration,  nonlinear  random  oscillation,  oscillation  control,  nonlinear  dynamical  systems.  He  has  published  more  than  100  scientific  articles  in  National  and  International  Journals.  He  was  also  the  author  of  two  monograph  books.  He  successfully taught and educated many masters and  doctors.  TẠP CHÍ PHÁT TRIỂN KH&CN, TẬP 20, SỐ K2-2017 75  Nguyen Nhu Hieu  was  born  in  Bac  Ninh  province,  Vietnam.  He  received  the  B.S.  and  M.S.  degrees  in  Mechanics  of  Solids  from  the  Hanoi  National  University  in  2008  and  2011,  respectively. At present, he works at the Institute of  Mechanics,  Vietnam  Academy  of  Science  and  Technology.  His  current  areas  of  interest  include  applied  mathematics  and  nonlinear  dynamical  systems.  He  has  published  more  than  twenty  scientific  papers  in  National  Conferences  and  International Journals.  76          SCIENCE & TECHNOLOGY DEVELOPMENT, Vol 20, No.K2- 2017 Tóm tắt - Trong bài báo này, đáp ứng nhiệt phi tuyến của một vệ tinh nhỏ trên quỹ đạo thấp của Trái đất được phân tích dựa trên mô hình nhiệt nhiều nút. Các thành phần kết cấu chính của một vệ tinh dạng hình hộp gồm có thân với sáu mặt hình chữ nhật và một cánh nối với thân. Các thành phần thân và cánh có thể được mô hình hóa trên cơ sở phương pháp nhiệt phân bổ, nghĩa là mỗi mặt của thân và cánh được đặc trưng bởi một nút nhiệt. Để ước lượng nhiệt độ cho các thành phần này, chúng ta có thể sử dụng mô hình nhiệt tám nút: sáu nút cho các mặt của thân và hai nút cho mặt trước và mặt sau của cánh. Các nút hấp thụ ba nguồn nhiệt chủ yếu từ môi trường không gian bao gồm bức xạ mặt trời, bức xạ albedo và bức xạ hồng ngoại Trái đất. Hệ phương trình cân bằng nhiệt xác lập cho các nút là hệ phương trình vi phân phi tuyến và được giải bằng một phương pháp số. Với mục đích mô phỏng, giả sử rằng vệ tinh luôn duy trì ở tư thế “Earth-pointing” trong suốt thời gian nó chuyển động trên quỹ đạo.Tiến triển nhiệt độ theo thời gian của các nút sẽ được nghiên cứu một cách chi tiết. Kết quả thu được chỉ ra giá trị nhiệt độ dự đoán của các nút nằm trong giới hạn nhiệt cho phép của vệ tinh. Từ khóa - vệ tinh nhỏ, Earth-pointing, đáp ứng nhiệt, giới hạn nhiệt.  Phân tích ứng xử nhiệt phi tuyến của vệ tinh nhỏ  trên quỹ đạo thấp sử dụng mô hình nhiều nút  Phạm Ngọc Chung, Nguyễn Đông Anh, Nguyễn Như Hiếu