symcad.parts.composite.YFormAirfoils

  1#!/usr/bin/env python3
  2# Copyright (C) 2022, Will Hedgecock
  3#
  4# This program is free software: you can redistribute it and/or modify
  5# it under the terms of the GNU General Public License as published by
  6# the Free Software Foundation, either version 3 of the License, or
  7# (at your option) any later version.
  8#
  9# This program is distributed in the hope that it will be useful,
 10# but WITHOUT ANY WARRANTY; without even the implied warranty of
 11# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 12# GNU General Public License for more details.
 13#
 14# You should have received a copy of the GNU General Public License
 15# along with this program.  If not, see <http://www.gnu.org/licenses/>.
 16
 17from __future__ import annotations
 18from PyFreeCAD.FreeCAD import FreeCAD, Part
 19from typing import Dict, Optional, Tuple, Union
 20from sympy import Expr, Symbol, sqrt, sin, cos
 21from . import CompositeShape
 22import math
 23
 24class YFormAirfoils(CompositeShape):
 25   """Model representing a set of Y-form parameteric airfoils.
 26
 27   By default, the airfoils are oriented in the following configuration:
 28
 29   ![YFormAirfoils](https://symbench.github.io/SymCAD/images/YFormAirfoils.png)
 30
 31   The `geometry` of this shape includes the following parameters:
 32
 33   - `max_thickness`: Maximum thickness (in `% of length`) of each airfoil
 34   - `chord_length`: Length (in `m`) of the chord of each airfoil
 35   - `span`: Width (in `m`) of each airfoil
 36   - `separation_radius`: Radius (in `m`) of separation between each airfoil
 37   - `curvature_tilt`: Amount of tilt (in `deg`) of each airfoil
 38
 39   Note that the above dimensions should be interpreted as if the airfoils are unrotated. In other
 40   words, any shape rotation takes place *after* the airfoil dimensions have been specified.
 41   """
 42
 43   # Constructor ----------------------------------------------------------------------------------
 44
 45   def __init__(self, identifier: str, material_density_kg_m3: Optional[float] = 1.0) -> None:
 46      """Initializes a Y-form parametric airfoil object.
 47
 48      Parameters
 49      ----------
 50      identifier : `str`
 51         Unique identifying name for the object.
 52      material_density_kg_m3 : `float`, optional, default=1.0
 53         Uniform material density in `kg/m^3` to be used in mass property calculations.
 54      """
 55      super().__init__(identifier, self.__create_cad__, None, material_density_kg_m3)
 56      setattr(self.geometry, 'max_thickness', Symbol(self.name + '_max_thickness'))
 57      setattr(self.geometry, 'chord_length', Symbol(self.name + '_chord_length'))
 58      setattr(self.geometry, 'span', Symbol(self.name + '_span'))
 59      setattr(self.geometry, 'separation_radius', Symbol(self.name + '_separation_radius'))
 60      setattr(self.geometry, 'curvature_tilt', Symbol(self.name + '_curvature_tilt'))
 61
 62
 63   # CAD generation function ----------------------------------------------------------------------
 64
 65   @staticmethod
 66   def z_val(x: float, max_thickness: float, chord_length: float) -> float:
 67      return 5.0 * max_thickness * chord_length * (
 68         (0.2969 * x**0.5) - (0.1260 * x) - (0.3516 * x**2) + (0.2843 * x**3) - (0.1036 * x**4))
 69
 70   @staticmethod
 71   def __create_cad__(params: Dict[str, float], _fully_displace: bool) -> Part.Solid:
 72      """Scripted CAD generation method for `YFormAirfoils`."""
 73      doc = FreeCAD.newDocument('Temp')
 74      max_thickness_percent = params['max_thickness']
 75      chord_length_mm = 1000.0 * params['chord_length']
 76      span_mm = 1000.0 * params['span']
 77      separation_radius_mm = 1000.0 * params['separation_radius']
 78      curvature_tilt = params['curvature_tilt']
 79      x, points_upper, points_lower = 0.0, [], []
 80      while x <= 1.009:
 81         points_upper.append(FreeCAD.Vector(x * chord_length_mm,
 82                                            YFormAirfoils.z_val(x, max_thickness_percent,
 83                                                                   chord_length_mm)))
 84         points_lower.append(FreeCAD.Vector(x * chord_length_mm,
 85                                            -YFormAirfoils.z_val(x, max_thickness_percent,
 86                                                                    chord_length_mm)))
 87         x += 0.01
 88      bodies = []
 89      for i in range(3):
 90         body = doc.addObject('PartDesign::Body','Airfoil' + str(i))
 91         sketch = doc.addObject('Sketcher::SketchObject', 'Sketch' + str(i))
 92         sketch.Support = doc.XZ_Plane
 93         sketch.MapMode = 'FlatFace'
 94         body.addObject(sketch)
 95         sketch.addGeometry(Part.BSplineCurve(points_upper, None, None,
 96                                              False, 3, None, False), False)
 97         sketch.addGeometry(Part.BSplineCurve(points_lower, None, None,
 98                                              False, 3, None, False), False)
 99         pad = doc.addObject('PartDesign::Pad', 'Pad' + str(i))
