Simulation method for calculating rigid-flexible-electric-thermal coupling dynamic response of spatial structure

Abstract

The invention discloses a simulation method for calculating rigid-flexible-electric-thermal coupling dynamic response of a space structure, which comprises the following steps of: simplifying the space structure into a central rigid body-cantilever slab model, describing and dispersing the model through a constructed absolute node method mixing unit, deducing a dynamic equation of the model and a temperature field and electrostatic field control equation, and calculating the dynamic response of the rigid-flexible-electric-thermal coupling dynamic response of the space structure. The method comprises the following steps of: obtaining a rigid-flexible-electric-thermal coupling equation set of the system by taking a dynamic response result of the structure into consideration, considering a bidirectional coupling effect between a displacement field and a physical field, solving the equation set by utilizing a solver based on a generalized alpha method, obtaining the dynamic response result of the structure, performing visualization processing, and displaying the change condition of physical quantities such as temperature, potential and displacement along with time. According to the method, the influence of physical field factors is comprehensively considered, the spatial structure dynamic behavior under multi-physical field coupling can be accurately simulated, and a certain basis is provided for the simulation design of the spatial structure.

Description

Technical Field

[0001] The present invention belongs to the field of multi-body system dynamics modeling, and specifically relates to a simulation method for calculating the rigid-flexible-electro-thermal coupled dynamic response of a spatial structure. Background Art

[0002] With the rapid development of space technology, spatial structures are increasingly widely used in spacecraft such as remote sensing satellites and optical satellites. These structures usually have the characteristics of large size, light weight, high flexibility and low damping, and are extremely likely to exhibit complex non-linear dynamic behaviors in the space environment, seriously affecting the attitude stability and mission accuracy of spacecraft. For example, the Chinese invention patent with the publication number CN116384198A discloses a "simulation method for the dynamic behavior of a rigid-flexible-thermal coupled thin plate system", which can capture complex dynamic responses such as thermally induced vibrations of the thin plate structure under the action of solar thermal radiation. However, this method only simplifies the spatial structure into a thin plate model and does not consider the influence of the spacecraft rigid body on the flexible appendage structure. This simplification limits its application scope in the modeling of complex rigid-flexible coupling systems and cannot comprehensively describe the dynamic behavior of spatial structures.

[0003] In addition, as a typical electronic equipment, the dynamic behavior of a spatial structure is not only significantly affected by the thermal field, but also by the electric field. Especially in spacecraft with high requirements for attitude stability, the influence of the electric field effect on the dynamic characteristics of the structure cannot be ignored. At the same time, considering the coupling effect between the electric field and the structure can also achieve precise control of the flexible spatial structure. However, existing research still lacks the unified modeling of the temperature field, displacement field and electric field within the same framework, and cannot take into account the coupling characteristics and non-linear dynamic behaviors of rigid-flexible-electro-thermal multi-physical fields at the same time. This limitation results in low calculation accuracy of existing analysis methods under multi-physical field coupling conditions and cannot meet the demand for accurate dynamic response prediction of spatial structures. Summary of the Invention

[0004] The purpose of the present invention is to propose a simulation method for calculating the rigid-flexible-electro-thermal coupled dynamic response of a spatial structure.

[0005] The technical solution for achieving the purpose of the present invention is as follows: A simulation method for calculating the rigid-flexible-electro-thermal coupled dynamic response of a spatial structure, comprising the following steps:

[0006] Step 1, simplify the spatial structure into a central rigid body-cantilever plate physical model, where the central rigid body-cantilever plate physical model includes a central rigid body simplified from a satellite and a flexible thin plate simplified from a flexible appendage; set the simulation parameters of the central rigid body-cantilever plate physical model;

[0007] Step 2: Describe and discretize the central rigid body - cantilever plate model based on the mixed element; obtain the element mass matrix, element generalized elastic force matrix of the physical model of the central rigid body - cantilever plate according to the theory of continuum mechanics, and obtain the element dynamic equation of the model; meanwhile, describe and discretize the temperature field control equation and the electrostatic field control equation based on the mixed element.

[0008] Step 3: Realize the bidirectional coupling between the temperature field and the displacement field by introducing the radiation angle coefficient, and realize the bidirectional coupling between the electrostatic field and the displacement field through the piezoelectric constitutive equation; finally, link the dynamic equation of the physical model of the central rigid body - cantilever plate, the temperature field control equation and the electrostatic field control equation through the bidirectional coupling effect, and then systematically assemble the element equations based on the finite element method to finally obtain the rigid - flexible - electro - thermal coupling system control equations of the physical model of the central rigid body - cantilever plate.

[0009] Step 4: Construct a multi - physical - field hybrid solver based on the generalized α - method to solve the rigid - flexible - electro - thermal coupling system control equations of the physical model of the central rigid body - cantilever plate, and obtain the dynamic response results of the position of the central rigid body - cantilever plate model and the physical information of the temperature and electric potential of the model through iterative calculation.

[0010] Step 5: Visualize the obtained dynamic response results and the physical information of the temperature and electric potential of the model to obtain the changes of the end displacement, lateral deformation, temperature, electric potential of the test points on the central rigid body - cantilever plate model and the angle of the central rigid body over time.

