From ccd3712f5e565ef71e3919ddd190a740d041a70d Mon Sep 17 00:00:00 2001 From: GiordanoNesci <153726969+GiordanoNesci@users.noreply.github.com> Date: Fri, 26 Jan 2024 15:35:29 +0100 Subject: [PATCH 1/7] Aero disturbance files --- aero_disturbance_test.py | 108 ++++++++++++++++++++++++++++++++ disturbance_calculations.py | 120 ++++++++++++++++++++++++++++++++++++ 2 files changed, 228 insertions(+) create mode 100644 aero_disturbance_test.py create mode 100644 disturbance_calculations.py diff --git a/aero_disturbance_test.py b/aero_disturbance_test.py new file mode 100644 index 0000000..11816c8 --- /dev/null +++ b/aero_disturbance_test.py @@ -0,0 +1,108 @@ +from paseos.attitude.disturbance_calculations import calculate_aero_torque +from paseos.utils.reference_frame_transfer import transformation_matrix_rpy_body +import trimesh +import numpy as np +""" +Getting the data to run the function +""" + + +""" Actual tests: + +Test zero: get the Cd of the single plates and weigh it against the Cd often used in missions: Cd=2 +Done modifying the function to print Cd value. + +Verified for a set 3 different altitudes and satellite's attitudes +Cd values obtained in [0.02 and 2.31], low Cd values correspond to low values of the angle of attack +""" + +""" +First test uses existing spacecraft data to check that with the same initial conditions the torque exerted on the +satellite is comparable to the real data. For example GOCE data reveal that the aerodynamic torque levels are in the +order of 10^-3 - 10^-4 for a axially symmetric satellite of 5x2 m. We can expect torques in the order of 10^-4 +in magnitude for a cuboid satellite 1x1 m at the same altitude as GOCE (255 km). +""" +vertices = [ +[-0.5, -0.5, -0.5], +[-0.5, -0.5, 0.5], +[-0.5, 0.5, -0.5], +[-0.5, 0.5, 0.5], +[0.5, -0.5, -0.5], +[0.5, -0.5, 0.5], +[0.5, 0.5, -0.5], +[0.5, 0.5, 0.5], +] +faces = [ +[0, 1, 3], +[0, 3, 2], +[0, 2, 6], +[0, 6, 4], +[1, 5, 3], +[3, 5, 7], +[2, 3, 7], +[2, 7, 6], +[4, 6, 7], +[4, 7, 5], +[0, 4, 1], +[1, 4, 5], +] +mesh = trimesh.Trimesh(vertices, faces) +# GOCE orbital altitude and speed +r = np.array([6626, 0, 0]) #6626 +v = np.array([0, 7756, 0]) +temperature_t = 1000 + +transformation_matrix = transformation_matrix_rpy_body(np.array([0, np.pi/4, np.pi/4])) + +T = calculate_aero_torque(r, v, mesh, temperature_t, transformation_matrix) +abs_T = np.linalg.norm(T)*10**3 + +# print("Total torque is ", abs_T, "mN") +# print("The expected torque is in the order of 10^-4") + +if 10**-4 < abs_T/10**3 < 10**-3: + print("Test 1 passed") +else: + print("Test 1 not passed") + + +""" +Second test aims to model simple cases and verify that the torque exerted matches the expected torque calculated by the +function, with particula emphasis on the torque vector direction. +""" + +# Base case: RPY frame is aligned with Body Frame. There should be no torque. +transformation_matrix = transformation_matrix_rpy_body(np.array([0, 0, 0])) # can change the euler angles if willing +T = calculate_aero_torque(r, v, mesh, temperature_t, transformation_matrix) +if np.linalg.norm(T) < 10**(-4): + print("Test 2.1 passed") +else: + print("Test 2.1 not passed") + +# Cases of simple rotations +# First case, satellite with an attitude that yields a torque vector parallel to [0, -1, 0]. +transformation_matrix = transformation_matrix_rpy_body(np.array([0, np.pi/6, 0])) +T = calculate_aero_torque(r, v, mesh, temperature_t, transformation_matrix) +similarity = np.dot(T,[0, 1, 0]) +if similarity == np.linalg.norm(T): + print("Test 2.2 passed") +else: + print("Test 2.2 not passed") + +# Second case, satellite with an attitude that yields zero torque +transformation_matrix = transformation_matrix_rpy_body(np.array([np.pi/6, 0, 0])) +T = calculate_aero_torque(r, v, mesh, temperature_t, transformation_matrix) +if all(T) == 0: + print("Test 2.3 passed") +else: + print("Test 2.3 not passed") + +# Third case, positive rotation of 30° about the z-axis, expected torque in positive z-axis direction +transformation_matrix = transformation_matrix_rpy_body(np.array([0, 0, np.pi/6])) +T = calculate_aero_torque(r, v, mesh, temperature_t, transformation_matrix) +similarity = np.dot(T, [0, 0, 1]) + +if similarity == np.linalg.norm(T): + print("Test 2.4 passed") +else: + print("Test 2.4 not passed") diff --git a/disturbance_calculations.py b/disturbance_calculations.py new file mode 100644 index 0000000..5df56f2 --- /dev/null +++ b/disturbance_calculations.py @@ -0,0 +1,120 @@ +# this functions could be implemented inside the attitude model class but I'd rather have it +# separately right now. + +# OUTPUT NUMPY ARRAYS +import numpy as np +import trimesh + +from paseos.utils.reference_frame_transfer import transformation_matrix_eci_rpy, transformation_matrix_rpy_body, \ + eci_to_rpy + + + +def calculate_aero_torque(r,v,mesh,temperature_spacecraft, transformation_matrix_rpy_body,): + """ Calculates the aerodynamic torque on the satellite. Possible further inputs are temperature_spacecraft from the + thermal model. + Returns: + T: torque vector in the spacecraft body frame + """ + + # Set data useful for aerodynamic coefficients calculation + temperature_spacecraft = 300 # [K] we could connect this with the thermal model + temperature_gas = 1000 # [K] valid for upper stratosphere, could be completely wrong higher up + altitude = np.linalg.norm(r)-6371 # [km] + density = 10**(-(altitude+1285)/151) # from data fit from thermospheric mass density + R = 8.314462 # [J/(K mol)] Universal Gas Constant + molecular_speed_ratio_t = np.linalg.norm(v)/np.sqrt(2*R*temperature_gas) # Used in the Cd and Cl calculation + molecular_speed_ratio_r = np.linalg.norm(v)/np.sqrt(2*R*temperature_spacecraft) # Used in the Cd and Cl calculation + accommodation_parameter = 0.85 + + # Get the normal vectors of all the faces of the mesh in the spacecraft body reference frame, and then they get + # translated in the Roll Pitch Yaw frame with a transformation from paseos.utils.reference_frame_transfer.py + face_normals_sbf = mesh.face_normals[0:12] + face_normals_rpy = np.dot(transformation_matrix_rpy_body, face_normals_sbf.T).T # Check transposition + + # Get the velocity and transform it in the Roll Pitch Yaw frame. Get the unit vector associated with the latter + v_rpy = eci_to_rpy(v, r, v) + unit_v_rpy = v_rpy/np.linalg.norm(v_rpy) + + # Get the normals, the areas, the angle of attack and the centroids of the faces with airflow, confronting them + # with the velocity vector direction. If the dot product between velocity and normal is positive, the plate + # corresponding to the normal is receiving airflow. + + # Initialize parameters and variables for the for loop. + j = 0 + normals_faces_with_airflow = np.zeros((6, 3)) # Maximum six faces have airflow + alpha = [0, 0, 0, 0, 0, 0] # angle of attack matrix, stores the angles of attack in radiant + area_all_faces_mesh = mesh.area_faces # ordered list of all faces' areas + area_faces_airflow = [0, 0, 0, 0, 0, 0] + centroids_faces_airflow = [0, 0, 0, 0, 0, 0] + + for i in range(12): + if np.dot(face_normals_rpy[i], v_rpy) > 0: + normals_faces_with_airflow[j] = face_normals_rpy[i] + area_faces_airflow[j] = area_all_faces_mesh[i] # get the area of the plate [i] which is receiving airflow + face_vertices = mesh.vertices[mesh.faces[j]] + centroids_faces_airflow[j] = np.mean(face_vertices, axis=0) # get the centroids of the face[i] with airflow + alpha[j] = np.arccos(np.dot(normals_faces_with_airflow[j], unit_v_rpy)) + j += 1 + + # Get the aerodynamic coefficient Cd and Cl for every plate and calculate the corresponding drag and lift forces + # for every plate [k] with airflow. Calculate the torque associated with the forces with the previously calculated + # centroid position of every plate. + + # Initialization of the variables + force_drag = [0, 0, 0, 0, 0, 0] + force_lift = [0, 0, 0, 0, 0, 0] + torque_aero = [[0, 0, 0], + [0, 0, 0], + [0, 0, 0], + [0, 0, 0], + [0, 0, 0], + [0, 0, 0]] + alpha = np.array(alpha) + + for k in range(j): + alpha_scalar = alpha[k] + C_d = 1/molecular_speed_ratio_t**2 * (molecular_speed_ratio_t / np.sqrt(np.pi) * \ + (4 * np.sin(alpha_scalar)**2 + 2 * accommodation_parameter * np.cos(2*alpha_scalar)) * \ + np.exp(-(molecular_speed_ratio_t*np.sin(alpha_scalar))**2) + np.sin(alpha_scalar) * \ + (1 + 2 * molecular_speed_ratio_t ** 2 + (1 - accommodation_parameter) * \ + (1 - 2 * molecular_speed_ratio_t ** 2 * np.cos(2 * alpha_scalar))) * \ + np.math.erf(molecular_speed_ratio_t * np.sin(alpha_scalar)) + accommodation_parameter * np.sqrt(np.pi)* \ + molecular_speed_ratio_t ** 2 / molecular_speed_ratio_r * np.sin(alpha_scalar) ** 2) + + C_l = np.cos(alpha[k]) / molecular_speed_ratio_t**2 * (2 / np.sqrt(np.pi) * (2-2*accommodation_parameter) * \ + molecular_speed_ratio_t*np.sin(alpha_scalar)*np.exp(-(molecular_speed_ratio_t*np.sin(alpha_scalar))**2)+\ + (2*(2-2*accommodation_parameter)*(molecular_speed_ratio_t*np.sin(alpha_scalar))**2 + \ + (2-accommodation_parameter)) * np.math.erf(molecular_speed_ratio_t*np.sin(alpha_scalar)) + \ + accommodation_parameter*np.sqrt(np.pi) * molecular_speed_ratio_t**2 / molecular_speed_ratio_r * \ + np.sin(alpha_scalar)) + + # Drag force on the plate [k]. Direction along the velocity vector. + force_drag[k] = - 0.5 * density * C_d * area_faces_airflow[k] * np.linalg.norm(v)**2 * unit_v_rpy + # Lift force on the plate [k]. Direction along the (v x n) x v vector. Intermediate step to get v x n. + v_x_n_vector = np.cross(unit_v_rpy, normals_faces_with_airflow[k]) + force_lift[k] = - 0.5 * density * C_l * area_faces_airflow[k] * np.linalg.norm(v)**2 * np.cross(v_x_n_vector, unit_v_rpy) + + # Torque calculated as the product between the distance of the centroid from the geometrica center of the + # satellite and the sum of the forces. + torque_aero[k] = np.cross(centroids_faces_airflow[k]-mesh.centroid, (np.array(force_drag[k]) + np.array(force_lift[k]))) + + # Compute aero torque in the timestep as the vector sum of all the torques + torque_aero = np.array(torque_aero) + T_rpy = torque_aero.sum(axis=0) + T = transformation_matrix_rpy_body @ T_rpy + return np.array(T) + + +def calculate_grav_torque(): + # calculations for torques + # T must be in actor body fixed frame (to be discussed) + T = [0, 0, 0] + return np.array(T) + + +def calculate_magnetic_torque(): + # calculations for torques + # T must be in actor body fixed frame (to be discussed) + T = [0, 0, 0] + return np.array(T) From aacaa42b3e76adaa8a8827d8e6cc77ee63b042e0 Mon Sep 17 00:00:00 2001 From: GiordanoNesci <153726969+GiordanoNesci@users.noreply.github.com> Date: Fri, 26 Jan 2024 15:38:39 +0100 Subject: [PATCH 2/7] Update and rename disturbance_calculations.py to paseos/disturbance_calculations.py --- disturbance_calculations.py => paseos/disturbance_calculations.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename disturbance_calculations.py => paseos/disturbance_calculations.py (100%) diff --git a/disturbance_calculations.py b/paseos/disturbance_calculations.py similarity index 100% rename from disturbance_calculations.py rename to paseos/disturbance_calculations.py From 93f63320633f8acb14f1da8852a7dae7b9f59fdb Mon Sep 17 00:00:00 2001 From: GiordanoNesci <153726969+GiordanoNesci@users.noreply.github.com> Date: Fri, 26 Jan 2024 15:40:02 +0100 Subject: [PATCH 3/7] Delete paseos/attitude/disturbance_calculations.py --- paseos/attitude/disturbance_calculations.py | 26 --------------------- 1 file changed, 26 deletions(-) delete mode 100644 paseos/attitude/disturbance_calculations.py diff --git a/paseos/attitude/disturbance_calculations.py b/paseos/attitude/disturbance_calculations.py deleted file mode 100644 index 43f6fe2..0000000 --- a/paseos/attitude/disturbance_calculations.py +++ /dev/null @@ -1,26 +0,0 @@ -"""this file contains functions that return attitude disturbance torque vectors expressed in the actor body frame""" - -# OUTPUT NUMPY ARRAYS -# OUTPUT IN BODY FRAME -import numpy as np - - -def calculate_aero_torque(): - # calculations for torques - # T must be in actor body fixed frame (to be discussed) - T = [0, 0, 0] - return np.array(T) - - -def calculate_grav_torque(): - # calculations for torques - # T must be in actor body fixed frame (to be discussed) - T = [0, 0, 0] - return np.array(T) - - -def calculate_magnetic_torque(): - # calculations for torques - # T must be in actor body fixed frame (to be discussed) - T = [0, 0, 0] - return np.array(T) From 3ce7dae63afa267d6caa03143f9457f628458ec8 Mon Sep 17 00:00:00 2001 From: GiordanoNesci <153726969+GiordanoNesci@users.noreply.github.com> Date: Fri, 26 Jan 2024 15:40:22 +0100 Subject: [PATCH 4/7] Update and rename paseos/disturbance_calculations.py to paseos/attitude/disturbance_calculations.py --- paseos/{ => attitude}/disturbance_calculations.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename paseos/{ => attitude}/disturbance_calculations.py (100%) diff --git a/paseos/disturbance_calculations.py b/paseos/attitude/disturbance_calculations.py similarity index 100% rename from paseos/disturbance_calculations.py rename to paseos/attitude/disturbance_calculations.py From b42eac97eeefd9fe28053d64d568a99af6a9435e Mon Sep 17 00:00:00 2001 From: GiordanoNesci <153726969+GiordanoNesci@users.noreply.github.com> Date: Fri, 26 Jan 2024 15:40:46 +0100 Subject: [PATCH 5/7] Update and rename aero_disturbance_test.py to paseos/attitude/aero_disturbance_test.py --- .../attitude/aero_disturbance_test.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename aero_disturbance_test.py => paseos/attitude/aero_disturbance_test.py (100%) diff --git a/aero_disturbance_test.py b/paseos/attitude/aero_disturbance_test.py similarity index 100% rename from aero_disturbance_test.py rename to paseos/attitude/aero_disturbance_test.py From f1c30379ede47117a0746b7c810aadc8dfa07d8a Mon Sep 17 00:00:00 2001 From: GiordanoNesci <153726969+GiordanoNesci@users.noreply.github.com> Date: Wed, 31 Jan 2024 20:30:34 +0100 Subject: [PATCH 6/7] Update aero_disturbance_test.py Reviewed tests and rewritten in pytest --- paseos/attitude/aero_disturbance_test.py | 158 ++++++++++++----------- 1 file changed, 81 insertions(+), 77 deletions(-) diff --git a/paseos/attitude/aero_disturbance_test.py b/paseos/attitude/aero_disturbance_test.py index 11816c8..5eae11d 100644 --- a/paseos/attitude/aero_disturbance_test.py +++ b/paseos/attitude/aero_disturbance_test.py @@ -2,26 +2,9 @@ from paseos.utils.reference_frame_transfer import transformation_matrix_rpy_body import trimesh import numpy as np -""" -Getting the data to run the function -""" -""" Actual tests: - -Test zero: get the Cd of the single plates and weigh it against the Cd often used in missions: Cd=2 -Done modifying the function to print Cd value. - -Verified for a set 3 different altitudes and satellite's attitudes -Cd values obtained in [0.02 and 2.31], low Cd values correspond to low values of the angle of attack -""" - -""" -First test uses existing spacecraft data to check that with the same initial conditions the torque exerted on the -satellite is comparable to the real data. For example GOCE data reveal that the aerodynamic torque levels are in the -order of 10^-3 - 10^-4 for a axially symmetric satellite of 5x2 m. We can expect torques in the order of 10^-4 -in magnitude for a cuboid satellite 1x1 m at the same altitude as GOCE (255 km). -""" +# create mesh vertices = [ [-0.5, -0.5, -0.5], [-0.5, -0.5, 0.5], @@ -47,62 +30,83 @@ [1, 4, 5], ] mesh = trimesh.Trimesh(vertices, faces) -# GOCE orbital altitude and speed -r = np.array([6626, 0, 0]) #6626 -v = np.array([0, 7756, 0]) -temperature_t = 1000 - -transformation_matrix = transformation_matrix_rpy_body(np.array([0, np.pi/4, np.pi/4])) - -T = calculate_aero_torque(r, v, mesh, temperature_t, transformation_matrix) -abs_T = np.linalg.norm(T)*10**3 - -# print("Total torque is ", abs_T, "mN") -# print("The expected torque is in the order of 10^-4") - -if 10**-4 < abs_T/10**3 < 10**-3: - print("Test 1 passed") -else: - print("Test 1 not passed") - - -""" -Second test aims to model simple cases and verify that the torque exerted matches the expected torque calculated by the -function, with particula emphasis on the torque vector direction. -""" - -# Base case: RPY frame is aligned with Body Frame. There should be no torque. -transformation_matrix = transformation_matrix_rpy_body(np.array([0, 0, 0])) # can change the euler angles if willing -T = calculate_aero_torque(r, v, mesh, temperature_t, transformation_matrix) -if np.linalg.norm(T) < 10**(-4): - print("Test 2.1 passed") -else: - print("Test 2.1 not passed") - -# Cases of simple rotations -# First case, satellite with an attitude that yields a torque vector parallel to [0, -1, 0]. -transformation_matrix = transformation_matrix_rpy_body(np.array([0, np.pi/6, 0])) -T = calculate_aero_torque(r, v, mesh, temperature_t, transformation_matrix) -similarity = np.dot(T,[0, 1, 0]) -if similarity == np.linalg.norm(T): - print("Test 2.2 passed") -else: - print("Test 2.2 not passed") - -# Second case, satellite with an attitude that yields zero torque -transformation_matrix = transformation_matrix_rpy_body(np.array([np.pi/6, 0, 0])) -T = calculate_aero_torque(r, v, mesh, temperature_t, transformation_matrix) -if all(T) == 0: - print("Test 2.3 passed") -else: - print("Test 2.3 not passed") - -# Third case, positive rotation of 30° about the z-axis, expected torque in positive z-axis direction -transformation_matrix = transformation_matrix_rpy_body(np.array([0, 0, np.pi/6])) -T = calculate_aero_torque(r, v, mesh, temperature_t, transformation_matrix) -similarity = np.dot(T, [0, 0, 1]) - -if similarity == np.linalg.norm(T): - print("Test 2.4 passed") -else: - print("Test 2.4 not passed") + + +def test_1(): + """ First test uses existing spacecraft data to check that with the same initial conditions the torque exerted on the + satellite is comparable to the real GOCE data. Goce has a cross-section of 1.1 m^2 and during it's drag free + phase experienced between 1 and 20 mN of drag force. It is reasonable to expect a torque in the same order of + magnitude for a satellite of 1x1x1 m^2 at the same altitude as GOCE. + """ + # GOCE orbital altitude and speed + r = np.array([6626, 0, 0]) # [km] + v = np.array([0, 7756, 0]) # [m/s] + + transformation_matrix = transformation_matrix_rpy_body(np.array([np.pi/4, np.pi/6, 0])) + + T = calculate_aero_torque(r, v, mesh, transformation_matrix) + abs_T = np.linalg.norm(T) + + assert 5*10**-4 < abs_T < 2*10**-2 + +def test_2a(): + """ Second tests aim to model simple cases and verify that the torque exerted matches the expected torque calculated + by the function, with particular emphasis on the torque vector direction. + + Base case: RPY frame is aligned with Body Frame. There should be no torque. + """ + # Circular orbit at GOCE altitude + r = np.array([6626, 0, 0]) # [km] + v = np.array([0, 7756, 0]) # [km] + + transformation_matrix = transformation_matrix_rpy_body(np.array([0, 0, 0])) # can change the euler angles if willing + T = calculate_aero_torque(r, v, mesh, transformation_matrix) + + assert np.linalg.norm(T) < 10**(-8) + +def test_2b(): + """ Second tests aim to model simple cases and verify that the torque exerted matches the expected torque calculated + by the function, with particular emphasis on the torque vector direction. + + Second case, satellite with an attitude that yields a torque vector parallel to [0, -1, 0] + """ + # Circular orbit at GOCE altitude + r = np.array([6626, 0, 0]) # [km] + v = np.array([0, 7756, 0]) # [km] + + transformation_matrix = transformation_matrix_rpy_body(np.array([0, np.pi/6, 0])) + T = calculate_aero_torque(r, v, mesh, transformation_matrix) + similarity = np.dot(T,[0, 1, 0]) + + assert similarity == np.linalg.norm(T) + +def test_2c(): + """ Second tests aim to model simple cases and verify that the torque exerted matches the expected torque calculated + by the function, with particular emphasis on the torque vector direction. + + Third case, satellite with an attitude that yields zero torque + """ + # Circular orbit at GOCE altitude + r = np.array([6626, 0, 0]) # [km] + v = np.array([0, 7756, 0]) # [km] + + transformation_matrix = transformation_matrix_rpy_body(np.array([np.pi/6, 0, 0])) + T = calculate_aero_torque(r, v, mesh, transformation_matrix) + + assert all(T) == 0 + +def test_2d(): + """ Second tests aim to model simple cases and verify that the torque exerted matches the expected torque calculated + by the function, with particular emphasis on the torque vector direction. + + Third case, positive rotation of 30° about the z-axis, expected torque in the positive z-axis direction. + """ + # Circular orbit at GOCE altitude + r = np.array([6626, 0, 0]) # [km] + v = np.array([0, 7756, 0]) # [km] + + transformation_matrix = transformation_matrix_rpy_body(np.array([0, 0, np.pi/6])) + T = calculate_aero_torque(r, v, mesh, transformation_matrix) + similarity = np.dot(T, [0, 0, 1]) + + assert similarity == np.linalg.norm(T) From ed4828f82733704bedf178479b09e4f8be606eb1 Mon Sep 17 00:00:00 2001 From: GiordanoNesci <153726969+GiordanoNesci@users.noreply.github.com> Date: Wed, 31 Jan 2024 20:33:40 +0100 Subject: [PATCH 7/7] Update disturbance_calculations.py Fixed version of the aero disturbances (there was a j instead of an i in a for loop) --- paseos/attitude/disturbance_calculations.py | 64 ++++++++++++++------- 1 file changed, 42 insertions(+), 22 deletions(-) diff --git a/paseos/attitude/disturbance_calculations.py b/paseos/attitude/disturbance_calculations.py index 5df56f2..397a666 100644 --- a/paseos/attitude/disturbance_calculations.py +++ b/paseos/attitude/disturbance_calculations.py @@ -1,7 +1,3 @@ -# this functions could be implemented inside the attitude model class but I'd rather have it -# separately right now. - -# OUTPUT NUMPY ARRAYS import numpy as np import trimesh @@ -10,19 +6,34 @@ -def calculate_aero_torque(r,v,mesh,temperature_spacecraft, transformation_matrix_rpy_body,): - """ Calculates the aerodynamic torque on the satellite. Possible further inputs are temperature_spacecraft from the - thermal model. +def calculate_aero_torque(r,v,mesh, transformation_matrix_rpy_body): + """ Calculates the aerodynamic torque on the satellite. + The model used is taken from "Roto-Translational Spacecraft Formation Control Using Aerodynamic Forces"; Ran. S, + Jihe W., et al.; 2017. The mass density of the atmosphere is calculated from the best linear fit of the data + extracted from "Thermospheric mass density: A review"; Emmert J.T.; 2015. This function only works for Earth + centered orbits and satellites defined by a cuboid mesh as in the geometrical model. + + Future addition to the model can be the temperature of the spacecraft as an input parameter. At the moment the + temperature of the spacecraft is assumed in the model. The sensitivity of the torque with respect to the + temperature of the spacecraft and of the gas is low. + + Args: + r (np.array): distance from the satellite and the Earth's center of mass [km] + v (np.array): velocity of the satellite in ECI reference frame [m/s] + mesh (trimesh): mesh of the satellite from the geometric model + transformation_matrix_rpy_body (np.array): transformation matrix from Roll-Pitch_Yaw frame to the Spacecraft + Body frame. From reference_frame_transfer.py . + Returns: - T: torque vector in the spacecraft body frame + T (np.array): torque vector in the spacecraft body frame """ - # Set data useful for aerodynamic coefficients calculation - temperature_spacecraft = 300 # [K] we could connect this with the thermal model - temperature_gas = 1000 # [K] valid for upper stratosphere, could be completely wrong higher up + # Constants for aerodynamic coefficients calculation + temperature_spacecraft = 300 # [K] + temperature_gas = 1000 # [K] + R = 8.314462 # Universal Gas Constant, [J/(K mol)] altitude = np.linalg.norm(r)-6371 # [km] - density = 10**(-(altitude+1285)/151) # from data fit from thermospheric mass density - R = 8.314462 # [J/(K mol)] Universal Gas Constant + density = 10**(-(altitude+1285)/151) # equation describing the best linear fit for the data, [kg/m^3] molecular_speed_ratio_t = np.linalg.norm(v)/np.sqrt(2*R*temperature_gas) # Used in the Cd and Cl calculation molecular_speed_ratio_r = np.linalg.norm(v)/np.sqrt(2*R*temperature_spacecraft) # Used in the Cd and Cl calculation accommodation_parameter = 0.85 @@ -30,15 +41,18 @@ def calculate_aero_torque(r,v,mesh,temperature_spacecraft, transformation_matrix # Get the normal vectors of all the faces of the mesh in the spacecraft body reference frame, and then they get # translated in the Roll Pitch Yaw frame with a transformation from paseos.utils.reference_frame_transfer.py face_normals_sbf = mesh.face_normals[0:12] - face_normals_rpy = np.dot(transformation_matrix_rpy_body, face_normals_sbf.T).T # Check transposition + face_normals_rpy = np.dot(transformation_matrix_rpy_body, face_normals_sbf.T).T # Get the velocity and transform it in the Roll Pitch Yaw frame. Get the unit vector associated with the latter v_rpy = eci_to_rpy(v, r, v) unit_v_rpy = v_rpy/np.linalg.norm(v_rpy) - # Get the normals, the areas, the angle of attack and the centroids of the faces with airflow, confronting them - # with the velocity vector direction. If the dot product between velocity and normal is positive, the plate - # corresponding to the normal is receiving airflow. + # Loop to get the normals, the areas, the angle of attack and the centroids of the faces with airflow, confronting them + # with the velocity vector direction. + # If the dot product between velocity and normal is positive, the plate corresponding to the normal is receiving + # airflow and is stored in the variable normals_faces_with_airflow. Similarly, the areas of the faces is stored in + # area_faces_airflow, and the angles of attack of the respective faces in alpha. The last three lines calculate the + # centroids of the faces that receive airflow, for later use in the torque calculation. # Initialize parameters and variables for the for loop. j = 0 @@ -52,9 +66,10 @@ def calculate_aero_torque(r,v,mesh,temperature_spacecraft, transformation_matrix if np.dot(face_normals_rpy[i], v_rpy) > 0: normals_faces_with_airflow[j] = face_normals_rpy[i] area_faces_airflow[j] = area_all_faces_mesh[i] # get the area of the plate [i] which is receiving airflow - face_vertices = mesh.vertices[mesh.faces[j]] - centroids_faces_airflow[j] = np.mean(face_vertices, axis=0) # get the centroids of the face[i] with airflow alpha[j] = np.arccos(np.dot(normals_faces_with_airflow[j], unit_v_rpy)) + face_vertices = mesh.vertices[mesh.faces[i]] + face_vertices_rpy = np.dot(np.linalg.inv(transformation_matrix_rpy_body), face_vertices.T).T + centroids_faces_airflow[j] = np.mean(face_vertices_rpy, axis=0) # get the centroids of the face[i] with airflow j += 1 # Get the aerodynamic coefficient Cd and Cl for every plate and calculate the corresponding drag and lift forces @@ -91,11 +106,14 @@ def calculate_aero_torque(r,v,mesh,temperature_spacecraft, transformation_matrix # Drag force on the plate [k]. Direction along the velocity vector. force_drag[k] = - 0.5 * density * C_d * area_faces_airflow[k] * np.linalg.norm(v)**2 * unit_v_rpy - # Lift force on the plate [k]. Direction along the (v x n) x v vector. Intermediate step to get v x n. + # Lift force on the plate [k]. Direction along the (v x n) x v direction, lift vector defined to be in that + # direction. Intermediate step to get v x n. v_x_n_vector = np.cross(unit_v_rpy, normals_faces_with_airflow[k]) - force_lift[k] = - 0.5 * density * C_l * area_faces_airflow[k] * np.linalg.norm(v)**2 * np.cross(v_x_n_vector, unit_v_rpy) + not_norm_lift_vector = np.cross(v_x_n_vector, unit_v_rpy) + lift_vector = not_norm_lift_vector/np.linalg.norm(not_norm_lift_vector) + force_lift[k] = - 0.5 * density * C_l * area_faces_airflow[k] * np.linalg.norm(v)**2 * lift_vector - # Torque calculated as the product between the distance of the centroid from the geometrica center of the + # Torque calculated as the product between the distance of the centroid from the geometric center of the # satellite and the sum of the forces. torque_aero[k] = np.cross(centroids_faces_airflow[k]-mesh.centroid, (np.array(force_drag[k]) + np.array(force_lift[k]))) @@ -103,6 +121,8 @@ def calculate_aero_torque(r,v,mesh,temperature_spacecraft, transformation_matrix torque_aero = np.array(torque_aero) T_rpy = torque_aero.sum(axis=0) T = transformation_matrix_rpy_body @ T_rpy + T[np.isnan(T)] = 0 # substitutes 0 to any NaN value + return np.array(T)