From b143c74d8fa66e6494b61528db7a4c9854018100 Mon Sep 17 00:00:00 2001 From: claudio perez <50180406+claudioperez@users.noreply.github.com> Date: Fri, 23 Aug 2024 10:21:35 -0700 Subject: [PATCH] cmp --- ...Concrete3D_Ex_CyclicUniaxialCompression.py | 125 ------------------ tests/openseespy/ASDConcrete3D_Ex_Surface.py | 113 ---------------- tests/openseespy/ASDConcrete3D_MakeLaws.py | 120 ----------------- 3 files changed, 358 deletions(-) delete mode 100644 tests/openseespy/ASDConcrete3D_Ex_CyclicUniaxialCompression.py delete mode 100644 tests/openseespy/ASDConcrete3D_Ex_Surface.py delete mode 100644 tests/openseespy/ASDConcrete3D_MakeLaws.py diff --git a/tests/openseespy/ASDConcrete3D_Ex_CyclicUniaxialCompression.py b/tests/openseespy/ASDConcrete3D_Ex_CyclicUniaxialCompression.py deleted file mode 100644 index 7807d2a656..0000000000 --- a/tests/openseespy/ASDConcrete3D_Ex_CyclicUniaxialCompression.py +++ /dev/null @@ -1,125 +0,0 @@ -from opensees import openseespy as ops -from ASDConcrete3D_MakeLaws import make as make_concrete -from matplotlib import pyplot as plt -import numpy as np - -# the 2D model -ops.wipe() -ops.model('basic', '-ndm', 2, '-ndf', 2) - -# the material (units = N, mm) -E = 30000.0 -v = 0.2 -fc = 30.0 -ft = fc/10.0 -ec = 2.0*fc/E -Gt = 0.073*fc**0.18 -Gc = 2.0*Gt*(fc/ft)**2 -Te, Ts, Td, Ce, Cs, Cd, lch_ref = make_concrete(E, ft, fc, ec, Gt, Gc) -ops.nDMaterial('ASDConcrete3D', 1, - E, v, # elasticity - '-Te', *Te, '-Ts', *Ts, '-Td', *Td, # tensile law - '-Ce', *Ce, '-Cs', *Cs, '-Cd', *Cd, # compressive law - '-autoRegularization', lch_ref # use auto regularization: the input Gt and Gc are NOT specific fracture energies - ) - -# the plane stress -ops.nDMaterial('PlaneStress', 2, 1) - -# a triangle with lch = 250 -lch = 250.0 -ops.node(1, 0, 0) -ops.node(2, lch, 0) -ops.node(3, 0, lch) -ops.element('tri31', 1, 1, 2, 3, 1.0, 'PlaneStress', 2) - -# fixity -ops.fix(1, 1, 1) -ops.fix(2, 0, 1) -ops.fix(3, 1, 0) - -# a simple ramp -ops.timeSeries('Linear', 1) - -# imposed strain -emax = 0.01 -cycles = np.linspace(0.0, emax, 10) - -# begin plot -SX = [0.0] -SY = [0.0] -PX = [0.0] -PY = [0.0] -plt.ion() -fig, ax = plt.subplots(1,1) -ax.grid(linestyle=':') -ax.set(xlim=[-emax*1.1, emax*0.1],ylim=[-fc*1.2, fc*0.1]) -ax.set_title('Cyclic uniaxial compression') -ax.set_xlabel('\N{GREEK SMALL LETTER EPSILON}\N{SUBSCRIPT ONE}\N{SUBSCRIPT ONE}') -ax.set_ylabel('\N{GREEK SMALL LETTER SIGMA}\N{SUBSCRIPT ONE}\N{SUBSCRIPT ONE}') -the_line, = ax.plot(SX, SY, '-k', linewidth=2.0) -the_tip, = ax.plot(PX, PY, 'or', fillstyle='full', markersize=8) - -# some default analysis settings -ops.constraints('Transformation') -ops.numberer('Plain') -ops.system('FullGeneral') -ops.test('NormDispIncr', 1.0e-8, 10, 0) -ops.algorithm('Newton') - -for icycle in range(1, len(cycles)): - # get current plastic strain (if any) - ep = -ops.nodeDisp(2)[0]/lch - # impose strain - current_strain = cycles[icycle] - ops.pattern('Plain', 1, 1) - ops.sp(2, 1, -current_strain*lch) - # load - # start from a percentage = ep/current_strain - time_start = ep/current_strain - ops.setTime(time_start) - num_incr = max(1, int((current_strain-ep)/emax*100.0)) - time_incr = (1.0-time_start)/float(num_incr) - ops.integrator('LoadControl', time_incr) - ops.analysis('Static') - for i in range(num_incr): - ok = ops.analyze(1) - if ok == 0: - strain = ops.eleResponse(1, 'material', 1, 'strain') - stress = ops.eleResponse(1, 'material', 1, 'stress') - SX.append(strain[0]) - SY.append(stress[0]) - the_line.set_xdata(SX) - the_line.set_ydata(SY) - PX = [SX[-1]] - PY = [SY[-1]] - the_tip.set_xdata(PX) - the_tip.set_ydata(PY) - fig.canvas.draw() - fig.canvas.flush_events() - else: - break - if ok != 0: - break - # unload - ops.remove('pattern', 1) - ops.setTime(0.0) - ops.integrator('LoadControl', 1.0) - ops.analysis('Static') - ok = ops.analyze(1) - # plot the unload point - strain = ops.eleResponse(1, 'material', 1, 'strain') - stress = ops.eleResponse(1, 'material', 1, 'stress') - SX.append(strain[0]) - SY.append(stress[0]) - the_line.set_xdata(SX) - the_line.set_ydata(SY) - PX = [SX[-1]] - PY = [SY[-1]] - the_tip.set_xdata(PX) - the_tip.set_ydata(PY) - fig.canvas.draw() - fig.canvas.flush_events() - -plt.ioff() -plt.show() diff --git a/tests/openseespy/ASDConcrete3D_Ex_Surface.py b/tests/openseespy/ASDConcrete3D_Ex_Surface.py deleted file mode 100644 index d8e839db64..0000000000 --- a/tests/openseespy/ASDConcrete3D_Ex_Surface.py +++ /dev/null @@ -1,113 +0,0 @@ -from opensees import openseespy as ops -from math import pi, sin, cos -from matplotlib import pyplot as plt - -def analyze_dir (dX, dY): - - # info - print("Analyze direction ({:8.3g}, {:8.3g})".format(dX, dY)) - - # the 2D model - ops.wipe() - ops.model('basic', '-ndm', 2, '-ndf', 2) - - # the isotropic material - E = 30000.0 - v = 0.2 - sig0 = 30.0 - # define a perfect bilinear behavior in tension and compression to record the failure surface - fc = sig0 - ec = fc/E - ft = fc/10.0 - et = ft/E - ops.nDMaterial('ASDConcrete3D', 1, E, v, - '-Ce', 0.0, ec, ec+1, - '-Cs', 0.0, fc, fc, - '-Cd', 0.0, 0.0, 0.0, - '-Te', 0.0, et, et+1, - '-Ts', 0.0, ft, ft, - '-Td', 0.0, 0.0, 0.0) - - # the plane stress - ops.nDMaterial('PlaneStress', 2, 1) - - # a triangle - ops.node(1, 0, 0) - ops.node(2, 1, 0) - ops.node(3, 0, 1) - ops.element('tri31', 1, 1, 2, 3, 1.0, 'PlaneStress', 2) - - # fixity - ops.fix(1, 1, 1) - ops.fix(2, 0, 1) - ops.fix(3, 1, 0) - - # a simple ramp - ops.timeSeries('Linear', 1, '-factor', 2.0*sig0) - - # imposed stresses in the current direction - ops.pattern('Plain', 1, 1) - ops.load(2, dX, 0.0) - ops.load(3, 0.0, dY) - - # analyze - ops.constraints('Transformation') - ops.numberer('Plain') - ops.system('FullGeneral') - ops.test('NormDispIncr', 1.0e-6, 10, 0) - ops.algorithm('Newton') - - dLambda = 0.1 - dLambdaMin = 0.0001 - Lambda = 0.0 - sX = 0.0 - sY = 0.0 - while True: - ops.integrator('LoadControl', dLambda) - ops.analysis('Static') - ok = ops.analyze(1) - if ok == 0: - stress = ops.eleResponse(1, 'material', 1, 'stress') - sX = stress[0] - sY = stress[1] - Lambda += dLambda - if Lambda > 0.9999: - break - else: - dLambda /= 2.0 - if dLambda < dLambdaMin: - break - - # done - return (sX, sY) - -# number of subdivisions -NDiv = 80 -NP = NDiv+1 -dAngle = 2.0*pi/NDiv -SX = [0.0]*NP -SY = [0.0]*NP - -plt.ion() -fig, ax = plt.subplots(1,1) -ax.set(xlim=[-45, 10],ylim=[-45, 10]) -ax.plot([-1000,1000],[0,0],color='black',linewidth=0.5) # S1 -ax.plot([0,0],[-1000,1000],color='black',linewidth=0.5) # S2 -ax.grid(linestyle=':') -ax.set_aspect('equal', 'box') - -the_line, = ax.plot(SX, SY, '-k', linewidth=2.0) -for i in range(NP): - angle = float(i)*dAngle - dX = cos(angle) - dY = sin(angle) - a,b = analyze_dir(dX, dY) - SX[i] = a - SY[i] = b - the_line.set_xdata(SX) - the_line.set_ydata(SY) - fig.canvas.draw() - fig.canvas.flush_events() - -plt.ioff() -plt.show() diff --git a/tests/openseespy/ASDConcrete3D_MakeLaws.py b/tests/openseespy/ASDConcrete3D_MakeLaws.py deleted file mode 100644 index 26cc9ed79f..0000000000 --- a/tests/openseespy/ASDConcrete3D_MakeLaws.py +++ /dev/null @@ -1,120 +0,0 @@ -import math - -def _bezier3(xi, x0, x1, x2, y0, y1, y2): - ''' - a quadratic bezier curve. gives the yi at xi, from 3 control points - ''' - A = x0 - 2.0 * x1 + x2 - B = 2.0 * (x1 - x0) - C = x0 - xi - if abs(A) < 1.0e-12: - x1 = x1 + 1.0E-6 * (x2 - x0) - A = x0 - 2.0 * x1 + x2 - B = 2.0 * (x1 - x0) - C = x0 - xi - if A == 0.0: - return 0.0 - D = B * B - 4.0 * A * C - t = (math.sqrt(D) - B) / (2.0 * A) - return (y0 - 2.0 * y1 + y2) * t * t + 2.0 * (y1 - y0) * t + y0 - -def _get_lch_ref(E,ft,Gt,fc,ec,Gc): - ''' - automatically computes a dummy reference characteristic length - when Gt and Gc are computed from equations - ''' - - # min lch for tension - et_el = ft/E # no damage at peak - Gt_min = ft*et_el/2.0 - hmin_t = Gt/Gt_min/100.0 - - # min lch for compression - ec1 = fc/E - ec_pl = (ec-ec1)*0.4 + ec1 - Gc_min = fc*(ec-ec_pl)/2.0 - hmin_c = Gc/Gc_min/100.0 - - # return the minimum - return min(hmin_t, hmin_c) - -def _make_tension(E, ft, Gt): - ''' - a trilinar hardening-softening law for tensile response - ''' - f0 = ft*0.9 - f1 = ft - e0 = f0/E - e1 = f1/E*1.5 - ep = e1-f1/E - f2 = ft/5.0 - f3 = 1.0e-3*ft - w2 = Gt/ft - w3 = 5.0*w2 - e2 = w2 + f2/E + ep - if e2 <= e1: e2 = e1*1.001 - e3 = w3 + f3/E + ep - if e3 <= e2: e3 = e2*1.001 - e4 = e3*10.0 - Te = [0.0, e0, e1, e2, e3, e4] # total strain points - Ts = [0.0, f0, f1, f2, f3, f3] # nominal stress points - Td = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0] # initialize damage list - Tpl = [0.0, 0.0, ep, e2*0.9, e3*0.8, e3*0.8] # desired values of equivalent plastic strains - for i in range(2, len(Te)): - xi = Te[i] - si = Ts[i] - xipl = Tpl[i] - xipl_max = xi-si/E - xipl = min(xipl, xipl_max) - qi = (xi-xipl)*E - Td[i] = 1.0-si/qi # compute damage - return (Te, Ts, Td) - -def _make_compression(E, fc, ec, Gc): - ''' - a quadratic hardening followed by linear softening for compressive response - ''' - fc0 = fc/2.0 - ec0 = fc0/E - ec1 = fc/E - fcr = fc*1.0e-1 - ec_pl = (ec-ec1)*0.4 + ec1 - Gc1 = fc*(ec-ec_pl)/2.0 - Gc2 = max(Gc1*1.0e-2, Gc-Gc1) - ecr = ec + 2.0*Gc2/(fc+fcr) - Ce = [0.0, ec0] # total strain points - Cs = [0.0, fc0] # nominal stress points - Cpl = [0.0, 0.0] # desired values of equivalent plastic strains - nc = 10 - dec = (ec-ec0)/(nc-1) - for i in range(nc-1): - iec = ec0+(i+1.0)*dec - Ce.append(iec) - Cs.append(_bezier3(iec, ec0, ec1, ec, fc0, fc, fc)) - Cpl.append(Cpl[-1]+(iec-Cpl[-1])*0.7) - # end of linear softening - begin residual plateau - Ce.append(ecr) - Cs.append(fcr) - Cpl.append(Cpl[-1] + (ecr-Cpl[-1])*0.7) - # extend to make a plateau - Ce.append(ecr+ec0) - Cs.append(fcr) - Cpl.append(Cpl[-1]) - # compute damage now - Cd = [0.0]*len(Ce) - for i in range(2, len(Ce)): - xi = Ce[i] - si = Cs[i] - xipl = Cpl[i] - xipl_max = xi-si/E - xipl = min(xipl, xipl_max) - qi = (xi-xipl)*E - Cd[i] = 1.0-si/qi # compute damage - # Done - return (Ce, Cs, Cd) - -def make(E, ft, fc, ec, Gt, Gc): - lch_ref = _get_lch_ref(E, ft, Gt, fc, ec, Gc) - Te, Ts, Td = _make_tension(E, ft, Gt/lch_ref) - Ce, Cs, Cd = _make_compression(E, fc, ec, Gc/lch_ref) - return (Te, Ts, Td, Ce, Cs, Cd, lch_ref) \ No newline at end of file