[0011] Compared with the prior art, the present invention has the following remarkable advantages: (1) Based on the absolute nodal coordinate method, a mixed element integrating temperature and electric potential is constructed, and the mixed element is used to discretize the central rigid body - cantilever plate system, realizing the unified modeling of the temperature field, displacement field and electric field under the same mesh division framework, and obtaining the dynamic response of the structure with the coupling characteristics of rigid - flexible - electro - thermal multi - physical fields. (2) In the process of modeling the space structure, the influence of the satellite rigid body on the dynamic behavior of the structure is considered, and more accurate simulation results that are more in line with engineering practice can be obtained. (3) The transient dynamic response of the structure under the action of multi - physical - field factors can be obtained, and the influence of different physical - field factors on the space structure can be analyzed. Brief Description of the Drawings

[0012] The present invention will be further described below in conjunction with the drawings and embodiments.

[0013] Figure 1 is the flow chart of the present invention.

[0014] Figure 2 is the schematic diagram of the central rigid body - cantilever plate model.

[0015] Figure 3It is a flow chart of a solver based on the generalized α method.

[0016] Figure 4 It is a visual interface for program parameter input.

[0017] Figure 5 It is the displacement, temperature, and electric potential changes of the end point on the central rigid body - cantilever plate model under the action of heat flow, as well as the change of the angle of the central rigid body. Figure 6 It is a schematic diagram of the change in the rotation angle of the central rigid body under the action of heat flow. Figure 7 It is a schematic diagram of the change in the temperature gradient of the end point on the central rigid body - cantilever plate under the action of heat flow. Detailed implementation manners

[0018] In order to make the objectives, technical solutions, and advantages of the present application clearer, the following further details the present application in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present application and are not used to limit the present application.

[0019] As Figure 1 shown, a simulation method for calculating the rigid - flexible - electro - thermal coupling dynamic response of a spatial structure includes the following steps:

[0020] Step 1: Simplify the spatial structure into a central rigid body - cantilever plate physical model, where the central rigid body - cantilever plate physical model includes a central rigid body simplified from a satellite and a flexible thin plate simplified from a flexible attachment; set the simulation parameters of the central rigid body - cantilever plate physical model, and the specific content is as follows:

[0021] (1) Geometric parameters

[0022] The geometric parameters of the central rigid body - cantilever plate system include: the radius of the central rigid body is R, and the moment of inertia of the central rigid body is J R ; the geometric parameters of the cantilever plate: the length is L, the width is W, and the thickness is H; the geometric parameters of the piezoelectric patch: the length is l, the width is w, and the thickness is h.

[0023] (2) Material parameters:

[0024] The material parameters of the cantilever plate: the density is ρ, the elastic modulus is E, and the Poisson's ratio is υ;

[0025] The material parameters of the piezoelectric patch: the density is ρ P and the elastic modulus is E P and the Poisson's ratio is υ P ;

[0026] (3) Mesh parameters:

[0027] The mesh parameters of the cantilever plate are as follows: the total number of mesh elements is N = N1×N2, the length of the mesh element is l e = L / N1, and the width of the mesh element is w e = W / N2; the same mesh parameters are adopted for the temperature field and the electrostatic field.

[0028] (4) Physical field parameters:

[0029] Parameters related to the temperature field: specific heat capacity с, thermal conductivity λ, coefficient of thermal expansion α, solar heat radiation density q T , surface heat absorption coefficient μ;

[0030] Parameters related to the electrostatic field: piezoelectric constant ζ, dielectric constant ξ.

[0031] Step 2: Describe and discretize the central rigid body - cantilever plate model based on the hybrid element; obtain the element mass matrix, element generalized elastic force matrix of the physical model of the central rigid body - cantilever plate according to the theory of continuum mechanics, and obtain the element dynamic equation of the model; at the same time, describe and discretize the temperature field control equation and the electrostatic field control equation based on the hybrid element, and the specific process is as follows:

[0032] Step 2.1: Describe and discretize the central rigid body - cantilever plate model based on the hybrid element, and the specific process is as follows:

[0033] In the inertial coordinate system, any point on the central rigid body - cantilever plate model can be described as follows:

[0034]

[0035] In the formula, the subscript h indicates that it is a hybrid element; r1, r2, r3 are the components of the position vector of this point in the inertial coordinate system O-XYZ in each direction; are the temperature and electric potential coordinates respectively; S h and q h are the element shape function and element generalized coordinate of the constructed hybrid element respectively.

[0036] When discretizing the central rigid body - cantilever plate model, each hybrid element on it consists of four nodes, and the specific form of the position coordinate array of each node is as follows:

[0037]

[0038] In the formula, are the position coordinate vectors, temperature and electric potential respectively. and are the gradient vectors in each direction respectively;

[0039] Step 2.2: Obtain the expressions of the element mass matrix, element generalized elastic force matrix, element generalized external force matrix, etc. of the model through derivation, and obtain the element dynamic equation of the model. The specific process is as follows:

[0040] According to the kinetic energy theorem, the kinetic energy T of the element e can be expressed as:

[0041]

