diff --git a/source/user/manual/material/uniaxialMaterial.rst b/source/user/manual/material/uniaxialMaterial.rst index 393f955d..ab7679d7 100644 --- a/source/user/manual/material/uniaxialMaterial.rst +++ b/source/user/manual/material/uniaxialMaterial.rst @@ -27,6 +27,7 @@ The following subsections contain information about **$matType** uniaxialMaterials/Steel01 uniaxialMaterials/Steel02 uniaxialMaterials/Steel4 + uniaxialMaterials/ASDSteel1D uniaxialMaterials/ReinforcingSteel uniaxialMaterials/DoddRestrepo uniaxialMaterials/RambergOsgoodSteel diff --git a/source/user/manual/material/uniaxialMaterials/ASDSteel1D.rst b/source/user/manual/material/uniaxialMaterials/ASDSteel1D.rst new file mode 100644 index 00000000..f0322e53 --- /dev/null +++ b/source/user/manual/material/uniaxialMaterials/ASDSteel1D.rst @@ -0,0 +1,173 @@ +.. _ASDSteel1D: + +.. |integer| replace:: ``integer`` +.. |string| replace:: ``string`` +.. |float| replace:: ``float`` + +ASDSteel1D Material +^^^^^^^^^^^^^^^^^^^ + +.. image:: figures/ASDSteel1D/ASDSteel1D_Ex_numerical_vs_experimental.gif + :align: center + +| The ASDSteel1D command is used to create a uniaxial material object that models the nonlinear response of reinforcing steel bars under monotonic and cyclic loading. +| The formulation integrates kinematic hardening plasticity, damage mechanics, buckling behavior and bond-slip mechanism. +| To improve numerical robustness and computational performance under severe nonlinearities an optional IMPL-EX integration scheme is available. +| Additionally, the entire model is regularized using the element's characteristic length to guarantee mesh-independent results. + +.. function:: + uniaxialMaterial ASDSteel1D $tag $E $sy $su $eu + <-implex> + <-auto_regularization> + <-buckling $lch <$r>> + <-fracture <$r>> + <-slip $matTag <$r>> + <-K_alpha $K_alpha> + <-max_iter $max_iter> + <-tolU $tolU> + <-tolR $tolR> + +.. csv-table:: + :header: "Argument", "Type", "Description" + :widths: 10, 10, 40 + + $tag, |integer|, "Unique tag identifying this material." + $E $sy $su $eu, 4 |float|, "Mandatory. Young's modulus, Yield stress, Ultimate stress and Ultimate strain." + -implex, |string|, "Optional. If defined, the IMPL-EX integration will be used, otherwise the standard implicit integration will be used (default)." + -auto_regularization, |string|, "Optional. Activates automatic regularization based on the characteristic length of the finite element." + -buckling $lch <$r>, |string| + |float| <+ |float|>, "Optional. Enables buckling simulation using an RVE-based approach. Requires characteristic length $lch and optionally a section radius $r." + -fracture <$r>, |string| <+ |float|>, "Optional. Activates fracture modeling. Optionally specify the section radius $r." + -slip $matTag <$r>, |string| + |integer| <+ |float|>, "Optional. Activates slip modeling with a secondary uniaxial material ($matTag). Optionally specify the section radius $r." + -K_alpha $K_alpha, |string| + |float|, "Optional. Defines the weight between the consistent elastoplastic tangent modulus and the purely elastic modulus (default = 0.5). Set to 1.0 for full consistent tangent, or 0.0 to use only the elastic modulus." + -max_iter $max_iter, |string| + |float|, "Optional. Maximum number of iterations for the global Newton-Raphson loop used in the RVE (default = 100)." + -tolU $tolU, |string| + |float|, "Optional. Tolerance on displacement increment convergence (default = 1e-6)." + -tolR $tolR, |string| + |float|, "Optional. Tolerance on residual force convergence (default = 1e-6)." + +Theory +"""""" + +| The ASDSteel1D material is formulated within an elasto-plastic framework enriched with damage, buckling, and slip mechanisms, specifically tailored to model the nonlinear behavior of reinforcing steel bars. + +Plasticity +"""""""""" + +| The plasticity model employs a two-term Chaboche kinematic hardening law. + +| The trial stress is computed as: + +.. math:: + \sigma_{trial} = \sigma_n + E \Delta \varepsilon + +| The yield function is defined as: + +.. math:: + f = | \sigma - \alpha_1 - \alpha_2 | - \sigma_y + +| The evolution of the backstresses follows: + +.. math:: + \dot{\alpha}_i = H_i n - \gamma_i \alpha_i, \quad i=1,2 + +.. math:: + n = \text{sign}(\sigma - \alpha_1 - \alpha_2) + +| The parameters :math:`H_1, \gamma_1` and :math:`H_2, \gamma_2` are calibrated to fit the target values of yield stress (:math:`\sigma_y`), ultimate stress (:math:`\sigma_u`), and ultimate strain (:math:`\varepsilon_u`). + + +Damage +"""""" + +| Fracture is controlled by a scalar damage variable :math:`d \in [0, 1]`, which progressively reduces the effective stress: + +.. math:: + \sigma = (1-d) \sigma_{eff} + +| Damage initiates once the equivalent plastic strain :math:`\varepsilon_{pl}` ​exceeds the ultimate plastic strain :math:`\varepsilon_{upl}`, and evolves as: + +.. math:: + d = 1 - \frac{\sigma_d}{\sigma_y}, \quad \sigma_d = \sigma_y \exp \left( -\frac{\sigma_y (\varepsilon_{pl} - \varepsilon_{upl})}{G_f} \right) + +Buckling +"""""""" + +| Buckling is modeled using a homogenization-inspired multiscale approach, in which the macroscopic behavior of the material is enriched by the nonlinear response of a localized Representative Volume Element (RVE). +| The RVE is defined at a lower scale and incorporates the minimal microstructural features required to represent local instability phenomena without modeling them explicitly across the entire structure. +| At this scale, we employ a generic configuration of corotational beam elements with geometric imperfections, capable of autonomously developing buckling patterns under increasing compressive loads. +| The multiscale coupling is carried out through a two-way exchange of information between the macroscopic and microscopic scales. +| In the downscaling step, the strain computed at the macroscopic level is transferred to the RVE and applied as a boundary condition. +| Once the RVE problem is solved, the upscaling step follows. The resulting stress field, stiffness degradation, or other internal variables obtained from the RVE are fed back into the macroscopic model. +| This bidirectional process ensures that the influence of microstructural effects is consistently reflected at the structural scale. + +.. image:: figures/ASDSteel1D/disegnohomog.png + :align: center + +Bond-Slip +""""""""" + +| A bond–slip material can be inserted in series with the steel to simulate anchor slip phenomena. The following conditions apply: + +.. math:: + \sigma_{\text{steel}} = \sigma_{\text{slip}} + +.. math:: + \varepsilon_{\text{total}} = \varepsilon_{\text{steel}} + \varepsilon_{\text{slip}} + +| The slip law (e.g., :math:`\tau–slip` curve) is defined by the user via an external uniaxial material and internally converted into a stress–strain relationship. + +Regularization +"""""""""""""" + +| Three of the four nonlinear mechanisms included in the model, damage, buckling, and bond-slip, can exhibit strain-softening behavior, which leads to mesh sensitivity when modeled without proper regularization. +| To ensure that the material response remains objective and independent of the element size, each of these mechanisms is regularized through tailored strategies based on the characteristic length of the finite element. + +| For the bond–slip mechanism, regularization must be handled by the secondary material used in series. At present, the only uniaxial material that supports this form of regularization is the :ref:`ASDConcrete1D`. + +IMPL-EX scheme +"""""""""""""" + +| An optional IMPL-EX integration scheme can be activated to improve convergence behavior. +| This scheme combines the efficiency of explicit integration with the stability of implicit correction, offering improved robustness in simulations involving highly nonlinear material behavior, such as softening or large strain evolution. +| At each time step, the algorithm first executes an explicit predictor step, estimating the evolution of internal variables without solving a full nonlinear system. +| This is followed by an implicit corrector step, which refines the predicted state to ensure consistency and numerical stability. + +.. math:: + \Delta \lambda_n \approx \frac{\Delta t_n}{\Delta t_{n-1}} (\lambda_n - \lambda_{n-1}) + + +Usage Notes +""""""""""" + +.. admonition:: Responses + + * All standard responses for a uniaxialMaterial object: **stress**, **strain**, **tangent**. + * **BucklingIndicator** or **BI**: The lateral displacement caused by buckling (available with the **-buckling** option). + * **EquivalentPlasticStrain** or **PLE**: The equivalent plastic strain. + * **Damage** or **D**: Current damage variable (available with the **-fracture** option). + * **SlipResponse**: 2 components (:math:`slip`, :math:`\tau`). Current bond-slip response (available with the **-slip** option). + +.. admonition:: Example 1 - Fracture + + A Python example to simulate tensile fracture: :download:`ASDSteel1D_Ex_Damage.py ` + + .. image:: figures/ASDSteel1D/ASDSteel1D_Ex_Damage_Output.gif + +.. admonition:: Example 2 - Buckling + + A Python example to simulate buckling: :download:`ASDSteel1D_Ex_Buckling.py ` + + .. image:: figures/ASDSteel1D/ASDSteel1D_Ex_Buckling_Output.gif + +.. admonition:: Example 3 - Bond-Slip + + A Python example to simulate rebars bond-slip: :download:`ASDSteel1D_Ex_Bond_Slip.py ` + + .. image:: figures/ASDSteel1D/ASDSteel1D_Ex_Bond_Slip_Output.gif + +Code Developed by: **Alessia Casalucci** at ASDEA Software, Italy. + +References +"""""""""" + + +.. [Kashani2015] Kashani, M. M., Lowes, L. N., Crewe, A. J., & Alexander, N. A. (2015). "Phenomenological hysteretic model for corroded reinforcing bars including inelastic buckling and low-cycle fatigue degradation" Computers & Structures, 156, 58-71 (`Link to article `_) +.. [Oliver2008] Oliver, J., Huespe, A. E., & Cante, J. C. (2008). "An implicit/explicit integration scheme to increase computability of non-linear material and contact/friction problems" Computer Methods in Applied Mechanics and Engineering, 197(21-24), 1865-1889 (`Link to article `_) diff --git a/source/user/manual/material/uniaxialMaterials/examples/ASDSteel1D/ASDSteel1D_Ex_Bond_Slip.py b/source/user/manual/material/uniaxialMaterials/examples/ASDSteel1D/ASDSteel1D_Ex_Bond_Slip.py new file mode 100644 index 00000000..f3cca9e3 --- /dev/null +++ b/source/user/manual/material/uniaxialMaterials/examples/ASDSteel1D/ASDSteel1D_Ex_Bond_Slip.py @@ -0,0 +1,279 @@ +import numpy as np +import math +import random +import matplotlib.pyplot as plt +from matplotlib.collections import LineCollection +from scipy.stats import norm +import time +import openseespy.opensees as ops +from PIL import Image + + +#Function for single run +def run_simulation(tim_noise, mod_noise, title, Npulse, Nsubsteps=50, do_plot=True): + # USER SETTINGS + # Steel + Es = 210000.0 + fy = 400.0 + fu = 520.0 + eu = 0.06 + # Rebars bond-slip + fc = 25.0 + good_bond = True + radius = 8.0 + # finite element length (just to test mesh-dependency) + lch = 200.0 + # Integration + implex = True + auto_regularization = True + bond_slip = True + # Compute bond-slip + alpha = 0.4 + s3 = max(25.4, 2.0*radius) + if good_bond: + beta = 2.5 + s1 = 1.0 + s2 = 2.0 + else: + beta = 1.25 + s1 = 1.8 + s2 = 3.6 + tau_max = beta * math.sqrt(fc) + tau_f = 0.4 * tau_max + # compute slip values + Te = np.concatenate((np.linspace(0.0, s1, 20), np.array([s2, s3]))) + # compute tau values + def bond_fun(s): + if s <= s1: + return tau_max * (s / s1) ** alpha + elif s <= s2: + return tau_max + elif s <= s3: + return tau_max - (tau_max - tau_f) * (s - s2) / (s3 - s2) + else: + return tau_f + + Ts = [bond_fun(slip) for slip in Te] + + # damage + Td = [0.2] * len(Te) + # Initial stiffness of the bond-slip material + bond_stiffness = Ts[1] / Te[1] + # Just to define an almost zero plateau + zero_tau = 1.0e-3 + zero_tau_slip = zero_tau / bond_stiffness + Ne = [0.0, zero_tau_slip, zero_tau_slip * 2.0] + Ns = [0.0, zero_tau, zero_tau ] + Nd = [1.0, 1.0, 1.0] + + # the 1D model + ops.wipe() + ops.model('basic', '-ndm', 1, '-ndf', 1) + + # the bond material + ops.uniaxialMaterial('ASDConcrete1D', 21, bond_stiffness, '-Te', *Te, '-Ts', *Ts, '-Td', *Td, + '-Ce', *Ne, '-Cs', *Ns, '-Cd', *Nd, '-secant', '-implex', '-autoRegularization',1.0) + ops.uniaxialMaterial('ASDConcrete1D', 22, bond_stiffness, '-Te', *Ne, '-Ts', *Ns, '-Td', *Nd, + '-Ce', *Te, '-Cs', *Ts, '-Cd', *Td, '-secant', '-implex', '-autoRegularization', 1.0) + ops.uniaxialMaterial('ASDConcrete1D', 23, bond_stiffness, '-Te', *Te, '-Ts', *Ts, '-Td', *Td, + '-Ce', *Te, '-Cs', *Ts, '-Cd', *Td, '-secant', '-implex', '-autoRegularization', 1.0) + + stiff_c = bond_stiffness*100.0 + Te = [0.0, zero_tau/stiff_c, zero_tau/stiff_c * 2.0] + Ts = [0.0, zero_tau, zero_tau ] + Td = [1.0, 1.0, 1.0] + Ce = [0.0, 1.0/stiff_c, 2.0/stiff_c] + Cs = [0.0, 1.0, 2.0 ] + ops.uniaxialMaterial('ASDConcrete1D', 41, stiff_c, '-Te', *Te, '-Ts', *Ts, '-Td', *Td, + '-Ce', *Ce, '-Cs', *Cs, '-Cd', *Td, '-secant', '-implex') + + ops.uniaxialMaterial('Parallel', 2, 21, 22, 23, 41, '-factors', 0.5, 0.5, 0.5, 1.0) + + # the steel material + # steel mandatory parameters + steel_params = [Es, fy, fu, eu] + if implex: + steel_params.append('-implex') + if auto_regularization: + steel_params.append('-auto_regularization') + if bond_slip: + steel_params.extend(['-slip', 2, radius]) + ops.uniaxialMaterial('ASDSteel1D', 1, *steel_params) + + # a truss + ops.node(1, 0) + ops.node(2, lch) + ops.element('truss', 1, 1, 2, 1.0, 1) + # fixity + ops.fix(1, 1) + + # a simple ramp + ops.timeSeries('Path', 1, '-time', *tim_noise , '-values', *mod_noise ) + ops.pattern('Plain', 1, 1) + ops.sp(2, 1, lch) + + # some default analysis settings + ops.constraints('Transformation') + ops.numberer('Plain') + ops.system('UmfPack') + ops.test('NormDispIncr', 1.0e-8, 10, 0) + ops.algorithm('Newton') + time_incr = 1.0 / (Npulse * Nsubsteps) + ops.integrator('LoadControl', time_incr) + ops.analysis('Static') + + + # begin plot + STRAIN = [0.0] + STRESS = [0.0] + SLIP = [0.0] + BOND = [0.0] + TIME = [0.0] + STRAIN_TIP = [0.0] + STRESS_TIP = [0.0] + TIME_TIP = [0.0] + + # plot settings + SMALL_SIZE = 8 + MEDIUM_SIZE = 10 + LARGE_SIZE = 12 + plt.rc('font', size=SMALL_SIZE) # controls default text sizes + plt.rc('axes', titlesize=MEDIUM_SIZE) # fontsize of the axes title + plt.rc('figure', titlesize=LARGE_SIZE) # fontsize of the figure title + plt.ion() + + # figure and axes + fig = plt.figure(figsize=(12, 6), dpi=100) + spec = plt.GridSpec(nrows=4, ncols=6,figure=fig) + ax = [fig.add_subplot(spec[0,:]), fig.add_subplot(spec[1:,0:3]), fig.add_subplot(spec[1:,3:])] + + ax[0].grid(linestyle=':') + ax[0].set_title('Strain history') + ax[0].set_xlabel('Time') + ax[0].set_ylabel('\N{GREEK SMALL LETTER EPSILON}') + + ax[1].grid(linestyle=':') + ax[1].set_title('Stress-strain response') + ax[1].set_xlabel('\N{GREEK SMALL LETTER EPSILON}') + ax[1].set_ylabel('\N{GREEK SMALL LETTER SIGMA}') + + ax[2].grid(linestyle=':') + ax[2].set_title('Bond-Slip') + ax[2].set_xlabel('Slip') #Tau vs slip + ax[2].set_ylabel('\N{GREEK SMALL LETTER TAU}') + + # the strain history + the_line_hist_bg, = ax[0].plot(tim_noise, mod_noise, '--k', linewidth=1.0) + the_line_hist, = ax[0].plot(TIME, STRAIN, 'k', linewidth=1.0) + the_tip_hist, = ax[0].plot(TIME_TIP, STRAIN_TIP, 'or', fillstyle='full', markersize=8) + + # the composite response + the_line, = ax[1].plot(STRAIN, STRESS, '-k', linewidth=1.0) + the_tip, = ax[1].plot(STRAIN_TIP, STRESS_TIP, 'or', fillstyle='full', markersize=8) + + # the slip response + the_line_slip, = ax[2].plot(SLIP, BOND, '-r', linewidth=1.0) + + + # done + fig.tight_layout(pad=2.0) + + # process each step + time_start = time.time() + frames = [] + for i in range(Npulse): + failed = False + for j in range(Nsubsteps): + ok = ops.analyze(1) + if not do_plot: + if ok != 0: + failed = True + break + continue + if ok == 0: + # time + TIME.append(ops.getTime()) + # composite response + STRAIN.append(ops.eleResponse(1, 'material', 1, 'strain')[0]) + STRESS.append(ops.eleResponse(1, 'material', 1, 'stress')[0]) + # bond-slip response + try: + bs = ops.eleResponse(1, 'material', 1, 'SlipResponse') + SLIP.append(bs[0]) + BOND.append(bs[1]) + except: + BOND.append(0.0) + + # update + the_line_hist.set_xdata(TIME) + the_line_hist.set_ydata(STRAIN) + the_line.set_xdata(STRAIN) + the_line.set_ydata(STRESS) + the_line_slip.set_xdata(SLIP) + the_line_slip.set_ydata(BOND) + + else: + failed = True + break + if failed: + break + STRAIN_TIP = [STRAIN[-1]] + STRESS_TIP = [STRESS[-1]] + TIME_TIP = [TIME[-1]] + if i == Npulse-1: + the_tip.remove() + the_tip_hist.remove() + else: + the_tip.set_xdata(STRAIN_TIP) + the_tip.set_ydata(STRESS_TIP) + the_tip_hist.set_xdata(TIME_TIP) + the_tip_hist.set_ydata(STRAIN_TIP) + for iax in ax: + iax.relim() + iax.autoscale_view() + + fig.canvas.draw() + fig.canvas.flush_events() + plt.pause(0.001) + # store frame + w, h = fig.canvas.get_width_height() + buf = np.frombuffer(fig.canvas.buffer_rgba(), dtype=np.uint8) + buf.shape = (h, w, 4) + frame = buf[:, :, :3].copy() + + frames.append(frame) + time_end = time.time() + + plt.ioff() + plt.close(fig) + ops.wipe() + + print(f'Saving GIF animation {title}.gif...') + images = [Image.fromarray(frame) for frame in frames] + images[0].save('ASDSteel1D_Ex_Bond_Slip_Output.gif', format='GIF', save_all=True, append_images=images[1:], loop=0, duration=100, disposal=2) + + +# To generate load history +Npulse = 50 +seed=2 +max_scale=4.0 +eu = 0.06 +random.seed(seed) +np.random.seed(seed) + +mod_noise = np.linspace(0, eu * max_scale, Npulse) +mod_noise = list(mod_noise) + +tim_noise = [0.0] * Npulse +for i in range(1, Npulse): + dy = abs(mod_noise[i] - mod_noise[i-1]) + tim_noise[i] = tim_noise[i-1] + dy + +tim_noise_max = tim_noise[-1] +if tim_noise_max == 0.0: + tim_noise_max = 1.0 +tim_noise = [t / tim_noise_max for t in tim_noise] + + +title = f'results__slip_run' +run_simulation(tim_noise, mod_noise, title, Npulse=Npulse) \ No newline at end of file diff --git a/source/user/manual/material/uniaxialMaterials/examples/ASDSteel1D/ASDSteel1D_Ex_Buckling.py b/source/user/manual/material/uniaxialMaterials/examples/ASDSteel1D/ASDSteel1D_Ex_Buckling.py new file mode 100644 index 00000000..57592885 --- /dev/null +++ b/source/user/manual/material/uniaxialMaterials/examples/ASDSteel1D/ASDSteel1D_Ex_Buckling.py @@ -0,0 +1,259 @@ +import numpy as np +import math +import random +import matplotlib.pyplot as plt +from matplotlib.collections import LineCollection +from scipy.stats import norm +import time +import openseespy.opensees as ops +from PIL import Image + +#Function for single run +def run_simulation(tim_noise, mod_noise, title, Npulse=100, Nsubsteps=300, do_plot=True): + # USER SETTINGS + # Steel + Es = 210000.0 + fy = 400.0 + fu = 520.0 + eu = 0.06 + # Rebars buckling + radius = 8.0 + lch_buckling = 200.0 + # finite element length (just to test mesh-dependency) + lch = 200.0 + # Integration + implex = True + auto_regularization = True + buckling = True + + # the 1D model + ops.wipe() + ops.model('basic', '-ndm', 1, '-ndf', 1) + + # the steel material + # steel mandatory parameters + steel_params = [Es, fy, fu, eu] + if implex: + steel_params.append('-implex') + if auto_regularization: + steel_params.append('-auto_regularization') + if buckling: + steel_params.extend(['-buckling', lch_buckling, radius]) + ops.uniaxialMaterial('ASDSteel1D', 1, *steel_params) + + # a truss + ops.node(1, 0) + ops.node(2, lch) + ops.element('truss', 1, 1, 2, 1.0, 1) + # fixity + ops.fix(1, 1) + + # a simple ramp + ops.timeSeries('Path', 1, '-time', *tim_noise , '-values', *mod_noise ) + ops.pattern('Plain', 1, 1) + ops.sp(2, 1, lch) #ops.load(2,0,1) + + # some default analysis settings + ops.constraints('Transformation') + ops.numberer('Plain') + ops.system('UmfPack') + ops.test('NormDispIncr', 1.0e-8, 10, 0) + ops.algorithm('Newton') + time_incr = 1.0 / (Npulse * Nsubsteps) + ops.integrator('LoadControl', time_incr) + ops.analysis('Static') + + # begin plot + STRAIN = [0.0] + STRESS = [0.0] + RVE = [0.0] + TIME = [0.0] + STRAIN_TIP = [0.0] + STRESS_TIP = [0.0] + TIME_TIP = [0.0] + + # plot settings + SMALL_SIZE = 8 + MEDIUM_SIZE = 10 + LARGE_SIZE = 12 + plt.rc('font', size=SMALL_SIZE) # controls default text sizes + plt.rc('axes', titlesize=MEDIUM_SIZE) # fontsize of the axes title + plt.rc('figure', titlesize=LARGE_SIZE) # fontsize of the figure title + plt.ion() + + # figure and axes + fig = plt.figure(figsize=(12, 6), dpi=100) + spec = plt.GridSpec(nrows=5, ncols=5,figure=fig) + ax = [fig.add_subplot(spec[0,:]), fig.add_subplot(spec[1:,0:3]), fig.add_subplot(spec[1:,3:])] + + ax[0].grid(linestyle=':') + ax[0].set_title('Strain history') + ax[0].set_xlabel('Time') + ax[0].set_ylabel('\N{GREEK SMALL LETTER EPSILON}') + + ax[1].grid(linestyle=':') + ax[1].set_title('Stress-strain response') + ax[1].set_xlabel('\N{GREEK SMALL LETTER EPSILON}') + ax[1].set_ylabel('\N{GREEK SMALL LETTER SIGMA}') + + strain = STRAIN[-1] + l_eff = lch * (1 + strain) + ax[2].grid(linestyle=':') + ax[2].set_title('RVE') + ax[2].set_xlabel('Lateral Displacement') #Lateral vs vertical displacement + ax[2].set_ylabel('Vertical Displacement') + ax[2].axvline(x=0, color='gray', linestyle = '--', linewidth=0.8) + ax[2].set_xlim(-lch,lch) + ax[2].set_ylim(-lch/20, l_eff + lch/20) + + # the strain history + the_line_hist_bg, = ax[0].plot(tim_noise, mod_noise, '--k', linewidth=1.0) + the_line_hist, = ax[0].plot(TIME, STRAIN, 'k', linewidth=1.0) + the_tip_hist, = ax[0].plot(TIME_TIP, STRAIN_TIP, 'or', fillstyle='full', markersize=8) + + # the composite response + the_line, = ax[1].plot(STRAIN, STRESS, '-k', linewidth=1.0) + the_tip, = ax[1].plot(STRAIN_TIP, STRESS_TIP, 'or', fillstyle='full', markersize=8) + + #RVE + y_shape = np.linspace(0,lch,100) + x_shape = np.zeros_like(y_shape) + points = np.array([x_shape, y_shape]).T.reshape(-1, 1, 2) + segments = np.concatenate([points[:-1], points[1:]], axis=1) + + # Initial uniform linewidth + linewidths = np.full(len(segments), radius) + + # Create LineCollection + the_line_rve = LineCollection(segments, colors='k', linewidths=linewidths) + ax[2].add_collection(the_line_rve) + # done + fig.tight_layout(pad=2.0) + + # process each step + time_start = time.time() + frames = [] + for i in range(Npulse): + failed = False + for j in range(Nsubsteps): + ok = ops.analyze(1) + if not do_plot: + if ok != 0: + failed = True + break + continue + if ok == 0: + # time + TIME.append(ops.getTime()) + # composite response + STRAIN.append(ops.eleResponse(1, 'material', 1, 'strain')[0]) + STRESS.append(ops.eleResponse(1, 'material', 1, 'stress')[0]) + # RVE response + try: + RVE.append(ops.eleResponse(1, 'material', 1, 'BucklingIndicator')[0]) + except: + RVE.append(0.0) + # update + the_line_hist.set_xdata(TIME) + the_line_hist.set_ydata(STRAIN) + the_line.set_xdata(STRAIN) + the_line.set_ydata(STRESS) + A = RVE[-1] + strain = STRAIN[-1] + l_eff = lch * (1 + strain) + + # coordinate Y aggiornate + y_half = np.linspace(0, l_eff/2, 100) + + def buckling_cubic(y, l, A): + a = -4 * A / l**3 + b = 6 * A / l**2 + return a * y**3 + b * y**2 + + x_half = buckling_cubic(y_half, l_eff/2, A) + y_shape = np.concatenate([y_half, l_eff - y_half[::-1]]) + x_shape = np.concatenate([x_half, +x_half[::-1]]) + points = np.array([x_shape, y_shape]).T.reshape(-1, 1, 2) + segments = np.concatenate([points[:-1], points[1:]], axis=1) + + # Update the line collection + the_line_rve.set_segments(segments) + + n_segments = len(segments) + mid = n_segments // 2 + + sigma = n_segments / 20 + + # smooth profile (Maximum at the center and gradually decreases toward the sides) + x = np.arange(n_segments) + gaussian_profile = np.exp(-((x - mid) ** 2) / (2 * sigma ** 2)) + gaussian_profile = gaussian_profile / gaussian_profile.max() + linewidths = radius * 2.0 + + the_line_rve.set_linewidths(linewidths) + else: + failed = True + break + if failed: + break + STRAIN_TIP = [STRAIN[-1]] + STRESS_TIP = [STRESS[-1]] + TIME_TIP = [TIME[-1]] + if i == Npulse-1: + the_tip.remove() + the_tip_hist.remove() + else: + the_tip.set_xdata(STRAIN_TIP) + the_tip.set_ydata(STRESS_TIP) + the_tip_hist.set_xdata(TIME_TIP) + the_tip_hist.set_ydata(STRAIN_TIP) + for iax in ax: + iax.relim() + iax.autoscale_view() + + fig.canvas.draw() + fig.canvas.flush_events() + plt.pause(0.001) + # store frame + w, h = fig.canvas.get_width_height() + buf = np.frombuffer(fig.canvas.buffer_rgba(), dtype=np.uint8) + buf.shape = (h, w, 4) + frame = buf[:, :, :3].copy() + + frames.append(frame) + time_end = time.time() + + plt.ioff() + plt.close(fig) + ops.wipe() + + print(f'Saving GIF animation {title}.gif...') + images = [Image.fromarray(frame) for frame in frames] + images[0].save('ASDSteel1D_Ex_Buckling_Output.gif', format='GIF', save_all=True, append_images=images[1:], loop=0, duration=100, disposal=2) + + +# To generate load history +Npulse = 100 +seed=2 +max_scale=4.0 +min_scale=4.0 +eu = 0.06 +random.seed(seed) +np.random.seed(seed) + +mod_noise = np.linspace(0, -eu * min_scale, Npulse) +mod_noise = list(mod_noise) + +tim_noise = [0.0] * Npulse +for i in range(1, Npulse): + dy = abs(mod_noise[i] - mod_noise[i-1]) + tim_noise[i] = tim_noise[i-1] + dy + +tim_noise_max = tim_noise[-1] +if tim_noise_max == 0.0: + tim_noise_max = 1.0 +tim_noise = [t / tim_noise_max for t in tim_noise] + + +title = f'results__buck_run' +run_simulation(tim_noise, mod_noise, title, Npulse=Npulse) diff --git a/source/user/manual/material/uniaxialMaterials/examples/ASDSteel1D/ASDSteel1D_Ex_Damage.py b/source/user/manual/material/uniaxialMaterials/examples/ASDSteel1D/ASDSteel1D_Ex_Damage.py new file mode 100644 index 00000000..ed7a9a37 --- /dev/null +++ b/source/user/manual/material/uniaxialMaterials/examples/ASDSteel1D/ASDSteel1D_Ex_Damage.py @@ -0,0 +1,283 @@ +import numpy as np +import math +import random +import matplotlib.pyplot as plt +from matplotlib.collections import LineCollection +from scipy.stats import norm +import time +import openseespy.opensees as ops +from PIL import Image + +#Function for single run +def run_simulation(tim_noise, mod_noise, max_scale, title, Npulse=100, Nsubsteps=300, do_plot=True): + # USER SETTINGS + # Steel + Es = 210000.0 + fy = 400.0 + fu = 520.0 + eu = 0.06 + # fracture + radius = 8.0 + # finite element length (just to test mesh-dependency) + lch = 100.0 + # Integration + implex = True + auto_regularization = True + fracture = True + + # the 1D model + ops.wipe() + ops.model('basic', '-ndm', 1, '-ndf', 1) + + # the steel material + # steel mandatory parameters + steel_params = [Es, fy, fu, eu] + if implex: + steel_params.append('-implex') + if auto_regularization: + steel_params.append('-auto_regularization') + if fracture: + steel_params.extend(['-fracture', radius]) + ops.uniaxialMaterial('ASDSteel1D', 1, *steel_params) + + # a truss + ops.node(1, 0) + ops.node(2, lch) + ops.element('truss', 1, 1, 2, 1.0, 1) + # fixity + ops.fix(1, 1) + + # a simple ramp + ops.timeSeries('Path', 1, '-time', *tim_noise , '-values', *mod_noise ) + ops.pattern('Plain', 1, 1) + ops.sp(2, 1, lch) + + # some default analysis settings + ops.constraints('Transformation') + ops.numberer('Plain') + ops.system('UmfPack') + ops.test('NormDispIncr', 1.0e-8, 10, 0) + ops.algorithm('Newton') + time_incr = 1.0 / (Npulse * Nsubsteps) + ops.integrator('LoadControl', time_incr) + ops.analysis('Static') + + # begin plot + STRAIN = [0.0] + STRESS = [0.0] + DAMAGE = [0.0] + RVE = [0.0] + EQP = [0.0] + TIME = [0.0] + STRAIN_TIP = [0.0] + STRESS_TIP = [0.0] + TIME_TIP = [0.0] + + # plot settings + SMALL_SIZE = 8 + MEDIUM_SIZE = 10 + LARGE_SIZE = 12 + plt.rc('font', size=SMALL_SIZE) # controls default text sizes + plt.rc('axes', titlesize=MEDIUM_SIZE) # fontsize of the axes title + plt.rc('figure', titlesize=LARGE_SIZE) # fontsize of the figure title + plt.ion() + + # figure and axes + fig = plt.figure(figsize=(12, 6), dpi=100) + spec = plt.GridSpec(nrows=5, ncols=7,figure=fig) + ax = [fig.add_subplot(spec[0,:]), fig.add_subplot(spec[1:,0:3]), fig.add_subplot(spec[1:3,5:]), fig.add_subplot(spec[3:5,5:]), fig.add_subplot(spec[1:,3:5])] + + ax[0].grid(linestyle=':') + ax[0].set_title('Strain history') + ax[0].set_xlabel('Time') + ax[0].set_ylabel('\N{GREEK SMALL LETTER EPSILON}') + + ax[1].grid(linestyle=':') + ax[1].set_title('Stress-strain response') + ax[1].set_xlabel('\N{GREEK SMALL LETTER EPSILON}') + ax[1].set_ylabel('\N{GREEK SMALL LETTER SIGMA}') + + ax[2].grid(linestyle=':') + ax[2].set_title('Damage') + ax[2].set_xlabel('Time') + ax[2].set_ylabel('D') + + ax[3].grid(linestyle=':') + ax[3].set_title('Equivalent plastic strain') + ax[3].set_xlabel('Time') #Equivalent plastic strain vs time + ax[3].set_ylabel('\N{GREEK SMALL LETTER EPSILON}\u209A\u2097') + + strain = STRAIN[-1] + l_eff = lch * (1 + strain) + ax[4].grid(linestyle=':') + ax[4].set_title('RVE') + ax[4].set_xlabel('Lateral Displacement') #Lateral vs vertical displacement + ax[4].set_ylabel('Vertical Displacement') + ax[4].axvline(x=0, color='gray', linestyle = '--', linewidth=0.8) + ax[4].set_xlim(-lch,lch) + ax[4].set_ylim(-lch/20, lch+ max_scale*eu* lch + lch/20) + + # the strain history + the_line_hist_bg, = ax[0].plot(tim_noise, mod_noise, '--k', linewidth=1.0) + the_line_hist, = ax[0].plot(TIME, STRAIN, 'k', linewidth=1.0) + the_tip_hist, = ax[0].plot(TIME_TIP, STRAIN_TIP, 'or', fillstyle='full', markersize=8) + + # the composite response + the_line, = ax[1].plot(STRAIN, STRESS, '-k', linewidth=1.0) + the_tip, = ax[1].plot(STRAIN_TIP, STRESS_TIP, 'or', fillstyle='full', markersize=8) + + # the fracture + the_line_damage, = ax[2].plot(TIME, DAMAGE, '-b', linewidth=1.0) + + #equivalent plastic strain + the_line_eqp, = ax[3].plot(TIME, EQP, '-b', linewidth=1.0) + + #rve + y_shape = np.linspace(0,lch,100) + x_shape = np.zeros_like(y_shape) + points = np.array([x_shape, y_shape]).T.reshape(-1, 1, 2) + segments = np.concatenate([points[:-1], points[1:]], axis=1) + + # Initial uniform linewidth + linewidths = np.full(len(segments), radius*2.0) + + # Create LineCollection + the_line_rve = LineCollection(segments, colors='k', linewidths=linewidths) + ax[4].add_collection(the_line_rve) + # done + fig.tight_layout(pad=2.0) + + + + # process each step + time_start = time.time() + frames = [] + for i in range(Npulse): + failed = False + for j in range(Nsubsteps): + ok = ops.analyze(1) + if not do_plot: + if ok != 0: + failed = True + break + continue + if ok == 0: + # time + TIME.append(ops.getTime()) + # composite response + STRAIN.append(ops.eleResponse(1, 'material', 1, 'strain')[0]) + STRESS.append(ops.eleResponse(1, 'material', 1, 'stress')[0]) + # damage response + try: + DAMAGE.append(ops.eleResponse(1, 'material', 1, 'Damage')[0]) + except: + DAMAGE.append(0.0) + # equivalent plastic strain response + try: + EQP.append(ops.eleResponse(1, 'material', 1, 'EquivalentPlasticStrain')[0]) + except: + EQP.append(0.0) + # RVE response + try: + RVE.append(ops.eleResponse(1, 'material', 1, 'BucklingIndicator')[0]) + except: + RVE.append(0.0) + # update + the_line_hist.set_xdata(TIME) + the_line_hist.set_ydata(STRAIN) + the_line.set_xdata(STRAIN) + the_line.set_ydata(STRESS) + the_line_damage.set_xdata(TIME) + the_line_damage.set_ydata(DAMAGE) + the_line_eqp.set_xdata(TIME) + the_line_eqp.set_ydata(EQP) + A = RVE[-1] + strain = STRAIN[-1] + l_eff = lch * (1 + strain) + y_half = np.linspace(0, l_eff/2, 100) + x_half = np.zeros_like(y_half) + y_shape = np.concatenate([y_half, l_eff - y_half[::-1]]) + x_shape = np.concatenate([x_half, +x_half[::-1]]) + points = np.array([x_shape, y_shape]).T.reshape(-1, 1, 2) + segments = np.concatenate([points[:-1], points[1:]], axis=1) + + # Update the line collection + the_line_rve.set_segments(segments) + + n_segments = len(segments) + mid = n_segments // 2 + + sigma = n_segments / 20 + + # smooth profile (Maximum at the center and gradually decreases toward the sides) + x = np.arange(n_segments) + gaussian_profile = np.exp(-((x - mid) ** 2) / (2 * sigma ** 2)) + gaussian_profile = gaussian_profile / gaussian_profile.max() + linewidth_min = radius*2.0 + linewidth_max = max(radius*2.0 * (1.0 - DAMAGE[-1]), 0.5) + linewidths = linewidth_min + (linewidth_max - linewidth_min) * gaussian_profile + + the_line_rve.set_linewidths(linewidths) + else: + failed = True + break + if failed: + break + STRAIN_TIP = [STRAIN[-1]] + STRESS_TIP = [STRESS[-1]] + TIME_TIP = [TIME[-1]] + if i == Npulse-1: + the_tip.remove() + the_tip_hist.remove() + else: + the_tip.set_xdata(STRAIN_TIP) + the_tip.set_ydata(STRESS_TIP) + the_tip_hist.set_xdata(TIME_TIP) + the_tip_hist.set_ydata(STRAIN_TIP) + for iax in ax: + iax.relim() + iax.autoscale_view() + + fig.canvas.draw() + fig.canvas.flush_events() + plt.pause(0.001) + # store frame + w, h = fig.canvas.get_width_height() + buf = np.frombuffer(fig.canvas.buffer_rgba(), dtype=np.uint8) + buf.shape = (h, w, 4) + frame = buf[:, :, :3].copy() + + frames.append(frame) + time_end = time.time() + + plt.ioff() + plt.close(fig) + ops.wipe() + + print(f'Saving GIF animation {title}.gif...') + images = [Image.fromarray(frame) for frame in frames] + images[0].save('ASDSteel1D_Ex_Damage_Output.gif', format='GIF', save_all=True, append_images=images[1:], loop=0, duration=100, disposal=2) + + +# To generate load history +Npulse = 50 + +max_scale= 4.0 +eu = 0.06 + +mod_noise = np.linspace(0, eu * max_scale, Npulse) +mod_noise = list(mod_noise) + +tim_noise = [0.0] * Npulse +for i in range(1, Npulse): + dy = abs(mod_noise[i] - mod_noise[i-1]) + tim_noise[i] = tim_noise[i-1] + dy + +tim_noise_max = tim_noise[-1] +if tim_noise_max == 0.0: + tim_noise_max = 1.0 +tim_noise = [t / tim_noise_max for t in tim_noise] + + +title = f'results__damage_run' +run_simulation(tim_noise, mod_noise, max_scale, title, Npulse=Npulse) \ No newline at end of file diff --git a/source/user/manual/material/uniaxialMaterials/figures/ASDSteel1D/ASDSteel1D_Ex_Bond_Slip_Output.gif b/source/user/manual/material/uniaxialMaterials/figures/ASDSteel1D/ASDSteel1D_Ex_Bond_Slip_Output.gif new file mode 100644 index 00000000..71651b11 Binary files /dev/null and b/source/user/manual/material/uniaxialMaterials/figures/ASDSteel1D/ASDSteel1D_Ex_Bond_Slip_Output.gif differ diff --git a/source/user/manual/material/uniaxialMaterials/figures/ASDSteel1D/ASDSteel1D_Ex_Buckling_Output.gif b/source/user/manual/material/uniaxialMaterials/figures/ASDSteel1D/ASDSteel1D_Ex_Buckling_Output.gif new file mode 100644 index 00000000..936277e2 Binary files /dev/null and b/source/user/manual/material/uniaxialMaterials/figures/ASDSteel1D/ASDSteel1D_Ex_Buckling_Output.gif differ diff --git a/source/user/manual/material/uniaxialMaterials/figures/ASDSteel1D/ASDSteel1D_Ex_Damage_Output.gif b/source/user/manual/material/uniaxialMaterials/figures/ASDSteel1D/ASDSteel1D_Ex_Damage_Output.gif new file mode 100644 index 00000000..f4e4139a Binary files /dev/null and b/source/user/manual/material/uniaxialMaterials/figures/ASDSteel1D/ASDSteel1D_Ex_Damage_Output.gif differ diff --git a/source/user/manual/material/uniaxialMaterials/figures/ASDSteel1D/ASDSteel1D_Ex_numerical_vs_experimental.gif b/source/user/manual/material/uniaxialMaterials/figures/ASDSteel1D/ASDSteel1D_Ex_numerical_vs_experimental.gif new file mode 100644 index 00000000..c12cb6c2 Binary files /dev/null and b/source/user/manual/material/uniaxialMaterials/figures/ASDSteel1D/ASDSteel1D_Ex_numerical_vs_experimental.gif differ diff --git a/source/user/manual/material/uniaxialMaterials/figures/ASDSteel1D/disegnohomog.png b/source/user/manual/material/uniaxialMaterials/figures/ASDSteel1D/disegnohomog.png new file mode 100644 index 00000000..37f6f379 Binary files /dev/null and b/source/user/manual/material/uniaxialMaterials/figures/ASDSteel1D/disegnohomog.png differ