In this notebook we investigate the effect of varying the nonlinear analysis parameters of MSC Natran SOL 106 on the load-displacement diagram of a box beam structure under a bending load.
import matplotlib.pyplot as plt # package for plotting
import tol_colors as tc # package for colorblind-friendly colors
plt.rc('axes', prop_cycle=plt.cycler('color', list(tc.tol_cset('bright'))))
FIRST_SUBCASE_ID = 1
We consider the same box beam geometry and material properties of our previous notebook.
AR = 9. # aspect ratio - 2*b/w (the length of the box beam corresponds to half the span of the CRM wing)
w = 1e3 # width [mm]
l = AR*w/2 # length [mm]
non_dimensional_height = 0.2 # h/w
h = w*non_dimensional_height # box height [mm]
non_dimensional_thickness = 1/50 # t/h
t = h*non_dimensional_thickness # shell thickness [mm]
print(f"""
Box beam dimensions:
- width: {w/1e3:.1f} m
- length: {l/1e3:.1f} m
- height: {h/1e3:.1f} m
- wall thickness: {t:.1f} mm""")
Box beam dimensions: - width: 1.0 m - length: 4.5 m - height: 0.2 m - wall thickness: 4.0 mm
rho = 2780e-12 # density [ton/mm^3]
E = 73.1e3 # Young's modulus [MPa]
nu = 0.3 # Poisson's ratio
We discretize our box beam using the function mesh_box_beam
from the box_beam_utils
module. We choose the target length of the shell elements based on the results of the mesh convergence study of our previous notebook. Then we use the function create_base_bdf_input
from the same module to create a base BDF
object including material and element properties, nodes, quadrilateral shell elements, boundary conditions and some settings for the output files. Finally, we visualize the cards of our bdf input with the method get_bdf_stats
.
from resources import box_beam_utils # module with useful functions to work with box beam problems
shell_element_length = 59.9 # [mm]
box_beam_mesh = box_beam_utils.mesh_box_beam(width=w, length=l, height=h, element_length=shell_element_length)
nodes_coordinates_array = box_beam_mesh.points
nodes_connectivity_matrix = box_beam_mesh.faces.reshape(-1, 5)[:, 1:]
box_beam_bdf_input = box_beam_utils.create_base_bdf_input(
young_modulus=E, poisson_ratio=nu, density=rho, shell_thickness=t, nodes_xyz_array=nodes_coordinates_array,
nodes_connectivity_matrix=nodes_connectivity_matrix) # create base BDF object
print(box_beam_bdf_input.get_bdf_stats())
---BDF Statistics--- SOL None bdf.spcs[1]: 1 SPC1: 1 bdf.params: 0 PARAM : 1 bdf.nodes: 0 GRID : 3388 bdf.elements: 0 CQUAD4 : 3344 bdf.properties: 0 PSHELL : 1 bdf.materials: 0 MAT1 : 1
subcase=0 already exists...skipping
To complete the setup of our numerical model, we need to define a load. We consider a concentrated bending load applied at the tip section, which is rigidly constrained. We choose this configuration among the different ones investigated in our last notebook because it's the one that showed the most similar response to a global beam-like behavior.
We define a master node at the center of the tip section and we rigidly constrain all the nodes along the tip section to the master node with a RBE2
element. Then we apply a unitary force along the $z$-axis at the master node using a FORCE
card.
import numpy as np # package to work with arrays and linear algebra
# Add RBE2 element to make tip section rigid
MASTER_NODE_ID = np.size(nodes_coordinates_array, 0) + 1 # define id of master node of tip section
box_beam_bdf_input.add_grid(MASTER_NODE_ID, [w/2, l, 0.]) # add GRID card to define master node of tip section
nodes_ids = np.arange(1, np.size(nodes_coordinates_array, 0) + 1) # find the id of all nodes
tip_nodes_ids = nodes_ids[nodes_coordinates_array[:,1]==l] # find id of the nodes located at the tip section
rbe2_eid = len(box_beam_bdf_input.elements) + 1 # define id of the RBE2 element
box_beam_bdf_input.add_rbe2(rbe2_eid, MASTER_NODE_ID, '123456', tip_nodes_ids) # add RBE2 card
# Add concentrated force
force_set_id = 11
concentrated_force_magnitude = 1.
concentrated_force_direction = [0., 0., 1.]
box_beam_bdf_input.add_force(sid=force_set_id, node=MASTER_NODE_ID, mag=concentrated_force_magnitude,
xyz=concentrated_force_direction)
FORCE 11 3389 1. 0. 0. 1.
Now we want to evaluate the influence of SOL 106's nonlinear analysis parameters on the results of the analysis. We are going to assess this by observing the change of the load-displacement curve of the box beam with the value of each considered parameter. The monitored displacement will be the $z$-displacement of the master node.
We set up the nonlinear analysis with the arc-length method using the function set_up_arc_length_method
from the pynastran_utils
module. We use Nastran's default parameters to run a reference analysis.
from resources import pynastran_utils # module with useful functions to work with pyNastran objects
pynastran_utils.set_up_arc_length_method(box_beam_bdf_input, kstep=-1) # set up SOL 106 with arc-length method setting KSTEP=-1 to prevent the solver from ignoring the displacement error function
We recall the buckling load predicted by SOL 105 from our last notebook and we create a subcase where we apply such load to our box beam. Finally, we run the analysis and we plot both the load-displacement diagram and the structural deformation at the end of the analysis. For the displacement we monitor the displacement along $z$ of the master node of the tip section, $u_{z, tip}$, and we nondimensionalize it with the length of the box beam $l$.
import os # package to work with directories and files
from pyNastran.op2.op2 import read_op2 # function to read Nastran op2 files
from resources import plot_utils # module with useful functions to plot results
# Create subcase where we apply the buckling load predicted by SOL 105
sol_105_buckling_load = 1654. # [N]
load_set_id = force_set_id + 1 # define id of load set
box_beam_bdf_input.add_load(sid=load_set_id, scale=1., scale_factors=[sol_105_buckling_load], load_ids=[force_set_id]) # add LOAD card to define load set
pynastran_utils.create_static_load_subcase(bdf=box_beam_bdf_input, subcase_id=1, load_set_id=load_set_id) # create subcase with defined load set
# Run analysis
analysis_directory_name = '05_Sensitivity_Study_of_SOL_106_Nonlinear_Analysis_Parameters'
analysis_directory_path = os.path.join(os.getcwd(), 'analyses', analysis_directory_name)
input_filename = 'nonlinear_analysis_reference'
pynastran_utils.run_analysis(directory_path=analysis_directory_path, bdf=box_beam_bdf_input,
filename=input_filename, run_flag=False)
# Read op2 file
op2_filepath = os.path.join(analysis_directory_path, input_filename + '.op2')
sol_105_op2 = read_op2(op2_filename=op2_filepath, debug=None, load_geometry=True)
# Find load and displacement history
_, loads, displacements = pynastran_utils.read_load_displacement_history_from_op2(
op2=sol_105_op2, node_ids=[MASTER_NODE_ID])
# Plot load-displacement diagram
plt.rcParams['figure.dpi'] = 120 # default dpi of figures
load_component_index = 2 # select load along z-axis
displacement_component_index = 2 # select displacement along z-axis
_, ax = plt.subplots() # create figure with one subplot
ax.plot(displacements[MASTER_NODE_ID][FIRST_SUBCASE_ID][:,displacement_component_index]/l*100,
loads[FIRST_SUBCASE_ID][:, load_component_index]/sol_105_buckling_load, 'o-') # plot load vs displacement
plt.xlabel('$u_{z, tip}/l$, %')
plt.ylabel('$P/P_{SOL 105}$')
plt.grid()
plt.show()
# Plot structural deformation at the last converged iteration
_, ax, cbar = plot_utils.plot_deformation(op2=sol_105_op2, displacement_component='tz', length_unit="mm")
ax.locator_params(axis='x', nbins=3) # number of x-axis ticks
ax.locator_params(axis='z', nbins=2) # number of z-axis ticks
ax.tick_params(axis='y', which='major', pad=20) # distance of y-axis tick labels
ax.tick_params(axis='z', which='major', pad=6) # distance of z-axis tick labels
ax.yaxis.labelpad = 60 # distance of y-axis label
ax.zaxis.labelpad = 10 # distance of z-axis label
cbar.ax.set_position(cbar.ax.get_position().shrunk(1.0, .66)) # decrease colorbar size
cbar.ax.set_position(cbar.ax.get_position().translated(0, .14)) # move colorbar upwards
plt.show()
The load-displacement diagram appears to be linear, suggesting that the nonlinear analysis parameters should not have any influence on the results for the load range investigated here. However, the final deformation of the structure indicates the presence of a slightly buckled-like shape, as a blob can be observed over the root of the box beam, suggesting that some nonlinear effects may come into play.
The first parameter that we are going to investigate is the error function, defined by the CONV field of the NLPARM
card. The choice of the error functions determines how convergence is checked for every load increment. Three different error functions are available in Nastran:
In the above expressions $| ... |$ represents the vector norm, that is to say $|\boldsymbol{x}|=\sum_{r=1}^n |x_r|$, $\boldsymbol{R}$ represents the vector of out-of-balance forces (residual), $\boldsymbol{\Delta u}$ represents the vector of incremental displacements at the current iteration, $E^*=\left|\boldsymbol{P^\prime}\cdot\boldsymbol{u}\right|$ is the reference energy, with $\boldsymbol{P^\prime}=\boldsymbol{\Delta P} + \boldsymbol{P}$. $\boldsymbol{u}$, and $\boldsymbol{\Delta P}$ are respectively the displacement vector and the incremental load vector at the current iteration, while $\boldsymbol{P}$ is the total load at the previous loading step. Finally, $\boldsymbol{w}$ is a weighting function formed by collecting the square root of the diagonal terms of the stiffness matrix and $q$ is a contraction factor defined as the ratio of displacement increments between two successive steps. For more details please refer to the NX Nastran Handbook of Nonlinear Analysis (Solutions 106 and 129) manual.
We are going to consider four different combinations of error functions:
Vector component checking consists in performing the convergence checking on the maximum vector component of all components in the model.
Let's define the list of the test flags corresponding to the mentioned combinations of functions.
error_functions = ['PW', 'PWV', 'PU', 'PUV', 'PWU', 'PWUV']
Since we need to plot the load-displacement diagram for different values of the nonlinear analysis parameters, we define a list of markers to use in the plots. We also define the function plot_load_displacement_diagram
to run the Nastran analysis, read the load and displacement history from the op2 file and plot the load-displacement curve on a pre-defined figure.
from matplotlib.lines import Line2D # class defining the characters for the marker styles
markers = list(Line2D.markers.keys())[2:] # store list of markers
# Define function to run analysis and plot load-displacement diagram
def plot_load_displacement_diagram(bdf_input, filename, axes, marker_style, line_label):
# Run analysis
pynastran_utils.run_analysis(directory_path=analysis_directory_path, bdf=bdf_input, filename=filename,
run_flag=False)
# Read load and displacement history from op2 file
op2_path = os.path.join(analysis_directory_path, filename + '.op2')
op2 = read_op2(op2_path, debug=None)
_, p, disp = pynastran_utils.read_load_displacement_history_from_op2(
op2=op2, node_ids=[MASTER_NODE_ID])
# Plot load-displacement curve on input axes
axes.plot(disp[MASTER_NODE_ID][FIRST_SUBCASE_ID][:,displacement_component_index]/l*100,
p[FIRST_SUBCASE_ID][:, load_component_index]/sol_105_buckling_load,
marker=marker_style, linestyle='-', label=line_label)
Now we can run the analyses for the different combinations of error functions and plot the results.
nlparm_id = next(iter(box_beam_bdf_input.nlparms)) # retrieve id of NLPARM card
_, ax = plt.subplots() # create figure with one subplot
# Run analysis for each set of error functions and plot load-displacement diagram
for count, fun in enumerate(error_functions):
box_beam_bdf_input.nlparms[nlparm_id].conv = fun
input_filename = "nonlinear_analysis_error_function_" + fun
plot_load_displacement_diagram(box_beam_bdf_input, input_filename, ax, markers[count], f"{fun}")
# Set plot appearance
plt.xlabel("$u_{z,\,tip}/l$, %")
plt.ylabel("$P/P_\mathrm{SOL\,105}$")
plt.legend()
plt.grid()
plt.show()
Nastran job nonlinear_analysis_error_function_PU.bdf completed Wall time: 9.0 s Nastran job nonlinear_analysis_error_function_PUV.bdf completed Wall time: 10.0 s Nastran job nonlinear_analysis_error_function_PWU.bdf completed Wall time: 9.0 s Nastran job nonlinear_analysis_error_function_PWUV.bdf completed Wall time: 9.0 s
The results do not appear to be influenced by the error functions. This suggests that the buckled-like shape observed earlier in the deformation plot develops only at the end of the load range investigated, so no nonlinear behavior can be observed. As a consequence, the response of the structure appears linear and the nonlinear parameters cannot have any effect on the equilibrium path found by the analysis. For this reason, we restore the default load and energy error functions combination.
box_beam_bdf_input.nlparms[nlparm_id].conv = 'PW'
The second parameter that we are going to consider is the convergence tolerance, which corresponds to the value of the error function that determines the convergence of the load increment iteration. The convergence tolerance is defined by the EPSP and the EPSW fields of the NLPARM
card for the load and energy function, respectively.
We are going to consider three different combinations of convergence tolerances, taken from table 21 of MSC Nastran Quick Reference Guide manual on the default tolerances for static nonlinear SOL 106 models :
where the default values are $EPSP=10^{-2}$ and $EPSW=10^{-2}$.
Let's define the list of convergence tolerances, run the analyses and visualize the results.
# Define error tolerances
load_tolerances = [1e-1, 1e-2, 1e-3]
work_tolerances = [1e-1, 1e-3, 1e-7]
# Create figure, run analysis for each couple of convergence tolerances and plot load-displacement diagram
_, ax = plt.subplots()
for count, epsp in enumerate(load_tolerances):
epsw = work_tolerances[count]
box_beam_bdf_input.nlparms[nlparm_id].eps_p = epsp
box_beam_bdf_input.nlparms[nlparm_id].eps_w = epsw
input_filename = f"nonlinear_analysis_load_tolerance_{epsp:.0e}".replace('.','_')
plot_load_displacement_diagram(box_beam_bdf_input, input_filename, ax, markers[count],
f"$EPSP={epsp:.0e},\,EPSW={epsw:.0e}$")
# Set plot appearance
plt.xlabel("$u_{z,\,tip}$, %")
plt.ylabel("$P/P_\mathrm{SOL\,105}$")
plt.legend()
plt.grid()
plt.show()
As expected, also in this case the equilibrium diagram does not show any dependece on the considered parameter, so we switch back to the default values.
box_beam_bdf_input.nlparms[nlparm_id].eps_p = 1e-2
box_beam_bdf_input.nlparms[nlparm_id].eps_w = 1e-2
Now we consider the initial load increment, which is related to the field NINC of the NLPARM
card. When using the arc-length method, as in our case, the initial load increment $\Delta\mu^1$ is defined as:
$$\Delta\mu^1=\frac{1}{NINC},$$
where $NINC$ represents the number of increments in a nonlinear analysis employing the Newton method with constant load increments.
We are going to consider three values for the initial load increment:
Let's define the list of initial load increments, run the analyses and visualize the results.
initial_load_increments = [.5, .1, .01] # define list of initial load increments
# Create figure, run analysis for each initial load increment and plot load-displacement diagram
_, ax = plt.subplots()
for count, delta_mu_1 in enumerate(initial_load_increments):
box_beam_bdf_input.nlparms[nlparm_id].ninc = int(1/delta_mu_1)
input_filename = f"nonlinear_analysis_initial_load_increment_{delta_mu_1:.9g}".replace('.','_')
plot_load_displacement_diagram(box_beam_bdf_input, input_filename, ax, markers[count],
f"$\Delta\mu^1={delta_mu_1:.9g}$")
# Set plot appearance
plt.xlabel("$u_{z,\,tip}/l$, %")
plt.ylabel("$P/P_\mathrm{SOL\,105}$")
plt.legend()
plt.grid()
plt.show()
Nastran job nonlinear_analysis_initial_load_increment_0_01.bdf completed Wall time: 11.0 s
We observe that a smaller initial load increment results in a better resolution of the load-displacement curve, but the equilibrium path is not affected by the change in the initial load increment. As a consequence we restore the default value of the NINC field.
box_beam_bdf_input.nlparms[nlparm_id].ninc = 10
Now we want to vary the desired number of iteration $I_d$, defined by the DESITER field of the NLPCI
card. This parameter controls how the arc-length $\Delta l$ is updated between successive load increments:
$$\Delta l_{new}=\sqrt{\frac{I_d}{I_{max}}}\Delta l_{old},$$
where $I_{max}$ is the number of iterations required for convergence at the preceding step.
We are going to consider three values of desired number of iterations:
Let's define the list of desired number of iterations, run the analyses and visualize the results.
no_iterations = [20, 12, 5] # define list of desired numbers of iterations
# Create figure, run analysis for each desired number of iterations, plot load-displacement diagram
_, ax = plt.subplots()
for count, desiter in enumerate(no_iterations):
box_beam_bdf_input.nlpcis[nlparm_id].desiter = desiter
input_filename = f"nonlinear_analysis_desiter_{desiter:d}"
plot_load_displacement_diagram(box_beam_bdf_input, input_filename, ax, markers[count], f"$\mathrm{{DESITER}}={desiter}$")
# Set plot appearance
plt.xlabel("$u_{z,\,tip}/l$, %")
plt.ylabel("$P/P_\mathrm{SOL\,105}$")
plt.legend()
plt.grid()
plt.show()
We observe that a smaller number of desired number of iterations produces a load-equilibrium curve with more resolution. However, once again the equilibrium path does not appear to have any notable change, so we restore the default value.
box_beam_bdf_input.nlpcis[nlparm_id].desiter = 12
Next, we investigate the effect of the minimum allowable arc-length adjustment ratio, defined by the MINALR field of the NLPCI
card. The arc-length adjustment ratio $\Delta l_{new}/\Delta l_{old}$ is bound in the following way:
$$MINALR<\frac{\Delta l_{new}}{\Delta l_{old}}<MAXALR.$$
We investigate three values of minimum allawable arc-length adjustment ratio:
Let's define the list of values, run the analyses and visualize the results.
min_arclength_adjustment_ratios = [.9, .25, .01] # define list of minimum allowable arc-length adjustment ratios
# Create figure, run analysis for each minimum allowable arc-length adjustment ratio and plot load-displacement diagram
_, ax = plt.subplots()
for count, minalr in enumerate(min_arclength_adjustment_ratios):
box_beam_bdf_input.nlpcis[nlparm_id].minalr = minalr
input_filename = f"nonlinear_analysis_minalr_{minalr:.9g}".replace('.','_')
plot_load_displacement_diagram(box_beam_bdf_input, input_filename, ax, markers[count], f"$\mathrm{{MINALR}}={minalr:.9g}$")
# Set plot appearance
plt.xlabel("$u_{z,\,tip}/l$, %")
plt.ylabel("$P/P_\mathrm{SOL\,105}$")
plt.legend()
plt.grid()
plt.show()
In this case we do not observe see any effect of $MINALR$ on the load-displacement curve, so we switch back to the default value.
box_beam_bdf_input.nlpcis[nlparm_id].minalr = .25
Finally, we consider the maximum allowable arc-length adjustment ratio, defined by the MAXALR field of the NLPCI
card. We consider three values:
max_arclength_adjustment_ratios = [10., 4., 1.01, 1.0001] # define list of maximum allowable arc-length adjustment ratios
# Create figure, run analysis for each maximum allowable arc-length adjustment ratio and plot load-displacement diagram
_, ax = plt.subplots()
for count, maxalr in enumerate(max_arclength_adjustment_ratios):
box_beam_bdf_input.nlpcis[nlparm_id].maxalr = maxalr
input_filename = f"nonlinear_analysis_maxalr_{maxalr:.9g}".replace('.','_')
plot_load_displacement_diagram(box_beam_bdf_input, input_filename, ax, markers[count], f"$\mathrm{{MAXALR}}={maxalr:.9g}$")
# Set plot appearance
plt.xlabel("$u_{z,\,tip}/l$, %")
plt.ylabel("$P/P_\mathrm{SOL\,105}$")
plt.legend()
plt.grid()
plt.show()
Nastran job nonlinear_analysis_maxalr_1_01.bdf completed Wall time: 12.0 s Nastran job nonlinear_analysis_maxalr_1_0001.bdf completed Wall time: 11.0 s
We observe that a lower $MAXALR$ results in a more resolved load-displacement curve, but all curves appear once again seem to describe the same equilibrium path.
In this notebook we have investigated the influence of SOL 106's nonlinear analysis parameters on the load-displacement diagram of a box beam loaded with a concentrated bending load at the tip. We applied a load equal to the buckling load predicted by SOL 105 and we found a linear load-displacement curve. As a consequence the variation of the nonlinear analysis parameters investigated did not show any effect on the curve.
Based on the mentioned evidence, we can conclude that the results obtained with SOL 106's nonlinear buckling method in our last notebook do not depend on the nonlinear analysis parameters. In our next notebook we are going to continue our investigation on the validity of SOL 106's nonlinear buckling method.