[0042] In the above formula, ρ is the material density, is the first derivative of the position vector, is the first derivative of the column vector of node coordinates, and M e1 is the element mass matrix.

[0043] Therefore, the mass matrix M of the thin plate element can be obtained e1 , and its expression is as follows:

[0044]

[0045] In the formula, S is the element shape function, h is the thickness of the thin plate, x and y are the element length and width respectively;

[0046] The elastic energy U of the element e is in the following specific form:

[0047]

[0048] In the formula, ε g is the Green-Lagrange strain of the thin plate element, κ is the curvature of the thin plate element, and E ε , E κ is the elastic coefficient matrix of the thin plate element;

[0049] The element generalized elastic force Q e vector can be obtained by taking the partial derivative of the element elastic energy U e with respect to the generalized node coordinate q. The generalized elastic force vector of the element can be expressed as:

[0050]

[0051] In addition, the flexible thin plate and the central rigid body are connected by a fixed support method, and the relationship between the rotation angle and the constraints of the left-end node of the thin plate is as follows:

[0052]

[0053] In the formula, R is the radius of the central rigid body, q i1 , q i3 , q i4 , q i6are the partial generalized coordinates of the i-th node at the left end of the thin plate, and θ h is the rotation angle of the central rigid body.

[0054] Furthermore, the relationship between the angular velocity ω of the central rigid body and the generalized nodal coordinates at the left end of the thin plate can be obtained as follows:

[0055]

[0056] For the central rigid body connecting plate system with large-range motion, when the central rigid body performs fixed-axis rotation, its kinetic energy can be expressed as:

[0057]

[0058] In the formula, ω is the angular velocity of the central rigid body, and J oh is the moment of inertia of the central rigid body;

[0059] The total mass matrix of the central rigid body cantilever plate system is:

[0060]

[0061] In addition, due to the rigid assumption of the central rigid body, the generalized elastic force of the system remains unchanged.

[0062] Only by assembling the generalized nodal arrays, mass matrices, elastic force matrices, and generalized external force matrices of each thin plate in the central rigid body - cantilever plate system can the dynamic equation of the entire central rigid body - cantilever plate system be obtained, and its specific expression form is as follows:

[0063]

[0064] In the formula, λ is the Lagrange multiplier, F is the element external force matrix, and Φ(q, t) is the system constraint, which includes the fixed support constraint between the central rigid body and the flexible thin plate. The system constraint can be processed by the Lagrange multiplier method. Step 2.3: Describe and discretize the temperature field control equation and the electrostatic field control equation based on the hybrid element, and the specific process is as follows:

[0065] Step 2.3.1: Discretize the element temperature field control equation based on the hybrid element, and the specific content is as follows: The control equation of the transient temperature field is as follows:

[0066]

[0067] In the formula, Q is the internal heat source density of the system, ρ is the material density, c is the specific heat capacity of the material, k is the thermal conductivity of the material, and T is the temperature at any point on the thin plate;

[0068] Assume that the temperature at any point on the thin plate is T, and it can be expressed as:

[0069] T = T09TN +zT M

[0070] Taking the variation of the temperature at any point on the thin plate, the variational equation of heat conduction in the Cartesian coordinate system can be obtained. Then, taking the functional extremum of it, the variational equation of heat conduction in the following form can be obtained:

[0071]

[0072] The variational control equation of the temperature field after discrete treatment by the hybrid element is as follows:

[0073]

[0074] For the temperature gradient T M , through a similar process, the control equation after its discrete treatment by the hybrid element can also be derived:

[0075]

[0076] In the formula, q1, q2, and q3 represent three thermodynamic boundary conditions of the given boundary temperature, the given boundary heat flux, and the boundary being a convective heat transfer boundary respectively, q N is the generalized coordinate of temperature, q M is the generalized coordinate of the temperature gradient, z is the thickness of the thin plate, and S is the element shape function.

[0077] After further arrangement, the matrix form of the temperature field control equation can be obtained as follows:

[0078]

[0079] Each item in the formula can be further expressed as:

[0080]

[0081] In the formula, q T is the generalized coordinate matrix of the temperature field, C is the heat capacity matrix, K c is the heat conduction matrix, Q T is the heat source matrix, and the specific form is as follows:

[0082]

[0083] Step 2.3.2: Discretize the electrostatic field control equation based on the hybrid element, and the specific content is as follows:

[0084] The control equation of the electrostatic field, and the specific form is as follows:

[0085]

[0086] In the formula, ρ is the charge density, ε is the permittivity, is the electric potential.

[0087] Through the functional of the Poisson equation and discretization using mixed elements, the following equation form is obtained:

[0088]

[0089] where ε r is the relative permittivity

[0090] Taking the extremum of the functional and arranging it, its matrix form can be obtained:

[0091]

[0092] where the specific representation forms of each term are as follows:

[0093]

[0094] where λ is the permittivity matrix and q a is the charge density.

[0095] Step 3: Achieve the bidirectional coupling between the temperature field and the displacement field by introducing the radiation view factor, and achieve the bidirectional coupling between the electrostatic field and the displacement field through the piezoelectric constitutive equation; finally, connect the dynamic equation, the temperature field control equation, and the electrostatic field control equation of the central rigid body - cantilever plate physical model through the bidirectional coupling effect, and then based on the finite element method, systematically assemble the element equations to finally obtain the rigid - flexible - electro - thermal coupling system control equations of the central rigid body - cantilever plate physical model;

