wing torsion box calculation implemented in python

Question

Is my wing torsion box calculation implemented in python correct? It would be great if someone would also kindly check if my calculation for torsion of a wing box is correctly done. I posted this on the aviation stack exchange but they somehow deemed it off topic, which is understandable but it's still aeronautically structures related. Maybe because it's because of the Python coding.

Anyways, for this, I have two main files: the torsion_box_calc.py file which contains all the necessary formulas and the torsion_test.py file where I add some dummy variables to test if the calculations work. I also added a constraint in the torsion_test.py file which states that if the wing tip deformation due to torsion is more than 4°, then it is above the limit which I imposed.

torsion_box_calc.py:

import numpy as np
class TorsionBox:
 def __init__(self, 
 enclosed_area: float, 
 segment_lengths: np.ndarray, 
 thickness: np.ndarray, 
 shear_modulus: float) -> None:
 self.enclosed_area = enclosed_area
 self.segment_lengths = segment_lengths
 self.thickness = thickness
 self.shear_modulus = shear_modulus
 self.torque = None 
 self.J_t = None
 def calculate_torsional_constant(self) -> float:
 """
 Using Bredt's formula to find J_t in m^4.
 But first, we should probs check if segments given is equal to thickness given.
 Gonna use mah boi raise ValueError here because it might continue calculating shiz.
 """
 if len(self.segment_lengths) != len(self.thickness):
 raise ValueError("The number of segment lengths and number of thickness values do not match!")
 perimeter_sum = np.sum(self.segment_lengths / self.thickness)
 self.J_t = (4 * (self.enclosed_area ** 2)) / perimeter_sum
 return self.J_t
 
 def calculate_twist_rate(self) -> float:
 """
 For d_theta / d_x based on given self.torque by baby Bredt
 """
 if self.J_t is None:
 self.calculate_torsional_constant()
 
 twist_rate = self.torque / (self.shear_modulus * self.J_t)
 return twist_rate
 
 def calculate_twist_angle(self, wing_length: float) -> float:
 twist_rate = self.calculate_twist_rate()
 twist_angle = twist_rate * wing_length
 return twist_angle
 
 def calculate_aerodynamic_torque(self,
 rho: float,
 velocity: float, 
 coefficient_of_lift: float,
 wing_span: float, 
 aerodynamic_center: float,
 center_of_pressure: float, 
 chord_length: float, 
 n_zp: float) -> float:
 area = wing_span * chord_length
 lift_force = 0.5 * rho * (velocity ** 2) * area * coefficient_of_lift * n_zp
 self.torque = lift_force * abs(center_of_pressure - aerodynamic_center)
 return self.torque

torsion_test.py:

import numpy as np
import json
from torsion_box_calc import TorsionBox
from data_transfer import dataObject
data_filename = "data.json"
with open(data_filename, 'r') as file:
 data_file = json.load(file)
data = dataObject(data_file)
def main():
 # From Performance Data
 rho = data.performance.rho.value
 #From Structures Data
 enclosed_area = 0.5 # m^2
 segment_lengths = np.array([0.025, 0.01, 0.025, 0.007]) # m
 thickness = np.array([0.006, 0.006, 0.006, 0.006]) # m
 shear_modulus = 5e9 # Pa CFRP UD and Fabric
 wing_span = 1.7 # m --> or from aerodynamics: sqrt(AR*area)
 chord_length = 0.12 # m
 wing_length = wing_span/2 # m
 n_zp = 6 # pre-defined in RFP; the positive maximum load factor
 #From Aerodynamics Data
 coefficient_of_lift = data.aerodynamics.CL_opt.value
 speed = 25 # m/s
 aerodynamic_center = chord_length / 4
 center_of_pressure = 0.15 # m
 # Creation of TorsionBox
 torsion_box = TorsionBox(enclosed_area, segment_lengths, thickness, shear_modulus)
 # Calculate aerodynamic torque
 torque = torsion_box.calculate_aerodynamic_torque(rho, speed, coefficient_of_lift, 
 wing_span, aerodynamic_center, 
 center_of_pressure, chord_length, n_zp)
 # Calculate torsional constant
 J_t = torsion_box.calculate_torsional_constant()
 print(f"Torsional constant (J_t): {J_t} m^4")
 # Calculate twist rate
 twist_rate = torsion_box.calculate_twist_rate()*(180/np.pi)
 print(f"Twist rate: {twist_rate}°/m")
 # Calculate twist angle
 twist_angle = torsion_box.calculate_twist_angle(wing_length)*(180/np.pi)
 print(f"Twist angle: {twist_angle}°")
 # Check if twist angle exceeds 4°
 max_angle = 4 # °
 if twist_angle > max_angle:
 print("Twist angle EXCEEDS 4°!")
 else:
 print("Twist angle is within the limit.")
