Simple cantilever for load-controlled sensitivity analysis¶
This document will teach you how to use opstool to post-process sensitivity analysis.
The source code see 14.8.1. Simple cantilever for load-controlled sensitivity analysis.
The numerical model with the associated analysis was described in detail by Prof. Michael Scott within OpenSees Days 2011.
[1]:
import time
import numpy as np
import matplotlib.pyplot as plt
import openseespy.opensees as ops
import opstool as opst
OpenSees Model¶
Use the unit system provided by opstool:
[2]:
UNIT = opst.pre.UnitSystem(length="m", force="kN", time="sec")
[3]:
# Angle
rad = 1.0
deg = np.pi / 180.0 * rad
g = 9.80665 * UNIT.m / (UNIT.sec**2)
Define model¶
[ ]:
# Create ModelBuilder
# -------------------
ops.wipe()
ops.model("basic", "-ndm", 2, "-ndf", 3)
# Create nodes
# ------------
L = 5 * UNIT.m
ops.node(1, 0.0, 0.0) # Fixed end
ops.node(2, L, 0.0) # Free end
# Fixed support
# -------------
ops.fix(1, 1, 1, 1)
# Define material
# ---------------
matTag = 1
Fy = 410 * UNIT.MPa # Yield stress
Es = 200 * UNIT.GPa # Modulus of Elasticity of Steel
b = 2 / 100 # 2% Strain hardening ratio
Hkin = b / (1 - b) * Es
# Sensitivity-ready steel materials: Hardening, Steel01, SteelMP, BoucWen, SteelBRB, StainlessECThermal, SteelECThermal, ...
# Hardening Sensitivity Params: sigmaY/fy/Fy, E, H_kin/Hkin, H_iso/Hiso
ops.uniaxialMaterial("Hardening", matTag, Es, Fy, 0, Hkin)
# ops.uniaxialMaterial("Steel01", matTag, Fy, Es, b) # Sensitivity Params: sigmaY/fy/Fy, E, b, a1, a2, a3, a4
# ops.uniaxialMaterial("SteelMP", matTag, Fy, Es, b) # Sensitivity Params: sigmaY/fy, E, b
# Define sections
# ---------------
# Sections defined with "canned" section ("WFSection2d"), otherwise use a FiberSection object (ops.section("Fiber",...))
beamSecTag = 1
beamWidth, beamDepth = 10 * UNIT.cm, 50 * UNIT.cm
# secTag, matTag, d, tw, bf, tf, Nfw, Nff
ops.section("WFSection2d", beamSecTag, matTag, beamDepth, beamWidth, beamWidth, 0, 20, 0) # Beam section
# Define elements
# ---------------
beamTransTag, beamIntTag = 1, 1
# Linear, PDelta, Corotational
ops.geomTransf("Corotational", beamTransTag)
nip = 5
# Lobatto, Legendre, NewtonCotes, Radau, Trapezoidal, CompositeSimpson
ops.beamIntegration("Legendre", beamIntTag, beamSecTag, nip)
# Beam elements
numEle = 5
meshSize = L / numEle # mesh size
eleType = "dispBeamColumn" # forceBeamColumn, "dispBeamColumn"
# tag, Npts, nodes, type, dofs, size, eleType, transfTag, beamIntTag
ops.mesh("line", 1, 2, *[1, 2], 0, 3, meshSize, eleType, beamTransTag, beamIntTag)
# Create a Plain load pattern with a Sine/Trig TimeSeries
# -------------------------------------------------------
# tag, tStart, tEnd, period, factor
ops.timeSeries("Trig", 1, 0.0, 1.0, 1.0, "-factor", 1.0) # "Sine", "Trig" or "Triangle"
ops.pattern("Plain", 1, 1)
P = 1710 * UNIT.kN
# Create nodal loads at node 2
# nd FX FY MZ
ops.load(2, 0.0, P, 0.0)
Define Sensitivity Parameters¶
[ ]:
# /// Each parameter must be unique in the FE domain, and all parameter tags MUST be numbered sequentially starting from 1! ///
ops.parameter(1) # Blank parameters
ops.parameter(2)
ops.parameter(3)
ops.parameter(4)
for ele in range(1, numEle + 1): # Only column elements
ops.addToParameter(1, "element", ele, "E") # E
# Check the sensitivity parameter names in *.cpp files ("sigmaY" or "fy" or "Fy")
# https://github.com/OpenSees/OpenSees/blob/master/SRC/material/uniaxial/HardeningMaterial.cpp
ops.addToParameter(2, "element", ele, "Fy") # "sigmaY" or "fy" or "Fy"
ops.addToParameter(3, "element", ele, "Hkin") # "H_kin" or "Hkin" or "b"
ops.addToParameter(4, "element", ele, "d") # "d"
ops.parameter(5, "node", 2, "coord", 1) # parameter for coordinate of node 2 in DOF "1" (PX=1, PY=2, MZ=3)
# Map parameter 6 to vertical load at node 2 contained in load pattern 1 (last argument is global DOF, e.g., in 2D PX=1, PY=2, MZ=3)
ops.parameter(6, "loadPattern", 1, "loadAtNode", 2, 2)
# Parameter Symbols and Variables values
ParamSym = {1: "E", 2: "F_y", 3: "H_{kin}", 4: "d", 5: "L", 6: "P"}
ParamVars = {1: Es, 2: Fy, 3: Hkin, 4: beamDepth, 5: L, 6: P}
Define Sensitivity Analysis¶
[ ]:
steps = 500
start_time = time.time()
ops.wipeAnalysis()
ops.system("BandGeneral")
ops.numberer("RCM")
ops.constraints("Transformation")
ops.test("NormDispIncr", 1.0e-12, 10, 3)
ops.algorithm("Newton") # KrylovNewton
ops.integrator("LoadControl", 1 / steps)
ops.analysis("Static")
ops.sensitivityAlgorithm("-computeAtEachStep") # automatically compute sensitivity at the end of each step
Define the output data base and save the data¶
[7]:
ODB = opst.post.CreateODB(odb_tag="sensitivity")
for i in range(steps):
ops.analyze(1)
ODB.fetch_response_step()
ODB.save_response()
print("Sensitvity Analysis Completed!")
print(f"Analysis elapsed time is {(time.time() - start_time):.3f} seconds.\n")
OPSTOOL :: All responses data with odb_tag = sensitivity saved in _OPSTOOL_ODB/RespStepData-sensitivity.nc!
Sensitvity Analysis Completed!
Analysis elapsed time is 2.334 seconds.
Post-processing¶
We want to process the response of node 2 on degree of freedom UY:
[ ]:
ctrlNode = 2
dof = "UY"
Loading node response¶
[9]:
node_resp = opst.post.get_nodal_responses(odb_tag="sensitivity")
print("data variables:", node_resp.data_vars)
print("-" * 100)
print("dimensions:", node_resp.dims)
print("-" * 100)
print("coordinates:", node_resp.coords)
OPSTOOL :: Loading all response data from _OPSTOOL_ODB/RespStepData-sensitivity.nc ...
data variables: Data variables:
disp (time, nodeTags, DOFs) float64 144kB 0.0 0.0 ... -0.8291
vel (time, nodeTags, DOFs) float64 144kB 0.0 0.0 ... 0.0 0.0
accel (time, nodeTags, DOFs) float64 144kB 0.0 0.0 ... 0.0 0.0
reaction (time, nodeTags, DOFs) float64 144kB 0.0 ... -3.176e-10
reactionIncInertia (time, nodeTags, DOFs) float64 144kB 0.0 ... -3.176e-10
rayleighForces (time, nodeTags, DOFs) float64 144kB 0.0 0.0 ... 0.0 0.0
pressure (time, nodeTags) float64 24kB 0.0 0.0 0.0 ... 0.0 0.0
----------------------------------------------------------------------------------------------------
dimensions: FrozenMappingWarningOnValuesAccess({'time': 501, 'nodeTags': 6, 'DOFs': 6})
----------------------------------------------------------------------------------------------------
coordinates: Coordinates:
* nodeTags (nodeTags) int32 24B 1 2 3 4 5 6
* DOFs (DOFs) <U2 48B 'UX' 'UY' 'UZ' 'RX' 'RY' 'RZ'
* time (time) float64 4kB 0.0 0.002 0.004 0.006 ... 0.994 0.996 0.998 1.0
We use the .sel method to retrieve the displacement of node 2 on UY and total reaction force:
[10]:
times = node_resp["time"].data
disp = node_resp["disp"].sel(nodeTags=ctrlNode, DOFs=dof)
force = -node_resp["reaction"].sel(DOFs=dof).sum(dim="nodeTags")
[11]:
# Plotting
fig, axs = plt.subplots(1, 2, figsize=(10, 3))
axs[0].plot(times, force, lw=2, color="blue")
axs[0].set_xlabel("Time (s)")
axs[1].plot(disp, force, lw=2, color="blue")
axs[1].set_xlabel("Displacement of Node 2 (m)")
for ax in axs:
ax.set_ylabel("Total Reaction Force (kN)")
# ax.grid()
plt.subplots_adjust(wspace=0.3)
plt.show()
Loading the sensitivity response data¶
[12]:
sens_resp = opst.post.get_sensitivity_responses(odb_tag="sensitivity")
print("data variables:", sens_resp.data_vars)
print("-" * 100)
print("dimensions:", sens_resp.dims)
print("-" * 100)
print("coordinates:", sens_resp.coords)
OPSTOOL :: Loading response data from _OPSTOOL_ODB/RespStepData-sensitivity.nc ...
data variables: Data variables:
disp (time, paraTags, nodeTags, DOFs) float64 866kB 0.0 ... 0.0007176
vel (time, paraTags, nodeTags, DOFs) float64 866kB 0.0 0.0 ... 0.0 0.0
accel (time, paraTags, nodeTags, DOFs) float64 866kB 0.0 0.0 ... 0.0 0.0
pressure (time, paraTags, nodeTags) float64 144kB 0.0 0.0 0.0 ... 0.0 0.0
lambda (time, paraTags, patternTags) float64 24kB 0.0 0.0 0.0 ... 0.0 0.0
----------------------------------------------------------------------------------------------------
dimensions: FrozenMappingWarningOnValuesAccess({'time': 501, 'paraTags': 6, 'nodeTags': 6, 'DOFs': 6, 'patternTags': 1})
----------------------------------------------------------------------------------------------------
coordinates: Coordinates:
* paraTags (paraTags) int32 24B 1 2 3 4 5 6
* nodeTags (nodeTags) int32 24B 1 2 3 4 5 6
* DOFs (DOFs) <U2 48B 'UX' 'UY' 'UZ' 'RX' 'RY' 'RZ'
* patternTags (patternTags) int32 4B 1
* time (time) float64 4kB 0.0 0.002 0.004 0.006 ... 0.996 0.998 1.0
Get all parameter tags:
[13]:
paraTags = sens_resp.paraTags.data
print("parameter tags:", paraTags)
parameter tags: [1 2 3 4 5 6]
Sensitivity of node 2 displacement to various parameters¶
Get the displacement sensitivity of the control node on the specified degree of freedom, which stores the sensitivity of each parameter.
paraTags refers to each parameter.
[14]:
sens_disp = sens_resp["disp"].sel(nodeTags=ctrlNode, DOFs=dof)
print(sens_disp)
<xarray.DataArray 'disp' (time: 501, paraTags: 6)> Size: 24kB
array([[ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
[-2.15417342e-11, 0.00000000e+00, 0.00000000e+00,
-2.58500809e-02, 2.58500937e-03, 2.00500802e-04],
[-4.30797758e-11, 0.00000000e+00, 0.00000000e+00,
-5.16957304e-02, 5.16958326e-03, 2.00499448e-04],
...,
[-6.30950787e-10, 4.00315434e-07, 2.93504888e-07,
7.06634167e+00, -1.09069129e+00, 1.89185144e-03],
[-6.06092026e-10, -7.33959131e-08, 3.09006005e-07,
6.70169179e+00, -1.07346954e+00, 2.06308002e-03],
[-5.55271745e-10, -6.10191329e-07, 3.25417143e-07,
6.28457886e+00, -1.05278005e+00, 2.26779623e-03]])
Coordinates:
* paraTags (paraTags) int32 24B 1 2 3 4 5 6
nodeTags int32 4B 2
DOFs <U2 8B 'UY'
* time (time) float64 4kB 0.0 0.002 0.004 0.006 ... 0.994 0.996 0.998 1.0
Plotting¶
The left column shows the sensitivity of displacement \(U\) to each parameter \(P\) and the product of the parameter value:
\[\left( \partial U/\partial {P} \right){P}\]
The right column shows the hysteresis diagram of displacement \(U\) and total reaction force \(F\), where the displacement includes a sensitivity change of 0.1 times, i.e., \(0.1 \left( \partial U/\partial {P} \right){P}\)
[ ]:
fig, axs = plt.subplots(6, 2, figsize=(5 * 2, 2 * 7))
for j, para_tag in enumerate(paraTags):
para_value = ParamVars[para_tag]
para_sym = ParamSym[para_tag]
sens = sens_disp.sel(paraTags=para_tag)
# Subplot #i
# ----------
axs[j, 0].plot(times, sens * para_value, c="#136ad5", linewidth=1.5, label="DDM")
axs[j, 0].fill_between(times, sens * para_value, 0.0, color="#fb8a2e", alpha=0.35)
axs[j, 0].set_ylabel(f"$(\\partial U/\\partial {para_sym}){para_sym}$ [m]")
axs[j, 0].legend(fontsize=9)
if j == 5:
axs[j, 0].set_xlabel(r"Time, $t$")
# Subplot #ii
# -----------
axs[j, 1].plot(disp, force, c="#136ad5", linewidth=2.0, label="$U$")
axs[j, 1].plot(
disp + sens * 0.1 * para_value,
force,
c="#136ad5",
ls="--",
linewidth=1.5,
label=f"$U + (\\partial U/\\partial {para_sym})0.1{para_sym}$",
)
axs[j, 1].plot(
disp - sens * 0.1 * para_value,
force,
c="#136ad5",
ls="-.",
linewidth=1.5,
label=f"$U - (\\partial U/\\partial {para_sym})0.1{para_sym}$",
)
axs[j, 1].fill_betweenx(
force, disp + sens * 0.1 * para_value, disp - sens * 0.1 * para_value, color="#fb8a2e", alpha=0.35
)
axs[j, 1].set_ylabel(r"Load, $P$ [kN]")
axs[j, 1].legend(fontsize=9)
if j == 5:
axs[j, 1].set_xlabel(r"Displacement, $U$ [m]")
for ax in axs.flat:
ax.tick_params(direction="in", length=5, colors="k", width=0.75)
ax.grid(True, color="silver", linestyle="solid", linewidth=0.75, alpha=0.75)
plt.subplots_adjust(wspace=0.3)
plt.show()