[0096] Step 3.1: Achieve the bidirectional coupling between the temperature field and the displacement field by introducing the radiation view factor, and the specific process is as follows:

[0097] In practical applications, when a spacecraft operates in a complex space environment, the temperature field and the displacement field do not act independently but are coupled with each other. The thermal load generated by the temperature field will cause structural deformation, and the structural deformation will in turn affect the boundary conditions of the temperature field.

[0098] The elastic energy caused by thermal stress The expression is as follows:

[0099]

[0100] where E is the elastic coefficient matrix of the thin plate, ε T is the thermal strain, σ is the element stress, and ε is the element strain; taking the derivative of the elastic energy caused by thermal stress with respect to the generalized coordinates can obtain the generalized temperature load, which is specifically expressed as:

[0101]

[0102] The radiation input in the space environment mainly comes from solar thermal radiation, which is usually assumed to irradiate the system structure in the form of parallel light rays; in order to describe the bidirectional coupling effect between the deformation of the flexible body and environmental radiation, the radiation angle coefficient Γ is introduced. α , and the specific form is as follows:

[0103]

[0104] In the formula, s is the unit vector of solar radiation heat flux; n is the normal vector of the thin plate element surface;

[0105] The radiation angle coefficient Γ α relates the solar radiation to the generalized nodal coordinates of the displacement field, thus realizing the bidirectional coupling between the temperature field and the displacement field.

[0106] Step 3.2: Realize the bidirectional coupling between the electrostatic field and the displacement field through the piezoelectric constitutive equation. The specific process is as follows: Based on the piezoelectric constitutive equation, the specific coupling relationship between the displacement field and the electrostatic field can be further derived; the common piezoelectric constitutive equation can be expressed as:

[0107]

[0108] In the formula, σ is the unit stress vector of the element, ε g is the unit strain vector of the element, c E is the elastic matrix, E is the electric field strength, D is the electric displacement, e is the piezoelectric constant matrix, and ζ is the dielectric constant matrix;

[0109] The strain energy of the system can be expressed as:

[0110]

[0111] Substitute the piezoelectric constitutive equation into the expression of the strain energy, and take the derivative of the strain energy of the system with respect to the displacement vector q to obtain the generalized elastic force as:

[0112]

[0113] In the formula, ε g is the unit strain of the element, c E is the elastic matrix, is the electric field strength, and e is the piezoelectric constant matrix; the electric field energy of the electrostatic field is expressed as follows:

[0114]

[0115] Take the derivatives of the electric field energy of the system with respect to the electric potential respectively, then the electric charge quantity of the electrostatic field can be obtained

[0116]

[0117] In the formula, Q c is the electric charge of the electric field; when there is no input voltage from the outside world, the electric charge Q of the electric field c is zero. At this time, solving the above formula can obtain the piezoelectric potential generated by the structural deformation.

[0118] Step 3.3: Connect the dynamic equation, temperature field control equation, and electrostatic field control equation of the system through the bidirectional coupling effect, and synthesize the rigid-flexible-electro-thermal coupling system control equations of the system through element assembly. The specific content is as follows:

[0119] Through element assembly, the rigid-flexible-electro-thermal coupling system control equations of the multi-body system based on the hybrid element are given. This set of equations consists of three parts: the dynamic equation, the temperature field control equation, and the electrostatic field control equation. The specific representation is as follows:

[0120]

[0121] In the formula, M qq is the system mass matrix; K qq , K qφ , K qT , F q are the system stiffness matrix, the force-electricity coupling matrix, the force-thermal coupling matrix, and the system external force matrix respectively; K φφ , F φ are the dielectric stiffness matrix and the charge load column vector respectively; C, K TT , F T are the system heat capacity matrix, the heat conduction matrix, and the heat flux load column vector respectively; Φ(q, t), Φ(q T , t), Φ(q φ , t) are the system constraints, the temperature field boundary conditions, and the electrostatic field boundary conditions respectively.

[0122] Step 4: Construct a multi-physics field hybrid solver based on the generalized α method to solve the rigid-flexible-electro-thermal coupling system control equations of the system, and obtain the dynamic response results such as the displacement, velocity, and acceleration of the central rigid body-cantilever plate model and the physical information of the model temperature and electric potential through iterative calculation. The specific calculation process is as follows:

[0123] Step 4.1: Construct a multi-physics field hybrid solver based on the generalized α method. The specific content is as follows:

[0124] The hybrid solver constructed in this paper based on the generalized-α method nests two series of inner loops in the iterative loop of a single integration time step, which are used to solve the physical field information and the system dynamic equation respectively;

[0125] Step 4.1.1: Perform iterative calculations on the system dynamics equations as follows:

[0126] The system dynamics equations are as follows:

[0127]

[0128] In the formula, M qq is the system mass matrix; K qq , K qφ , K qT , F q are the system stiffness matrix, the force - electric coupling matrix, the force - thermal coupling matrix, and the system external force matrix respectively. q, q T, q e are the displacement coordinate array, the temperature coordinate array, and the electric potential coordinate array respectively. λ is the Lagrange multiplier, and Φ(q, t) is the system constraint;

