G2 Example 5

G2Example5.py

# ------------------------------------------------------------------------
# The following Python code is implemented by Professor Terje Haukaas at
# the University of British Columbia in Vancouver, Canada. It is made
# freely available online at terje.civil.ubc.ca together with notes,
# examples, and additional Python code. Please be cautious when using
# this code; it may contain bugs and comes without warranty of any kind.
# ------------------------------------------------------------------------

from G2AnalysisNonlinearStatic import *
from G2Model import *

# Input [N, mm]
elementType = 2            # Nonlinear truss element with path-dependent material
materialType = 'BoucWen'   # 'Bilinear' / 'Plasticity' / 'BoucWen'
L = 1000.0                 # Length of element
A = 30*30                  # Cross-section area
F = 400e3                  # Maximum load
N0 = 0.0                   # Initial axial force
E = 200e3                  # Modulus of elasticity
fy = 350                   # Yield stress
alpha = 0.02               # Second-slope stiffness
eta = 3                    # Bouc-Wen sharpness
gamma = 0.5                # Bouc-Wen parameter
beta = 0.5                 # Bouc-Wen parameter
H = alpha*E/(1-alpha)      # Kinematic hardening parameter
K = 0                      # Linear isotropic hardening parameter
delta = 0                  # Saturation isotropic hardening parameter
fy_inf = 100.0             # Asymptotic yield stress for saturation isotropic hardening
dt = 0.2                   # Delta-t
nsteps = 20                # Number of pseudo-time steps, each of length dt
KcalcFrequency = 1         # 0=initial stress method, 1=Newton-Raphson, maxIter=Modified NR
maxIter = 100              # Maximum number of equilibrium iterations in the Newton-Raphson algorithm
tol = 1e-5                 # Convergence tolerance for the Newton-Raphson algorithm
trackNode = 2              # Node to be plotted
trackDOF = 1               # DOF to be plotted

# Nodal coordinates
NODES = [[0.0, 0.0],
         [L,   0.0]]

# Boundary conditions (0=free, 1=fixed, sets #DOFs per node)
CONSTRAINTS = [[1, 1],
               [0, 1]]

# Element
ELEMENTS = [[elementType, N0, 1, 2]]

# Section
SECTIONS = [['Truss', A]]

# Material
if materialType == 'Bilinear':
    MATERIALS = [['Bilinear', E, fy, alpha]]
elif materialType == 'Plasticity':
    MATERIALS = [['Plasticity', E, fy, H, K, delta, fy_inf]]
elif materialType == 'BoucWen':
    MATERIALS = [['BoucWen', E, fy, alpha, eta, beta, gamma]]
else:
    print('\n'"Cannot understand the material type")
    import sys
    sys.exit()

# Nodal loads
LOADS = [[0.0, 0.0],
         [F,   0.0]]

# Lumped mass
MASS = [[0, 0],
        [0, 0]]

# Create the model object
a = [NODES, CONSTRAINTS, ELEMENTS, SECTIONS, MATERIALS, LOADS, MASS]
m = model(a)

# Do a DDM check
selectedDDMparameter = 'E'
perturbationFraction = 1e-6
theValue = 'value'
exec("%s = %s" % (theValue, selectedDDMparameter))
DDMparameters = [['Element', selectedDDMparameter, range(1, 2)]]
loadFactor, responseOrig, ddmSensitivity = nonlinearStaticAnalysis(m, nsteps, dt, maxIter, KcalcFrequency, tol, trackNode, trackDOF, DDMparameters)
m = model(a)
m.setParameters([['Element', selectedDDMparameter, range(1, 2), value*(1.0+perturbationFraction)]])
loadFactor, responsePert, blank = nonlinearStaticAnalysis(m, nsteps, dt, maxIter, KcalcFrequency, tol, trackNode, trackDOF, [])
fdmSensitivity = 1.0/(value*perturbationFraction) * np.subtract(responsePert, responseOrig)
print('\n'"FDM-DDM difference:", np.max(np.abs(np.subtract(fdmSensitivity, ddmSensitivity[:, 0]))))

# Plot the DDM check
plt.ion()
plt.figure()
plt.autoscale(True)
plt.title("Response Sensitivity for $\\theta$ = " + selectedDDMparameter)
t = np.linspace(1, len(ddmSensitivity), len(ddmSensitivity))
plt.plot(t, fdmSensitivity, 'ro-', label='FDM')
plt.plot(t, ddmSensitivity[:, 0], 'ko-', label='DDM', markersize=3)
plt.xlabel("Load Increment")
plt.ylabel("$\partial u \; / \; \partial \\theta$")
plt.legend(loc='lower right')
print('\n'"Click somewhere in the plot to continue...")
plt.waitforbuttonpress()