Status: Please note thet this article is still a work in progress. The reason why it's published before completion is that to prevent me from getting into the perfectionism trap and dragging the article for months!! Please bear with me.
Updates to expect:
- Creating a loft path and using it for lofting,
- Similar approach for VAWT - Darrieus
from __future__ import print_function, division import numpy as np import matplotlib.pyplot as plt from matplotlib.style import use %matplotlib inline from IPython.core.display import Image use("ggplot")
Finite Element (FE) is an indispensable tool for structural analysis but there are few caveats, one of them is the model development time. Especially for complex geometries like wind turbine blades, where the Aerodynamic surfaces are based on non-planar surfaces, Getting a valid geometry to work with is still a challenge. Even with all the bi-directional support between CAD packages(CATIA, SolidWorks, ..) and modeling capabilities of CAD tools!
There are tonnes of journal articles which gives the bird's eye view of developing Abaqus model for wind turbine modeling, but for starters, it's time-consuming to figure out. I've faced a similar situation when I started my master thesis, so I thought it'll be helpful to others if I write it down.
Although the aerodynamic shape of the wind turbine blade is complex, complete geometry can be described with fairly simple topology As shown in Figure 1. Like region specific airfoil profile, distance from the root, pitch axis location, twist angle, chord length, etc. These 2D profiles can be drawn as wires in Abaqus and then the complete blade geometry can be created by lofting the successive 2D profiles.
The chord and alignment of airfoil section is fairly straight forward; chord is nothing but scaling of the airfoil in both the axis, whereas alignment of airfoil with pitch axis is subtraction along the x-axis. Twisting of the airfoil should be done about both the axis (2D rotation), which is nothing but a product of rotation matrix and coordinates .
The location of negative - in sin(theta) depends the direction of rotation which we prefer, for clockwise rotation the element (2, 1) is negative.
def transformCoords(airfoilFile, twist=0.0, chord=1.0, PA_Loc=0.25, **kwargs): """ Returns the twisted Coordinates Parameters ---------- airfoilFile, str: file name which contains the coordinates of airfoil used, twist, float: twist in degrees, default=0 chord, float: chord in metres, default=1 (unit chord) rad_Loc, float: radial Location of the airfoil in blade, default=0 PA_Loc, float: Pitch axis location in percent of chord, default=0.25, Quarter-Chord location. (0 to 1) Returns ---------- nd array: Coordinates rotated about pitch axis location """ rotationMatrix = lambda theta: np.array([ [np.cos(theta), np.sin(theta)], [-np.sin(theta), np.cos(theta)] ]) twist = np.deg2rad(twist) coords = np.loadtxt(airfoilFile, dtype=float, **kwargs) PA_Loc *= chord # % P.A*chord coords[:, 0] -= PA_Loc # Translate the x/c coordinates to the pitch axis twistedCoords = np.zeros((coords.shape + 1, 2)) coords_dash = np.dot(rotationMatrix(twist), coords.T).T x_dash, y_dash = coords_dash[:, 0], coords_dash[:, 1] twistedCoords[:-1, 0], twistedCoords[:-1, 1] = x_dash, y_dash #Join the first and last nodes to form a closed profile, very important for Abaqus twistedCoords[-1, 0], twistedCoords[-1, 1] = x_dash, y_dash twistedCoords *= chord # Scale x/c and y/c to the chord return twistedCoords
fig, ax = plt.subplots() for angle in range(-20, 30, 10): twistedCoords_1 = transformCoords("S822.txt", twist=angle, chord=1.0, PA_Loc=0.25, skiprows=2) ax.plot(twistedCoords_1[:, 0], twistedCoords_1[:, 1], label="%2d" % angle); ax.legend(loc=2)
The implementation of the above mentioned 2D rotation and an example with NREL S822 airfoil profile is shown in Figure 2.
In order to draw the airfoil in Abaqus, it needs one more piece of information which is z location of the point. In our case z is the radial location of the cross-section. One more thing to keep in mind is that Abaqus expects the each data point of coordinates to be in tuple type. Which can be done easily with list comprehension.
coordTuple = tuple([coordinates[i] for i in range(coordinates.shape)])
Another important thing to verify is that no two points are same (i.e., the distance between two consecutive points should be greater than 1E-16_). Yes, I removed the repeating (1, 0) in S822 file.
airfoilFileName = "S822.txt" twist, chord, radLoc, paLoc = 10.0, 1.0, 0.5, 0.25 # Arbitary information coordinates = transformCoords(airfoilFileName, twist, chord, paLoc, skiprows=2) tempCoord_storage = np.zeros((coordinates.shape, 3)) # Add radLoc to the coords tempCoord_storage[:, :2] = coordinates tempCoord_storage[:, 2] = radLoc # Each coord dp as tuple coordTuple_i = tuple([tempCoord_storage[i] for i in range(tempCoord_storage.shape)])
# Only in Abaqus CLI or Abaqus PDE or Executing script inside Abaqus from abaqus import * from abaqusConstants import * partObject = mdb.models['Model-1'].Part(name="wt_blade") # Abaqus part object partObject.WirePolyLine(coordTuple_i)
After creating a new part object, just passing this coordTuple_i to Abaqus part object will do the magic of drawing the profile as wire! Yeah, its simple as that. Now we can do all sort of pythonic stuff like,
* Providing a master blade file with the necessary information for each cross-section and iterating, * Renaming wires to a more sensible name on the fly and so on and so forth.
Now these cross-sectional profiles can be lofted manually or automatically using part.ShellLoft(loftsections=...), depends on what level of automation the designer wishes. Automatic lofting is a bit of overkill unless some sort of optimisation which involves parameters that change the blade profile.