[0129] For simplicity of expression, denote: Q = K qq q + K qφ q e + K qT q T - F q

[0130] First, introduce the following parameters and give the iterative initial values at t = 0, including the generalized displacement, generalized velocity, generalized acceleration, and the initial Lagrange multiplier. The specific forms are as follows:

[0131]

[0132] After experiencing a time t = n×h, use the dynamic response of the previous moment to estimate the n+1 at the next time step t = t The discrete format is as follows:

[0133]

[0134] In the formula, are the derivatives of the unit node displacement coordinates and the mixed unit node displacement coordinates with respect to time at the previous time step respectively; are the unit displacement node coordinates and the derivatives of the unit node displacement coordinates with respect to time at the current time step. a n , a n+1 are the accelerations at the previous time step and the current time step respectively. h is the integration time step of the generalized - α method; α m , α f , β, γ are the numerical solution parameters of the generalized - α method;

[0135] Substitute Substituting into the kinetic equation gives:

[0136]

[0137] Solve using Newton - Raphson iteration:

[0138]

[0139] Wherein, is the Jacobi matrix of q at the k - th iteration step, n+1 and is the iteration increment; in the generalized - α method, a correction factor also needs to be introduced. The iteration result of the (k + 1)-th after correction is as follows:

[0140]

[0141] Substitute the result after (k + 1) iterations and calculate If it means that the accuracy of the iteration result meets the set requirements. Take as the initial condition for the next time step and update the parameter a n+1 and start the loop iteration for the next time step until the numerical simulation time ends;

[0142]

[0143] Step 4.1.2: Perform iterative calculation on the temperature - field control equation, and the specific content is as follows:

[0144] The specific forms of the temperature - field control equation and its constraint conditions are as follows:

[0145]

[0146] Wherein, C, K TT , F T are the system heat - capacity matrix, heat - conduction matrix, and heat - flux load column vector respectively, and Φ(q T , t) is the temperature - field boundary condition.

[0147] Give the initial iteration values of the temperature q T of the temperature field, the temperature change rate and the Lagrange multiplier λ T as follows:

[0148]

[0149] After experiencing time t = n×h, use the dynamic response at the previous moment to estimate the n+1 at the next time step t = t The discrete format is as follows:

[0150]

[0151] Take the estimated iterative value as the next iterative value as the initial value, and substitute it into the temperature field control equation to obtain

[0152]

[0153] Furthermore, it can be obtained that:

[0154]

[0155] In the formula, is the Jacobi matrix of the k-th for , iteration increment;

[0156] Introduce a correction factor The iterative result of the (k + 1)-th time after correction is as follows:

[0157]

[0158] Substitute the result after (k + 1)-th iteration and calculate If it means that the accuracy of the iterative result meets the set requirements. Take as the initial condition for the next time step, and start the loop iteration of the next time step until the end of the numerical simulation time;

[0159] In view of the fact that the iterative process of the electrostatic field control equation is the same as that of the temperature field control equation, no further explanation will be given here.

[0160] Step 4.2: According to the iterative steps in Step 4.1, write a solution program, and obtain the dynamic response results such as displacement, velocity, and acceleration of the central rigid body - cantilever plate model and the physical information of the model temperature and electric potential through iterative calculation.

[0161] Step 5, perform visualization processing on the obtained dynamic response results such as displacement, velocity, and acceleration and the physical information of the model temperature and electric potential to obtain the variation of the end displacement, lateral deformation, temperature, electric potential of the test points on the model, and the angle of the central rigid body with time. The specific operations are as follows:

[0162] Step 5.1: Visualize the iterative results of the displacements, temperatures, electric potentials, etc. of the grid nodes obtained in Step 4, and use software such as MATLAB and Origin to plot the variations of the end displacements, lateral deformations, temperatures, electric potentials, and central rigid body angles of the test points over time.

[0163] Embodiment

[0164] To verify the effectiveness and accuracy of the present invention, this example simulates the dynamic response of a central rigid body cantilever plate with a piezoelectric structure under solar heat flux.

[0165] In this example, the radius of the central rigid body is 1 m, and its moment of inertia is 50 kg·m 2 . In addition, the mesh division method in this example is to divide 20 elements along the x direction and 1 element along the y direction.

[0166] The geometric parameters, material parameters, and physical field parameters of the system required in Step 1 are shown in Table 1.

[0167] Table 1 Geometric parameters, material parameters, and physical field parameters of the central rigid body - cantilever plate system

[0168]

[0169] The present invention simplifies the spatial structure into a central rigid body - cantilever plate model, and then realizes the unified description and modeling of the rigid - flexible - electro - thermal coupling of the spatial structure through the constructed absolute nodal coordinate method hybrid element integrating temperature and electric potential, obtaining the system's rigid - flexible - electro - thermal coupling equations. For this system of equations, a solver is designed and solved based on the generalized α method, and the corresponding program is written using MATLAB, and the corresponding program interface is made. As shown in the figure, after inputting the corresponding parameter information and clicking the run button, the solution can be carried out. After the solution is completed and entering the post - processing interface, the dynamic response of the model can be visually displayed.

