If we want to make a shear/moment diagram using FEMAP, we need to make repeated Interface Load resultants. This can be tedious, but is an option. Alternatively, this can be created using pyNastran. PyNastran enables analysts using Nastran to efficiently create, manipulate, and extract data from models. It handles the underlying details so you get models that will run smoothly, without worrying about field formatting in the process.
Challenges:
- Organizing and analyzing large result files
- Post-processing results outside of commercial software
Values:
- Results extraction from binary data file
- Calculation of interface loads across many load cases
- Creation of Shear/Moment diagram
- Customizable post-processing
An Overview of pyNastran
pyNastran is an open source python interface library to the various Nastran file formats (BDF, OP2, OP4) that began development back in 2011 and has had continuous development since. Using the BDF interface, you can read/edit/write Nastran geometry without worrying about field formatting. Many checks are also performed to verify that your model is correct. Using the OP2 interface, you can read large result files quickly and efficiently. Additionally, you can extract a subset of the result data and write OP2 and F06 result files. Rather than replacing FEMAP, Simcenter3D, or PATRAN, the goal of pyNastran is to make engineers more efficient and avoid some of the common pitfalls with model creation and analysis.
In this post, pyNastran will be used in conjunction with FEMAP w/ NX Nastran to run a static trim analysis on the BWB half model, and to post-process results. The following will be explored:
- Extracting interface loads using pyNastran and FEMAP. This will help validate the process.
- Creating a Shear-Moment diagram for the vehicle.
Levering pyNastran effectively enables engineers to automate the extraction of results (e.g., interface loads) and create simplified results visualizations to better understand a model. This understanding is essential to the development of the next generation of vehicles.
Blended Wing Body Model showing Aerodynamic Model
Blended Wing Body Aircraft
The blended wing body (BWB) is a future aircraft concept that has significant aerodynamic benefits. Regardless of whether a freight transport aircraft with overbody engines (see below) or passenger variant with distributed propulsion, the non-tubular design offers additional design challenges. This is fundamentally due to the fact that it’s not a pressure vessel, but rather amore traditional “ship design” (see below) or more exotic concept such as a “double-bubble” design.
This particular BWB is an all-composite concept with additional masses for payload (modeled as non-structural/concentrated masses). Skin thicknesses were pre-optimized under a design 2.5g pullup maneuver load case.
Blended Wing Body – Engine Model and Rigid Elements
The aerodynamic model uses a series aero boxes (shown above) which are splined to the structural model at key locations. For flutter, it’s important to capture the modes. However, for static aero, it’s important to pick robust structural locations (e.g., rib/spar intersections).
Blended Wing Body – Spline Points
What is an Interface Load?
Interface Loads are simply the internal loads on a model. NX/Simcenter Nastran calculates them when the user requests the following as part of the Case Control Deck:
GPFORCE(PLOT) = ALL
This will enable both FEMAP and pyNastran to extract the Grid Point Force table. This table lists the load contributions (forces/moments) on each node and lists what the load is caused by (e.g., element 200, SPC/MPC force, applied load, etc.) and is used to calculate interface loads.
Defining an Interface Load Cut
By exposing the loads (see below), we can expose the internal forces and sum them about a desired point (e.g., the lower left point of the mesh).
Revealing the Interface Loads
This accomplished by (ideally) grabbing all the elements on one side of the cut (including the elements on the cut). We also need the nodes on one side of the cut.
Fastest Way to define an Interface Loads Cut
Alternatively, we can include a few extra nodes to make pulling the nodes a bit easier (at the slight cost of runtime). There are two classifications of “extra nodes”:
- Nodes in front of the cut face not associated with elements (the rightmost nodes): will be removed automatically as there are no loads that affect them. These nodes do not affect runtime as they are quickly filtered.
- Nodes behind the cut face (left side of the cut): all elements that contribute to them are included and thus the loads on those nodes sum to 0. While these nodes have no effect on results, they do affect the runtime.
Method of Creating an Interface Loads cut used in the pyNastran Shear/Moment Plotter
Shear/Moment Diagram Overview
What is simple and straightforward for a 2D model is not surprisingly a bit more complicated in 3D. An important component is the decoupling of the direction that cuts should be taken vs. the orientation of the output coordinate system. Furthermore, it’s not just one output coordinate system; it’s many. We want to march down an axis and sum forces and moments at each station.
We’ll start off by creating a series of 50 coordinate systems about which interface loads should be summed. This is a bit of extra work (for an interface cut). However, we need our summation points anyways. If you know that apriori, it’s less work.
This is done by selecting the starting (p1) and end points (p3). Of importance, the coordinate axis that will be marched down is different than the coordinate axis where loads will be summed. The axial/torque direction of the wing is in the global -y direction. The vertical direction (z) will be consistent with the global to be consistent with the lift direction.
As 50 cutting planes are being considered (cid=110001 to cid=110050), the origin for cid=110001 will be located at p1, while the origin for cid=110050 will be at p3. We’ll discuss exactly how these coordinate systems are determined soon. Before that though, it’s important to note that coordinate system 110028 will be used for validation.
Blended Wing Body – Top View
Blended Wing Body – Side View
Interface Loads
Interface Loads – Model Setup
This is done in pyNastran with the following code. First, we load imports:1. import os
2. import unittest
3. import getpass
4. from pathlib import Path
5.
6. import numpy as np
7. from cpylog import get_logger
8. import matplotlib.pyplot as plt
9.
10. import pyNastran
11. from pyNastran.bdf.bdf import BDF, read_bdf, CORD2R
12. from pyNastran.op2.op2_geom import read_op2_geom
13. from pyNastran.bdf.mesh_utils.cut_model_by_plane import get_element_centroids, get_stations
14. from pyNastran.op2.tables.ogf_gridPointForces.smt import smt_setup, plot_smt
15. from pyNastran.utils.nastran_utils import run_nastran
16. from pyNastran.bdf.utils import parse_patran_syntax_dict 1. PKG_PATH = Path(pyNastran.__path__[0])
2. REL_MODEL_PATH = PKG_PATH / '..' / 'models'
3. MODEL_PATH = REL_MODEL_PATH.resolve()
4. BWB_PATH = MODEL_PATH / 'bwb'
5.
6. op2_filename = BWB_PATH / 'bwb_saero.op2'
7. bdf_filename_out = BWB_PATH / 'bwb_saero_smt.bdf'
1. log = get_logger(level='info')
2. model = read_op2_geom(
3. op2_filename,
4. include_results='grid_point_forces',
5. validate=True, xref=True,
6. build_dataframe=False,
7. log=log)
- nids: node ids in the model (GRIDs and SPOINTs)
- nid_cd: the node id and CD coordinate frame
- xyz_cid0: node locations in the basic coordinate system (coord=0 frame)
- icd_transform: a map of the GRID’s CD/output coordinate frame to node ids
- eids: element ids
- element_centroids_cid0: the element centroids in the basic coordinate system
1. gpforce = model.grid_point_forces[1]
2.
3. (nids, nid_cd, xyz_cid0, icd_transform,
4. eids, element_centroids_cid0) = smt_setup(model) Interface Loads – Extract Results using pyNastran
Relating the BWB model back to the interface loads, the elements that will be included have a centroid with a negative x-value in the output coordinate system (e.g., cid=11028 in the bwb_saero_smt.bdf model, which we’ll write later). The nodes that will be included have an x-location ≥ to -0.2% of the wingspan (-25 inches). This keeps the number of “active” nodes in any given interface loads cut relatively constant.
The coordinate system at the wing-body junction corresponds to coordinate system ID 110028 and is selected for validation.
Selecting Nodes and Elements in the pyNastran GUI using the Highlight Menu
1. data_str = (
2. 'Node 1827 1828 4288 4289 4290 4291 1201827 1202604 1887:1891 1893:1898 2588:2611
4311:4315 1202605:1202611 '
3. 'Elm 1396 1397 1398 1399 1418 1419 1749 1750 1751 1752 2010 2011 2012 2620 2621
2639 2640 2641 1247:1251 1344:1363 1372:1380 1526:1536 1766:1774 1842:1851 2141:2152
2310:2321 2342:2365 2569:2577 2801:2956 3081:3246 3683:3742 3855:3920 4506:4603
4968:5047 5070:5175 5298:5469 5494:5565 5837:5954'
4. )1. nodes_elems_dict = parse_patran_syntax_dict(data_str)
2. cut_nids = nodes_elems_dict['Node']
3. cut_eids = nodes_elems_dict['Elm']
1. origin = np.array([1125.529, 566.6673, 38.76809])
2. zaxis = origin + np.array([0., 0., 1.])
3. xzplane = origin + np.array([0., -1., 0.])
4. summation_point = origin
5.
6. cut_coord = CORD2R(2, origin, zaxis, xzplane) 1. force_sumi, moment_sumi = gpforce.extract_interface_loads(
2. cut_nids, cut_eids,
3. cut_coord, model.coords,
4. nid_cd, icd_transform,
5. xyz_cid0,
6. summation_point=summation_point)
7. force_sumi /= 1000. # lbf -> kips
8. moment_sumi /= 1000. # in-lbf -> in-kips
9. print('force_sum (kips) = ', force_sumi)
10. print('moment_sum (inch-kips) = ', moment_sumi)
force_sum (kips) = [-3.9062502e-6 2.4414064e-7 -126.28]
moment_sum (inch-kips) = [-10558.986 -43460.098 0.03] Note that the model is in English units (inch-lbf-second) and we want the output forces in kips (1 kip=1000 lbf) vs. lbf, and the moments in inch-kips vs. inch-lbf. Thus, we divide the force and moments by 1000.
Interface Loads – Extract Results using FEMAP
Interface Loads can be calculated in FEMAP by selecting a new freebody diagram as shown below by going to the PostProcessing Toolbox menu and Creating a Freebody Diagram. Simply follow the menus to setup the cut.
Creating a Freebody Diagram – Interface Load
The most critical piece is to select the nodes and elements at which to create the cut. Again, we’re selecting coordinate system 110028 for validation.
Creating Interface Loads – Selecting Nodes & Elements
Summed Interface Loads – Forces and Moments
Shear/Moment Diagram
Creating a shear/moment diagram requires the following steps:- Coordinate System Setup
- Creation of the Shear/Moment Diagram
- Writing Geometry for Validation
- Plotting the Shear Moment Diagram
Coordinate System Setup
So where did those coordinate systems come from? Simply define the starting (p1) and end points (p2) to define the march axis. We’re choosing a line that starts from (25% down the wing tip chord/length) to p3 (the projection of the 25% line down the wing root chord). We choose this 25% or ¼ chord line as the lift acts approximately here in subsonic flow. Thus, the torsion of the wing will be more reasonable.1. # defines the starting point for the shear, moment, torque plot
2. p1 = np.array([1388.9, 1262.0, 86.3471]) # 1/4c wing_tip; point A
3.
4. # defines the end point for the shear, moment, torque plot
5. p3 = np.array([910.93, 0.1, 0.0]) # 1/4c wing_root
6.
7. method = 'CORD2R'
8.
9. zaxis = p1 + np.array([0., 0., 1.]) # point B
10. # defines the XZ plane for the shears/moments
11. p2 = p1 + np.array([0., -1., 0.]) # point C
12.
13. xyz1, xyz2, xyz3, i, k, coord_out, iaxis_march, stations = get_stations(
14. model, p1, p2, p3, zaxis,
15. method=method, nplanes=50)
1. nstations = len(stations)
2. coord_origins = np.zeros((nstations, 3))
3. for istation, station in enumerate(stations):
4. coord_origin = coord_out.origin + iaxis_march * station
5. coord_origins[istation, :] = coord_origin
6. ylocations = coord_origins[:, 1]
Creation of the Shear/Moment Diagram
Finally, we can calculate the shear/moment diagram! As discussed previously, we can use a few extra rows of nodes to ensure all the critical nodes are included. This is set as nodes_tol=25.0 inches, which corresponds to about 2 rows of elements. As this model is fairly uniform in side, this is adequate.1. force_sum, moment_sum, new_coords, nelems, nnodes = gpforce.shear_moment_diagram(
2. nids, xyz_cid0, nid_cd, icd_transform,
3. eids, element_centroids_cid0,
4. stations, model.coords, coord_out,
5. iaxis_march=iaxis_march, # defaults to coord_out.i
6. itime=0,
7. nodes_tol=25., log=model.log) Writing Geometry for Validation
We’d like to see where the coordinate systems were actually created in order to validate our model, so we can write those out. We can then view them in the pyNastran GUI or FEMAP:1. for cid, coord in new_coords.items():
2. model.coords[cid] = coord
3. print(model.coords[110028])
4. model.sol = 144 # not saved to the OP2
5. model.write_bdf(bdf_filename_out)
Plotting the Shear/Moment Diagram
Finally, we plot the shear/moment diagram (SMT=shear moment torque).
1. plot_smt(ylocations, force_sum/1000, moment_sum/1000, nelems, nnodes,
2. plot_force_components=False, plot_moment_components=False,
3. show=True,
4. xtitle='Y', xlabel='Spanwise Location, Y (in)',
5. force_unit='kip', moment_unit='in-kip')
Force (1 kip = 1000 lb) is plotted vs. span location (in inches). Aerodynamic load is introduced to the structural model at the spline points and thus, we see large jumps at these locations. This is normal and the aerodynamic moment is correct. Primarily we see this outboard (where stresses/strains are low and minimum gage is typically in effect). Inboard, we see a loss of lift (Z force below 200 inches), as the inboard planform (the top view of the wing) is less efficient. Additionally, there is a large load X Force (y=172 in) corresponding to the engine.
Note: Loads are shown as negative as they’re reactions.
Blended Wing Body – Span Location vs. Force
When looking at the moment, we see a smooth wing bending (Y) moment distribution with the engine-induced torque (X) being clearly visible.
Blended Wing Body – Span Location vs. Moment
For validation purposes, we can also look at the total number of nodes and elements at each station of the cut. This is helpful in making sure we didn’t accidently march downstream (in the direction of the wind) vs. towards the root of the wing. Many shear/moment plots have been ruined this way! If this is constant, you did something very wrong!
Outboard, we see a smooth distribution, while inboard, there are large jumps caused by crossing centerbody ribs that have a large number of nodes/elements.
Blended Wing Body – Span Location vs. Fraction of Included Nodes/Elements
Conclusion
In this post, we showed how to extract interface loads in pyNastran and FEMAP. A shear/moment plot was also created using pyNastran for a BWB configuration, which cannot easily be created using FEMAP. As mentioned before, pyNastran allows engineers to implement and use custom tools in order to perform analyses more efficiently.
The code for this example can be found at: https://github.com/SteveDoyle2/pyNastran/tree/master/models/bwb/shear_moment_torque.ipynb.
If you are interested in a pyNastran blog post on flutter analysis:
https://www.m4-engineering.com/flutter-analysis-with-pynastran/
How to Obtain pyNastran
The GUI is available for download at: https://sourceforge.net/projects/pynastran/files/?source=navbar
The source code can be found at: https://github.com/SteveDoyle2/pyNastran.
Or installed by running:
>>> pip install pyNastran
The documentation can be found at: https://pynastran-git.readthedocs.io/en/1.4/quick_start/index.html.
Note: pyNastran is an open-source project; bug reports, feature requests and help with development is appreciated and encouraged
Written by Steve Doyle
Steve Doyle is a Senior Aerospace Engineer with a background in structural analysis, CFD, and programming. He’s also the author of pyNastran. Outside of dorky projects, he likes to hike and rock climb!