if __name__ == "__main__":
 main()

These are the two files. It would be great if anyone had any suggestions or remarks for improvement and correction.

• Reinderien: Can you please edit your question to include sample JSON data? Up to perhaps a ~50-row excerpt.

Answer 1

Yay, good job with type annotation! Thank you. I imagine you're using mypy --strict.

@dataclass

Your box has half a dozen parameters, and we somewhat tediously assign self.x = x, self.y = y, .... This makes it a good candidate for being a dataclass. And then we get the constructor "for free"!

citations

It's always important to cite your references. I imagine you found several formulas in a text on your bookshelf. But it's not on my bookshelf and I have no idea what you're describing here. After consulting the global default for such things, I only found Bredt's rule which seems non-responsive and is probably about some different guy. Absent reference material, I am unable to "kindly check if [the] calculation for torsion of a wing box is correctly done."

Strike that third docstring line. Yes, it is absolutely bog standard to raise ValueError() when caller has supplied inappropriate inputs. We wouldn't want to trust some bogus intermediate result and keep calculating with it.

side effects

This is slightly weird:

 self.J_t = ...
 return self.J_t

I would rather see the caller do the assignment:

 if self.J_t is None:
 self.J_t = self.calculate_torsional_constant()

We have several calculate_FOO() helpers, and the usual expectation would be that they are "pure" functions, without any side effects that we should also be thinking about.

Even better would be to handle such defaulting up in the constructor. If caller supplying a zero value for J_t is not physically meaningful, then this would suffice:

 self.J_t = J_t or self.calculate_torsional_constant()

units

return twist_angle

I imagine we're getting some number of radians back, rather than degrees. But I wouldn't mind seeing a """docstring""" spell that out. No one likes to see things like Mars Climate Orbiter encountering a planet "at a lower-than-anticipated altitude".

All maintenance engineers who interact with this code over the course of its lifetime will need to pull a text from your bookshelf, or consult an internet resource, or read the source code, to understand details such as the units being used. As written the OP code fails to make this crystal clear.

 shear_modulus = 5e9 # Pa CFRP UD and Fabric

Ooohh, look at that. Wonderful! We apply SI units to carbon-fiber reinforced polymer uni-directional braid.

But I would much prefer to see the constructor pin down such details.

Similarly, I never guessed that calculate_twist_rate() produces degrees per meter, though looking back I suppose I might have done dimensional analysis of its caller.

__main__ guard

At bottom of the test module we see the conventional invocation of main(), which is very nice. But earlier the import of that module attempted to execute

with open(data_filename, 'r') as file:
 data_file = json.load(file)

If I cd someplace which does not have a readable data.json file, then the import fails, which is very regrettable. We call this a "side effect", raising FileNotFoundError during import. A codebase littered with such annoyances is one which will prove difficult to recruit team members to and to maintain.

Prefer to put such logic within the guard, perhaps within def main().

assertions

 print(f"Twist angle: {twist_angle}°")

To make this module a self-documenting automated unit test, please self.assertEqual() or assert ... == ... so we know that twist_angle is e.g. three degrees for those specified inputs. As written, the OP code merely says "tada!" without telling us whether the computed figure is known to be accurate.

---

It this code correct? The OP code does not reveal enough details for me to pronounce it either correct or incorrect.