[0170] The technical features of the above embodiments can be combined arbitrarily. For the sake of brevity of description, not all possible combinations of the technical features in the above embodiments are described. However, as long as there is no contradiction in the combination of these technical features, it should be considered as the scope recorded in this specification.

[0171] The above - described embodiments only represent several implementation manners of the present application. The description is relatively specific and detailed, but it cannot be understood as a limitation to the scope of the present application. It should be noted that for those of ordinary skill in the art, without departing from the concept of the present application, several deformations and improvements can still be made, and these all belong to the protection scope of the present application. Therefore, the protection scope of the present application should be subject to the appended claims.

Claims

1. A simulation method for calculating the rigid-flexible-electrical-thermal coupled dynamic response of a spatial structure, characterized in that: The following steps are involved: Step 1, simplifying the space structure into a central rigid body-cantilever plate physical model, wherein the central rigid body-cantilever plate physical model includes a central rigid body simplified from a satellite and a flexible thin plate simplified from a flexible attachment; setting simulation parameters of the central rigid body-cantilever plate physical model; Step 2, based on the hybrid unit, the central rigid body-cantilever plate model is described and the unit discretization is performed; according to the theory of continuum mechanics, the unit mass matrix and the unit generalized elastic force matrix of the central rigid body-cantilever plate physical model are obtained to obtain the unit dynamic equation of the model; at the same time, the temperature field control equation and the electrostatic field control equation are described and discretized based on the hybrid unit; Step 3, bidirectional coupling between temperature field and displacement field is realized by introducing radiation angle coefficient, and bidirectional coupling between electrostatic field and displacement field is realized by piezoelectric constitutive equation; finally, the bidirectional coupling effect is used to link the dynamic equation, temperature field control equation and electrostatic field control equation of the central rigid body-cantilever plate physical model, and then the unit equations are systematically assembled based on the finite element method, and finally the rigid-flexible-electric-thermal coupling system control equation group of the central rigid body-cantilever plate physical model is obtained; Step 4: Based on the generalized α method, a multi-physics hybrid solver is constructed to solve the rigid-flexible-electric-thermal coupling system control equations of the central rigid body-cantilever plate physical model, and the position dynamic response results and model temperature and potential physical information of the central rigid body-cantilever plate model are obtained through iterative calculation; Step 5: Visualize the obtained dynamic response results and the model temperature and potential physical information to obtain the change of the end displacement, lateral deformation, temperature, potential and central rigid body angle of the test point on the central rigid body-cantilever plate model over time.

2. The simulation method for calculating the rigid-flexible-electrical-thermal coupled dynamic response of a spatial structure according to claim 1, characterized in that: The simulation parameters of the central rigid body-cantilever plate physical model set in step 1 include geometric parameters, material parameters, meshing parameters and physical field parameters, which are specifically: The geometric parameters of the central rigid body-cantilever plate system include: the radius of the central rigid body is R, the moment of inertia of the central rigid body is J oh ;Cantilever plate geometric parameters: length L, width W, thickness H; piezoelectric sheet geometric parameters: length l, width w, thickness h; Material parameters of the cantilever plate: density ρ, elastic modulus E, Poisson's ratio v; Piezoelectric material parameters: density is ρ P , the elastic modulus is E P , Poisson's ratio is v P ; The grid parameters of the cantilever plate are: the total number of grid cells is N = N1 × N2, the length of the grid cell is le = L / N1, and the width of the grid cell is we = W / N2; the temperature field and the electrostatic field use the same grid parameters; Temperature field related parameters: specific heat capacity с, thermal conductivity λ, thermal expansion coefficient α, solar thermal radiation density θ T , surface heat absorption coefficient μ; Electrostatic field related parameters: piezoelectric constant ζ, dielectric constant ξ.

3. The simulation method for calculating the rigid-flexible-electric-thermal coupled dynamic response of a spatial structure according to claim 1, characterized in that: Step 2 The specific process of describing and discretizing the central rigid body-cantilever plate model based on the hybrid unit is as follows: In the inertial coordinate system, any point on the central rigid body-cantilever plate model is described as follows: Where, the subscript h indicates a hybrid unit; r1, r2, r3 are the components of the position vector of the point in the inertial coordinate system O-XYZ in each direction; T, are the temperature and potential coordinates respectively; S h and q h are the unit shape function and unit generalized coordinates of the constructed hybrid unit respectively; when the central rigid body-cantilever plate model is discretized, each hybrid unit on the central rigid body-cantilever plate model consists of four nodes, and the specific form of the position coordinate array of each node is as follows: In the formula, r i ,T i , are the node position coordinate vector, temperature and electric potential respectively, and are the gradient vectors in each direction respectively.