100         body.addObject(pad)
101         pad.Length = int(span_mm)
102         pad.Profile = sketch
103         y_offset = sin(math.pi / 3.0) * separation_radius_mm
104         z_offset = cos(math.pi / 3.0) * separation_radius_mm
105         placement_vector = \
106            FreeCAD.Vector(0, y_offset if i == 1 else (-y_offset if i == 2 else 0),
107                              -z_offset if i == 1 or i == 2 else separation_radius_mm)
108         body.Placement = FreeCAD.Placement(placement_vector,
109                          FreeCAD.Rotation(-curvature_tilt if i == 1 else
110                                           (curvature_tilt if i == 2 else 0),
111                                           curvature_tilt if i == 0 else 0, -90.0 - (i * 120.0)))
112         bodies.append(body)
113      fusion = doc.addObject('Part::MultiFuse', 'Fusion')
114      fusion.Shapes = bodies
115      doc.recompute()
116      airfoils = fusion.Shape
117      FreeCAD.closeDocument(doc.Name)
118      return airfoils
119
120
121   # Geometry setter ------------------------------------------------------------------------------
122
123   def set_geometry(self, *, max_thickness_percent: Union[float, None],
124                             chord_length_m: Union[float, None],
125                             span_m: Union[float, None],
126                             separation_radius_m: Union[float, None],
127                             curvature_tilt_deg: Union[float, None]) -> YFormAirfoils:
128      """Sets the physical geometry of the current `YFormAirfoils` object.
129
130      See the `YFormAirfoils` class documentation for a description of each geometric
131      parameter.
132      """
133      self.geometry.set(max_thickness=max_thickness_percent,
134                        chord_length=chord_length_m,
135                        span=span_m,
136                        separation_radius=separation_radius_m,
137                        curvature_tilt=curvature_tilt_deg)
138      return self
139
140   def get_geometric_parameter_bounds(self, parameter: str) -> Tuple[float, float]:
141      parameter_bounds = {
142         'max_thickness': (0.01, 0.90),
143         'chord_length': (0.1, 2.0),
144         'span': (0.0, 2.0),
145         'separation_radius': (0.0, 1.5),
146         'curvature_tilt': (0.0, 45.0)
147      }
148      return parameter_bounds.get(parameter, (0.0, 0.0))
149
150
151   # Geometric properties -------------------------------------------------------------------------
152
153   @property
154   def material_volume(self) -> Union[float, Expr]:
155      return self.displaced_volume
156
157   @property
158   def displaced_volume(self) -> Union[float, Expr]:
159      x, points_upper = 0.0, []
160      while x <= 1.009:
161         points_upper.append((x * self.geometry.chord_length,
162                              YFormAirfoils.z_val(x, self.geometry.max_thickness,
163                                                     self.geometry.chord_length)))
164         x += 0.01
165      area = 0.0
166      for i in range(1, len(points_upper) - 1):
167         area += points_upper[i][1]
168      area = points_upper[0][1] + (2.0 * area) + points_upper[len(points_upper)-1][1]
169      h = (points_upper[len(points_upper)-1][0] - points_upper[0][0]) / (len(points_upper) - 1)
170      return 3.0 * h * area * self.geometry.span
171
172   @property
173   def surface_area(self) -> Union[float, Expr]:
174      x, points_upper = 0.0, []
175      while x <= 1.009:
176         points_upper.append((x * self.geometry.chord_length,
177                              YFormAirfoils.z_val(x, self.geometry.max_thickness,
178                                                     self.geometry.chord_length)))
179         x += 0.01
180      area = 0.0
181      for i in range(1, len(points_upper)):
182         x_diff = (points_upper[i][0] - points_upper[i-1][0])
183         y_diff = (points_upper[i][1] - points_upper[i-1][1])
184         area += sqrt(x_diff**2 + y_diff**2)
185      area = (2.0 * area) * self.geometry.span
186      return 3.0 * area
187
188   @property
189   def unoriented_center_of_gravity(self) -> Tuple[Union[float, Expr],
190                                                   Union[float, Expr],
191                                                   Union[float, Expr]]:
192      cosine = cos(self.geometry.curvature_tilt * math.pi / 180.0)
193      sine = sin(self.geometry.curvature_tilt * math.pi / 180.0)
194      return ((0.41 * cosine * self.geometry.chord_length) + (0.5 * sine * self.geometry.span),
195              0.5 * self.unoriented_width,
196              self.unoriented_height / 3.0)
197
198   @property
199   def unoriented_center_of_buoyancy(self) -> Tuple[Union[float, Expr],
200                                                    Union[float, Expr],
201                                                    Union[float, Expr]]:
202      return self.unoriented_center_of_gravity
203
204   @property
205   def unoriented_length(self) -> Union[float, Expr]:
206      cosine = cos(self.geometry.curvature_tilt * math.pi / 180.0)
207      sine = sin(self.geometry.curvature_tilt * math.pi / 180.0)
208      return (cosine * self.geometry.chord_length) + (sine * self.geometry.span)
209
210   @property
211   def unoriented_width(self) -> Union[float, Expr]:
212      cosine = cos(self.geometry.curvature_tilt * math.pi / 180.0)
213      y_offset = sin(math.pi / 3.0) * self.geometry.separation_radius
214      y_width = cos(math.pi / 6.0) * (cosine * self.geometry.span)
215      return 2.0 * (y_offset + y_width)
216
217   @property
218   def unoriented_height(self) -> Union[float, Expr]:
219      cosine = cos(self.geometry.curvature_tilt * math.pi / 180.0)
220      z_offset = cos(math.pi / 3.0) * self.geometry.separation_radius
221      z_width = sin(math.pi / 6.0) * (cosine * self.geometry.span)
222      return z_offset + z_width + self.geometry.separation_radius + self.geometry.span
223
224   @property
225   def oriented_length(self) -> Union[float, Expr]:
226      # TODO: Implement this
227      return 0
228
229   @property
230   def oriented_width(self) -> Union[float, Expr]:
231      # TODO: Implement this
232      return 0
233
234   @property
235   def oriented_height(self) -> Union[float, Expr]:
236      # TODO: Implement this
237      return 0