4. The simulation method for calculating the rigid-flexible-electric-thermal coupled dynamic response of a spatial structure according to claim 3, characterized in that: According to the theory of continuum mechanics, the unit mass matrix, unit generalized elastic force matrix and unit generalized external force matrix expressions of the central rigid body-cantilever plate physical model are obtained, and the specific method for obtaining the unit dynamics equation of the model is as follows: Determine the kinetic energy T of the mixing unit e , specifically: In the above formula, ρ is the material density, is the first-order derivative of the position vector, is the first-order derivative of the node coordinate column vector, M e1 is the unit mass matrix; The unit mass matrix M of the thin plate is obtained according to the kinetic energy of the hybrid unit e1 , specifically: Where S is the element shape function, h is the thickness of the thin plate, and x and y are the element length and width respectively; Determine the unit elastic energy U e , the specific form is as follows: In the formula, ε g is the Green Lagrangian strain of the thin plate element, κ is the curvature of the thin plate element, E ε ,E κ is the elastic coefficient matrix of the thin plate element; The unit elastic energy U e Taking partial derivatives of the generalized node coordinates q, we can obtain the unit generalized elastic force vector Q e , the generalized elastic force vector of the unit is specifically: For the moving central rigid body-cantilever plate physical model, when the central rigid body rotates along a fixed axis, the kinetic energy is expressed as: Where R is the radius of the central rigid body, ω is the angular velocity of the central rigid body, and J oh is the moment of inertia of the central rigid body; the generalized mass matrix of the central rigid body-cantilever plate system unit is: The unit generalized mass matrix and unit generalized elastic force matrix in the central rigid body-cantilever plate physical model are assembled to obtain the unit dynamics equation of the entire central rigid body-cantilever plate system. The specific expression is as follows: In the formula, Φ(q, t) is the system constraint, λ is the Lagrange multiplier, and F is the unit external force matrix. The temperature field control equation and the electrostatic field control equation are described and discretized to obtain the unit temperature field control equation and the unit electrostatic field control equation based on the hybrid unit discretization.

5. The simulation method for calculating the rigid-flexible-electrical-thermal coupled dynamic response of a spatial structure according to claim 4, characterized in that: The specific process of describing and discretizing the temperature field control equation and the electrostatic field control equation based on the hybrid unit is as follows: The unit temperature field control equation is discretized based on the hybrid unit, specifically: The governing equation for the transient temperature field is as follows: In the formula, Q is the heat source density inside the system, ρ is the material density, c is the material specific heat capacity, k is the material thermal conductivity, and T is the temperature of any point on the thin plate; The temperature field control equation after hybrid unit discretization is as follows: For the temperature gradient, determine the temperature gradient control equation after the mixed unit is discretized: In the formula, q1, q2, and q3 represent three thermodynamic boundary conditions: given boundary temperature, given boundary heat flux, and boundary as heat exchange boundary. N is the generalized coordinate of temperature, q M is the generalized coordinate of temperature gradient, z is the thickness of the thin plate, and S is the element shape function; The matrix form of the temperature field control equation is as follows: Where C is the heat capacity matrix, K c is the thermal conductivity matrix, Q T is the heat source matrix, q T is the generalized coordinate matrix of temperature field; The electrostatic field control equations are discretized based on the hybrid unit. The specific contents are as follows: The governing equations for the electrostatic field are in the following form: Where ρ is the charge density, ε is the dielectric constant, is the electric potential; Through the functional of Poisson's equation and discretization using mixed units, the following equation form is obtained: In the formula, ε r is the relative dielectric constant; Taking the functional extreme value and sorting it out, we get the discretized matrix form of the electrostatic field control equation:

6. The simulation method for calculating the rigid-flexible-electric-thermal coupled dynamic response of a spatial structure according to claim 1, characterized in that: The specific method of step 3 is: to achieve bidirectional coupling between temperature field and displacement field by introducing radiation angle coefficient. The specific process is as follows: Elastic energy caused by thermal stress The expression is as follows: Where E is the elastic coefficient matrix of the thin plate, ε T for thermal strain; The elastic energy caused by thermal stress is derived from the generalized coordinates to obtain the generalized temperature load, which is specifically expressed as: Introducing the radiation angle coefficient Γ α , the specific form is as follows: Where s is the unit vector of solar radiation heat flux; n is the surface normal vector of the thin plate unit; Radiation angle coefficient Γ α The solar radiation is linked to the generalized nodal coordinates of the displacement field to achieve bidirectional coupling between the temperature field and the displacement field; The bidirectional coupling between the electrostatic field and the displacement field is realized through the piezoelectric constitutive equation. The specific process is as follows: Based on the constitutive equation of piezoelectric materials, the specific coupling relationship between the displacement field and the electrostatic field is derived; the constitutive equation of piezoelectric materials is expressed as: Where σ is the element stress vector, ε g is the unit strain vector, c E is the elasticity matrix, is the electric field intensity, D is the electric displacement, e is the piezoelectric constant matrix, is the dielectric constant matrix; The strain energy is expressed as: Substitute the constitutive equation of the piezoelectric material into the expression of strain energy, and then differentiate the strain energy of the central rigid body-cantilever plate physical model with respect to the displacement vector q to obtain the generalized elastic force: In the formula, ε g is the unit strain, c E is the elasticity matrix, is the electric field intensity, e is the piezoelectric constant matrix; The electric field energy of the electrostatic field can be expressed as follows: The electric field energy of the central rigid body-cantilever plate physical model is converted to the electric potential Taking the derivatives respectively, we get the electrostatic field charge: In the formula, Q c is the electric field charge; When there is no external input voltage, the electric field charge Q c is zero, and solving the above equation gives the piezoelectric potential generated by the structural deformation.

7. The simulation method for calculating the rigid-flexible-electrical-thermal coupled dynamic response of a spatial structure according to claim 6, characterized in that: The dynamic equations, temperature field control equations and electrostatic field control equations of the central rigid body-cantilever plate physical model are linked through the bidirectional coupling effect, and then the unit equations are systematically assembled based on the finite element method, and finally the rigid-flexible-electric-thermal coupling system control equations of the central rigid body-cantilever plate physical model are obtained. The specific process is as follows: Through the assembly of finite element units, the rigid-flexible-electric-thermal coupling system control equations of the multi-body system based on the hybrid unit are given. The rigid-flexible-electric-thermal coupling system control equations consist of three parts: dynamic equations, temperature field control equations and electrostatic field control equations. The equation form is as follows: Where M qq is the system mass matrix; K qq ,K qφ ,K qT ,F q They are the system stiffness matrix, force-electric coupling matrix, force-thermal coupling matrix and system external force matrix respectively; K φφ ,F φ are the dielectric stiffness matrix and charge load column vector respectively; C, K TT ,F T They are the system heat capacity matrix, heat conduction matrix and heat flux load column vector respectively; Φ(q,t),Φ(q T ,t),Φ(q φ ,t) are system constraints, temperature field boundary conditions and electrostatic field boundary conditions respectively.

8. The simulation method for calculating the rigid-flexible-electrical-thermal coupled dynamic response of a spatial structure according to claim 1, characterized in that: The specific process of step 4 is as follows: Based on the generalized α method, a multi-physics hybrid solver is constructed. The multi-physics hybrid solver is used to solve the system dynamics equations and process the temperature field control equations and the electrostatic field control equations at the same time. The specific method is as follows: The system dynamics equation is iterated and calculated as follows: The system dynamics equation is shown below: Where M qq is the system mass matrix; K qq ,K qφ ,K qT ,F q They are the system stiffness matrix, force-electric coupling matrix, force-thermal coupling matrix and system external force matrix, q,q T, q e are the displacement coordinate array, temperature coordinate array and potential coordinate array respectively, λ is the Lagrange multiplier, Φ(q,t) is the system constraint; To simplify the expression, remember: Q = K qq q+K qφ q e +K qT q T -F q First, the following parameters are introduced and the initial values ​​of the iteration at time t = 0 are given, including the generalized displacement q0, the generalized velocity Generalized acceleration And the initial Lagrange multiplier λ0, the specific form is as follows: After a time of t = n × h, the next time step t = t is estimated using the dynamic response of the previous moment. n+1 Q at the moment n+1 , λ n+1 ,a n+1 The specific discrete format is as follows: In the formula, q n , are the derivatives of the displacement coordinates of the unit node and the hybrid unit node with respect to time at the previous time step; q n+1 , is the unit displacement node coordinates of the current time step and the derivative of the unit node displacement coordinates with respect to time; a n ,a n+1 are the acceleration of the previous time step and the acceleration of the current time step respectively; h is the integration time step of the generalized-α method; α m ,α f ,β,γ are the numerical solution parameters of the generalized-α method; Q n+1 , λ n+1 Substituting into the kinetic equation we get: Solve using Newton-Raphson iteration: In the formula, for In the kth step, iterate on q n+1 The Jacobi matrix of is the iterative increment; Introducing correction factors into the generalized-α method The corrected k+1 iteration results are as follows: Substitute the result after k+1 iterations and calculate like This means that the accuracy of the iteration result meets the set requirements; As the initial condition for the next time step, and update the parameter a n+1 Start the next time step loop iteration until the numerical simulation time ends; The temperature field control equation is iteratively calculated, and the specific contents are as follows: The specific forms of the temperature field control equation and its constraints are as follows: In the formula, C, K TT ,F T are the system heat capacity matrix, heat conduction matrix and heat flux load column vector, Φ(q T ,t) is the boundary condition of temperature field; The temperature q of the temperature field is given T , temperature change rate Lagrange multiplier λ T The initial value of the iteration is as follows: After a time of t = n × h, the next time step t = t is estimated using the dynamic response of the previous moment. n+1 Q at the moment T(n+1) , λ T(n+1) , the specific discrete format is as follows: The estimated iteration value is the next iteration value q T(n+1) , λ T(n+1) The initial value of q T(n+1) , λ T(n+1) Substitute into the temperature field control equation and we get get: In the formula, For the kth time right The Jacobi matrix of Iterative increment; Introducing correction factors The corrected k+1 iteration results are as follows: Substitute the result after k+1 iterations and calculate like This means that the accuracy of the iteration result meets the set requirements; As the initial condition of the next time step, the next time step loop iteration begins until the numerical simulation time ends; the iteration process of the electrostatic field control equation is consistent with the iteration process of the temperature field control equation; Finally, the displacement, velocity, acceleration dynamic response results of the central rigid body-cantilever plate model and the model temperature and electric potential physical information are obtained through iterative calculation.

Images