You will need Aladdin 2.0 to run this problem.
/* * ===================================================================== * Analysis of 4 span isolated bridge, 3 piers, 2 abutments, 5 isolators * * Modeling Assumptions: * * -- Bottoms of the abuntments' isolators and piers are fixed. * -- Deck is modeled as FRAME_3D elements. * -- Pier columns are modeled as FRAME_3D elements. * -- Isolators are modeled as FIBER_3D elements. * -- Perfect elastic-plastic for isolators. * -- Peak earthquake = 0.5*g. (rec.6) * -- 5% of structural damping is added, C=A0*M+A1+K. * * Written By : Wane-Jang Lin October 1997 * ===================================================================== */ NDimension = 3; NDofPerNode = 6; MaxNodesPerElement = 2; GaussIntegPts = 2; StartMesh(); H2 = 14 m; H3 = 7 m; H4 = 21 m; iso_h = 20 cm; /* Assemble mesh of nodes and elements */ total_node = 0; x = 0 m; y = 0 m; z = 0 m; for( node_no=total_node+1 ; x<=200m ; node_no=node_no+1 ) { AddNode( node_no, [x,y,z] ); x = x + 50 m; } total_node = node_no-1; x = 0 m; y = -iso_h; z = 0 m; for( node_no=total_node+1 ; x<=200m ; node_no=node_no+1 ) { AddNode( node_no, [x,y,z] ); x = x + 50 m; } total_node = node_no-1; x = 50 m; y = -iso_h-3.5m; z = 0 m; for( node_no=total_node+1 ; y>=-(iso_h+H2) ; node_no=node_no+1 ) { AddNode( node_no, [x,y,z] ); y = y - 3.5 m; } total_node = node_no-1; fix2 = total_node; x = 100 m; y = -iso_h-3.5m; z = 0 m; for( node_no=total_node+1 ; y>=-(iso_h+H3) ; node_no=node_no+1 ) { AddNode( node_no, [x,y,z] ); y = y - 3.5 m; } total_node = node_no-1; fix3 = total_node; x = 150 m; y = -iso_h-3.5m; z = 0 m; for( node_no=total_node+1 ; y>=-(iso_h+H4) ; node_no=node_no+1 ) { AddNode( node_no, [x,y,z] ); y = y - 3.5 m; } total_node = node_no-1; fix4 = total_node; total_elmt = 0; for( elmt_no=total_elmt+1 ; elmt_no<=4 ; elmt_no=elmt_no+1 ) { AddElmt( elmt_no, [elmt_no, elmt_no+1], "deck_attr" ); } AddElmt( 5, [1,6], "isolator_attr1" ); AddElmt( 6, [2,7], "isolator_attr2" ); AddElmt( 7, [3,8], "isolator_attr3" ); AddElmt( 8, [4,9], "isolator_attr4" ); AddElmt( 9, [5,10], "isolator_attr1" ); total_elmt = 9; node_no = fix2 - H2/(3.5m) + 1; AddElmt( total_elmt+1, [7, node_no], "column_attr2" ); total_elmt = total_elmt+1; for( elmt_no=total_elmt+1 ; node_no < fix2 ; elmt_no=elmt_no+1 ) { AddElmt( elmt_no, [node_no, node_no+1], "column_attr2" ); node_no = node_no+1; } total_elmt = elmt_no-1; bot_elmt2 = total_elmt; node_no = fix3 - H3/(3.5m) + 1; AddElmt( total_elmt+1, [8, node_no], "column_attr3" ); total_elmt = total_elmt+1; for( elmt_no=total_elmt+1 ; node_no < fix3 ; elmt_no=elmt_no+1 ) { AddElmt( elmt_no, [node_no, node_no+1], "column_attr3" ); node_no = node_no+1; } total_elmt = elmt_no-1; bot_elmt3 = total_elmt; node_no = fix4 - H4/(3.5m) + 1; AddElmt( total_elmt+1, [9, node_no], "column_attr4" ); total_elmt = total_elmt+1; for( elmt_no=total_elmt+1 ; node_no < fix4 ; elmt_no=elmt_no+1 ) { AddElmt( elmt_no, [node_no, node_no+1], "column_attr4" ); node_no = node_no+1; } total_elmt = elmt_no-1; bot_elmt4 = total_elmt; ElementAttr("deck_attr") { type = "FRAME_3D"; section = "deck_sec"; material = "deck_mat"; } SectionAttr("deck_sec") { Iyy = 103.77 m^4; Izz = 5.26 m^4; J = 16.13 m^4; area = 6.88 m^2; unit_weight = 200 kN/m; } MaterialAttr("deck_mat") { poisson = 0.25; E = 2.3e7 kN/m^2; } b = 2.2 m; h = 4.0 m; Ix = 7.39 m^4; Iy = 0.67 m^4; K2 = 30700 kN/m; K3 = 124400 kN/m; K4 = 13650 kN/m; E2 = K2*H2*H2*H2/3/Ix; E3 = K3*H3*H3*H3/3/Ix; E4 = K4*H4*H4*H4/3/Ix; ElementAttr("column_attr2") { type = "FRAME_3D"; section = "col_sec"; material = "col_mat2"; } ElementAttr("column_attr3") { type = "FRAME_3D"; section = "col_sec"; material = "col_mat3"; } ElementAttr("column_attr4") { type = "FRAME_3D"; section = "col_sec"; material = "col_mat4"; } SectionAttr("col_sec") { area = 4.16 m^2; width = h; depth = b; Iyy = Ix; Izz = Iy; unit_weight = 104 kN/m; shear_factor = 1.0; } MaterialAttr("col_mat2") { poisson = 0.25; E = E2; } MaterialAttr("col_mat3") { poisson = 0.25; E = E3; } MaterialAttr("col_mat4") { poisson = 0.25; E = E4; } h = 50 cm; b = 50 cm; E_lead = 356200 kN/m^2; G1 = 4240 kN/m^2; G2 = 12880 kN/m^2; G3 = 9200 kN/m^2; G4 = 37600 kN/m^2; fvy = 8400 kN/m^2; ElementAttr("isolator_attr1") { type = "FIBER_3D"; section = "iso_sec"; material = "iso_mat1"; fiber = "iso_fib"; } ElementAttr("isolator_attr2") { type = "FIBER_3D"; section = "iso_sec"; material = "iso_mat2"; fiber = "iso_fib"; } ElementAttr("isolator_attr3") { type = "FIBER_3D"; section = "iso_sec"; material = "iso_mat3"; fiber = "iso_fib"; } ElementAttr("isolator_attr4") { type = "FIBER_3D"; section = "iso_sec"; material = "iso_mat4"; fiber = "iso_fib"; } SectionAttr("iso_sec") { area = b*h; width = b; depth = h; unit_weight = 0.0001 kN/m; shear_factor = 1.0; } MaterialAttr("iso_mat1") { poisson = 0.25; G = G1; Gt = 0.001*G1; shear_yield = fvy/2; } MaterialAttr("iso_mat2") { poisson = 0.25; G = G2; Gt = 0.001*G2; shear_yield = fvy; } MaterialAttr("iso_mat3") { poisson = 0.25; G = G3; Gt = 0.001*G3; shear_yield = fvy; } MaterialAttr("iso_mat4") { poisson = 0.25; G = G4; Gt = 0.001*G4; shear_yield = fvy; } no_fiber = 4; iso_fcoord = Matrix([ 2, no_fiber ]); iso_farea = Matrix([ 1, no_fiber ]); iso_fmap = Matrix([ 1, no_fiber ]); iso_fcoord[1][1] = b/4; iso_fcoord[2][1] = h/4; iso_fcoord[1][2] = -b/4; iso_fcoord[2][2] = h/4; iso_fcoord[1][3] = b/4; iso_fcoord[2][3] = -h/4; iso_fcoord[1][4] = -b/4; iso_fcoord[2][4] = -h/4; for( i=1 ; i <= no_fiber ; i=i+1 ) { iso_farea[1][i] = b*h/no_fiber; iso_fmap[1][i] = 1; } iso_fattr = Matrix([ 3, 1 ]); iso_fattr[1][1] = E_lead; iso_fattr[2][1] = E_lead; iso_fattr[3][1] = fvy*1e6; FiberAttr( no_fiber, "iso_fib" ) { FiberMaterialAttr = iso_fattr; FiberCoordinate = iso_fcoord; FiberArea = iso_farea; FiberMaterialMap = iso_fmap; } /* Apply boundary conditions */ FixNode( 1, [1,0,0,1,0,0]); FixNode( 5, [1,0,0,1,0,0]); FixNode( 6, [1,1,1,1,1,1]); FixNode( 10, [1,1,1,1,1,1]); FixNode(fix2, [1,1,1,1,1,1]); FixNode(fix3, [1,1,1,1,1,1]); FixNode(fix4, [1,1,1,1,1,1]); /* Compile and print f.e. mesh */ EndMesh(); PrintMesh(); /* Compute initial stiffness and mass matrices */ stiff = Stiff(); mass = Mass([1]); mass_inv = Inverse( mass ); /* Compute and print eigenvalues */ no_eigen = 2; eigen = Eigen(stiff, mass, [no_eigen]); eigenvalue = Eigenvalue(eigen); for(i=1 ; i <= no_eigen ; i=i+1) { w = sqrt( eigenvalue[i][1] ); print "Mode",i,",\tw=",w,",\tperiod=",2*PI/w,"\n"; } eigenvector = Eigenvector(eigen); PrintMatrix(eigenvector); wa = sqrt( eigenvalue[1][1] ); wb = sqrt( eigenvalue[2][1] ); /* Setup Rayleigh damping and damping matrix */ rdamping = 0.20; A0 = 2*rdamping*wa*wb/(wa+wb); A1 = 2*rdamping/(wa+wb); damp = A0*mass + A1*stiff; /* El Centro Earthquake record */ Elcentro1 = ColumnUnits( [ 4.56; 13.77; 6.13; 3.73; 1.32; -6.81; -16.22; -22.41; -26.67; -24.10; -21.86; -14.02; -1.62; 1.67; 18.12; 39.26; 28.40; 16.72; 16.06; 12.54; 0.15; -8.60; -14.00; -12.01; -3.34; 5.14; 14.04; 16.36; 29.90; 43.09; 32.13; 17.28; 18.44; 22.31; 13.24; -0.91; -6.81; -15.20; -14.88; -9.54; -13.32; -10.60; 11.80; 10.20; -12.90; -22.88; 1.42; 21.73; 7.92; 6.47; 11.39; -8.29; -33.44; -26.25; -23.99; -31.05; -25.42; -26.97; -38.48; -35.01; -27.78; -33.35; -29.56; -13.60; -13.97; -11.63; -1.31; 13.86; 35.71; 27.62; 25.33; 46.45; 44.12; 28.02; 19.13; 12.57; 11.60; 12.29; 9.73; 10.38; 18.93; 29.43; 30.25; 26.23; 23.15; 22.29; 21.11; 14.84; 11.02; 14.74; 21.65; 26.66; 28.45; 30.34; 30.21; 28.82; 31.31; 26.66; 15.29; 9.10; 3.58; -6.65; -15.33; -18.05; -17.46; -11.68; -3.48; -5.25; -15.43; -25.38; -32.32; -29.67; -22.01; -17.97; -22.34; -30.51; -38.89; -45.29; -55.26; -67.84; -76.74; -84.58; -86.83; -86.63; -80.87; -69.83; -59.81; -54.51; -48.94; -39.94; -38.67; -40.13; -38.04; -34.78; -31.30; -21.29; -8.85; 3.58; 13.97; 30.88; 35.09; 8.84; -10.41; -16.30; -20.10; -3.69; 33.71; 51.69; 43.52; 27.04; 34.70; 38.30; 36.60; 29.85; -9.22; -23.98; -17.28; -14.20; 11.58; 47.99; 64.26; 70.36; 68.93; 62.67; 45.30; 28.43; 23.73; 28.51; 34.46; 43.12; 43.39; 26.64; 21.23; 31.23; 40.30; 16.76; -10.80; -22.88; -0.51; 12.51; -10.01; -12.52; -6.05; -24.18; -43.17; 7.01; 50.06; 54.19; 68.83; 69.26; 57.53; 8.73; -38.98; -43.59; -51.12; -41.86; -19.51; -15.78; -11.54; -2.18; -2.19; 20.84; 23.40; -0.91; -11.45; -22.94; -40.84; -53.17; -44.02; -21.96; -3.06; -0.49; -9.59; -11.83; -9.86; -7.02; 0.88; 24.00; 50.57; 56.01; 51.75; 52.30; 47.79; 16.77; -6.09; -20.13; -20.42; -4.35; -0.33; -6.33; 1.74; 19.90; 37.97; 53.36; 57.68; 53.32; 28.37; 7.54; 6.90; 4.73; -4.97; -21.97; -28.16; -23.25; -28.43; -18.57; -3.93; -16.60; -31.92; -35.08; -34.77; -30.85; -30.37; -6.83; 24.57; 27.39; 23.72; 29.30; 5.79; -35.35; -54.31; -57.00; -59.38; -19.00; 32.14; 7.01; -21.44; -21.26; -26.01; -35.66; -33.68; -24.09; -21.39; -26.29; -27.73; -23.91; -14.27; -6.41; -0.59; 4.21; 1.06; -4.59; -15.16; -18.12; -7.09; 2.41; 6.08; 1.20; -10.53; -8.69; -4.70; -13.05; -9.51; -0.67; -7.10; -7.19; 5.36; 9.80; 7.87; 12.63; 20.28; 22.44; 24.86; 24.31; 15.25; 6.54; 6.86; 11.08; 18.12; 23.64; 22.47; 23.76; 21.57; 19.80; 25.47; 27.94; 14.81; -6.94; -13.57; -7.80; -2.65; -3.97; -1.01; 0.20; 5.14; 2.31; -16.89; -18.18; -11.55; -12.10; -3.07; 6.16; 8.87; 19.72; 21.31; 16.45; 12.70; 4.15; -1.59; 0.57; 11.04; 21.39; 27.42; 15.82; -9.25; -23.67; -17.34; 0.67; 8.69; 11.01; 12.01; 11.12; 13.62; 11.31; 9.71; 17.68; 15.83; 5.59; 0.60; -2.53; -4.46; -4.09; -7.81; -15.06; -19.43; -15.49; -8.93; -5.44; -1.33; -1.35; 7.19; 12.52; 0.21; -4.89; 0.22; 11.51; 20.17; 18.62; 16.72; 13.74; 10.62; 11.14; 12.74; 13.34; 14.16; 14.59; 13.81; 10.94; 5.91; 5.74; 14.53; 19.56; 9.27; -6.75; -14.69; -14.25; -11.48; -3.96; 6.54; 10.59; 2.26; -12.97; -24.01; -25.73; -18.75; -6.99; -2.22; 0.17; 3.63; -1.70; -4.73; 6.13; 20.26; 22.36; 19.58; 9.14; -0.09; -0.23; -1.00; 4.82; 15.39; 19.75; 7.45; -0.99; 0.77; -1.64; -9.34; -14.34; -13.86; -6.61; 5.49; 6.30; 2.08; -2.88; -2.25; 2.99; 3.30; -1.57; -6.86; -8.44; -11.03; -16.68; -17.43; -12.39; -12.20; -15.43; -16.28; -17.99; -20.21; -21.58; -23.05; -23.79; -23.64; -21.60; -15.32; -5.72; 11.14; 28.30; 33.22; 30.27; 25.49; 18.45; 7.62; -1.66; -7.86; -10.04; -5.71; -5.16; -12.56; -19.39; -22.16; -18.18; -12.00; -3.60; 3.99; 3.89; 0.85; 4.01; 6.05; 0.68; -0.14; 1.50; 3.17; 4.96; -0.03; -7.20; -4.85; -0.60; -3.05; -6.01; -6.08; -6.40; -6.70; -5.70; -6.14; -8.41; -10.05; -12.35; -15.72; 0.00 ], [cm/sec/sec] ); dimen = Dimension(Elcentro1); x = dimen[1][1]; /* Interpolate accelerations */ divid_no = 16; Elcentro = ColumnUnits( Matrix([ divid_no*x, 1 ]),[in/sec/sec] ); Elcentro[1][1] = Elcentro1[1][1]/16; Elcentro[2][1] = Elcentro1[1][1]/8; Elcentro[3][1] = Elcentro1[1][1]*3/16; Elcentro[4][1] = Elcentro1[1][1]/4; Elcentro[5][1] = Elcentro1[1][1]*5/16; Elcentro[6][1] = Elcentro1[1][1]*3/8; Elcentro[7][1] = Elcentro1[1][1]*7/16; Elcentro[8][1] = Elcentro1[1][1]/2; Elcentro[9][1] = Elcentro1[1][1]*9/16; Elcentro[10][1]= Elcentro1[1][1]*5/8; Elcentro[11][1]= Elcentro1[1][1]*11/16; Elcentro[12][1]= Elcentro1[1][1]*3/4; Elcentro[13][1]= Elcentro1[1][1]*13/16; Elcentro[14][1]= Elcentro1[1][1]*7/8; Elcentro[15][1]= Elcentro1[1][1]*15/16; Elcentro[divid_no*x][1] = Elcentro1[x][1]; for( i=1 ; i < x; i=i+1 ) { Elcentro[divid_no*i][1] = Elcentro1[i][1]; Elcentro[divid_no*i+1][1] = ( Elcentro1[i][1]*15+Elcentro1[i+1][1] )/16; Elcentro[divid_no*i+2][1] = ( Elcentro1[i][1]*7 +Elcentro1[i+1][1] )/8; Elcentro[divid_no*i+3][1] = ( Elcentro1[i][1]*13+Elcentro1[i+1][1]*3 )/16; Elcentro[divid_no*i+4][1] = ( Elcentro1[i][1]*3 +Elcentro1[i+1][1] )/4; Elcentro[divid_no*i+5][1] = ( Elcentro1[i][1]*11+Elcentro1[i+1][1]*5 )/16; Elcentro[divid_no*i+6][1] = ( Elcentro1[i][1]*5 +Elcentro1[i+1][1]*3 )/8; Elcentro[divid_no*i+7][1] = ( Elcentro1[i][1]*9 +Elcentro1[i+1][1]*7 )/16; Elcentro[divid_no*i+8][1] = ( Elcentro1[i][1] +Elcentro1[i+1][1] )/2; Elcentro[divid_no*i+9][1] = ( Elcentro1[i][1]*7 +Elcentro1[i+1][1]*9 )/16; Elcentro[divid_no*i+10][1] = ( Elcentro1[i][1]*3 +Elcentro1[i+1][1]*5 )/8; Elcentro[divid_no*i+11][1] = ( Elcentro1[i][1]*5 +Elcentro1[i+1][1]*11 )/16; Elcentro[divid_no*i+12][1] = ( Elcentro1[i][1] +Elcentro1[i+1][1]*3 )/4; Elcentro[divid_no*i+13][1] = ( Elcentro1[i][1]*3 +Elcentro1[i+1][1]*13 )/16; Elcentro[divid_no*i+14][1] = ( Elcentro1[i][1] +Elcentro1[i+1][1]*7 )/8; Elcentro[divid_no*i+15][1] = ( Elcentro1[i][1] +Elcentro1[i+1][1]*15 )/16; } /* divid_no = 8; Elcentro = ColumnUnits( Matrix([ divid_no*x, 1 ]),[in/sec/sec] ); Elcentro[1][1] = Elcentro1[1][1]/8; Elcentro[2][1] = Elcentro1[1][1]/4; Elcentro[3][1] = Elcentro1[1][1]*3/8; Elcentro[4][1] = Elcentro1[1][1]/2; Elcentro[5][1] = Elcentro1[1][1]*5/8; Elcentro[6][1] = Elcentro1[1][1]*3/4; Elcentro[7][1] = Elcentro1[1][1]*7/8; Elcentro[divid_no*x][1] = Elcentro1[x][1]; for( i=1 ; i < x; i=i+1 ) { Elcentro[divid_no*i][1] = Elcentro1[i][1]; Elcentro[divid_no*i+1][1] = ( Elcentro1[i][1]*7+Elcentro1[i+1][1] )/8; Elcentro[divid_no*i+2][1] = ( Elcentro1[i][1]*3+Elcentro1[i+1][1] )/4; Elcentro[divid_no*i+3][1] = ( Elcentro1[i][1]*5+Elcentro1[i+1][1]*3 )/8; Elcentro[divid_no*i+4][1] = ( Elcentro1[i][1] +Elcentro1[i+1][1] )/2; Elcentro[divid_no*i+5][1] = ( Elcentro1[i][1]*3+Elcentro1[i+1][1]*5 )/8; Elcentro[divid_no*i+6][1] = ( Elcentro1[i][1] +Elcentro1[i+1][1]*3 )/4; Elcentro[divid_no*i+7][1] = ( Elcentro1[i][1] +Elcentro1[i+1][1]*7 )/8; } */ x = Max( [ Max(Elcentro), abs(Min(Elcentro)) ] ); Elcentro = (0.5*9.81/x)*Elcentro; print "\n*** Time History Analysis :"; print "\n*** Use Average Acceleration Method :\n\n"; /* ----------------------------------------------------- */ /* Setup Time History Analysis (and initial conditions) */ /* ----------------------------------------------------- */ /* Setup initial external load, internal load to zeros */ NodeLoad( 1, [ 0 kN, 0 kN, 0 kN, 0 kN*m, 0 kN*m, 0 kN*m] ); P_ext = ExternalLoad(); P_old = P_ext; /* Setup initial displacement, velocity and acceleration to zeros */ displ = Solve( stiff, P_ext ); velocity = displ/(1 sec); accel = velocity/(1 sec); new_displ = displ; new_velocity = velocity; new_accel = accel; displ_i = displ; /* Setup initial internal force and damping force */ Fs_old = InternalLoad( displ ); Fd_old = damp*velocity; Fs_i = Fs_old; /* Setup the influence vector */ r = displ/(1 m); accel_dir = 3; for( i=1 ; i<=total_node ; i=i+1 ) { dof = GetDof([i]); if( dof[1][accel_dir] > 0 ) { r[ dof[1][accel_dir] ][1] = 1; } } /* Setup time interval and analysis time */ dt = (0.02 sec)/divid_no; time = 0.0 sec; tol = 0.0001; stepno = 0; dimen = Dimension(Elcentro); total_stepno = dimen[1][1] + 250*divid_no; total_time = total_stepno * dt; quake_time = dimen[1][1] * dt; print " dt =",dt,"\n"; print "total step =",total_stepno," time =",total_time,"\n"; print "quake step =", dimen[1][1]," quake time =",quake_time,"\n"; /* Allocate response matrices */ deck_displ = ColumnUnits(Zero([total_stepno,5]),[m]); pier_displ = ColumnUnits(Zero([total_stepno,3]),[m]); bottom_shear = ColumnUnits(Zero([total_stepno,3]),[N]); isolator_shear = ColumnUnits(Zero([total_stepno,5]),[N]); energy = ColumnUnits(Zero([total_stepno,4]),[Jou]); /* Time-History Analysis */ while(stepno < total_stepno) { /* Update time and step no */ time = time + dt; stepno = stepno + 1; print "*** Start at step ", stepno, " : TIME = ", time, "\n"; /* Compute effective incremental loading */ if( stepno <= dimen[1][1] ) then { P_ext = -mass*r*Elcentro[stepno][1]; } else { P_ext = -mass*r*(0.0 m/sec/sec); } dPeff = P_ext - P_old + ((4/dt)*mass + 2*damp)*velocity + 2*mass*accel; /* while loop to check converge, Keff*U = Peff */ dp = displ - displ; err = 1 + tol; ii = 1; while( err > tol ) { /* Compute effective stiffness from tangent stiffness */ Keff = stiff + (2/dt)*damp + (4/dt/dt)*mass; /* Solve for d_displacement, d_velocity */ dp_i = Solve( Keff, dPeff); dp = dp + dp_i; dv = (2/dt)*dp - 2*velocity; /* Compute displacement, velocity */ new_displ = displ + dp; new_velocity = velocity + dv; /* Compute incremental displacement and internal load using old stiffness */ dFs = stiff*dp_i; Fs_i = Fs_i + dFs; if( ii==1 ) { x = L2Norm(dFs); if( x == 0 ) { x = 1; } } /* Check material yielding and compute new stiffness */ ElmtStateDet( dp_i ); stiff = Stiff(); /* Compute new internal load, damping force, and acceleration */ Fs = InternalLoad( new_displ ); Fd = damp*new_velocity; new_accel = mass_inv*( P_ext-Fs-Fd ); /* Calculate the unbalance force, and error percentage */ P_err = Fs_i - Fs; y = L2Norm(P_err); err = y/x; /* Assign new effective incremental load */ dPeff = P_err; displ_i = new_displ; Fs_i = Fs; print "*** In While Loop" ,ii, ", dFs = ", L2Norm(dFs), ", P_err =" , y , ", err =" ,err, "\n"; ii = ii+1; if( ii > 50 ) { flag = 1; err = tol; } } print "*** FINISHED CONVERGENCE AT STEP", stepno, "\n"; /* Calculate the energy balance, Wint(n+1) + T(n+1) = Wext(n+1) */ /* */ /* [1] column 1 is internal work by stiffness */ /* [2] column 2 is internal work by damping */ /* [3] column 3 is kinetic energy */ /* [4] column 4 is external work */ trans_vel = Trans(velocity); /* V(i-1) */ trans_vel_new = Trans(new_velocity); /* V(i) */ trans_dis_new = Trans(new_displ); /* D(i) */ if( stepno == 1 ) then { Wint_fs = dt/2*(trans_vel_new*Fs); Wint_fd = dt/2*(trans_vel_new*Fd); T = (trans_vel_new*mass*new_velocity)/2; Wext = dt/2*(trans_vel_new*P_ext); energy[1][1] = Wint_fs[1][1]; energy[1][2] = Wint_fd[1][1]; energy[1][3] = T[1][1]; energy[1][4] = Wext[1][1]; } else { Wint_fs = dt/2*(trans_vel*Fs_old + trans_vel_new*Fs); Wint_fd = dt/2*(trans_vel*Fd_old + trans_vel_new*Fd); T = (trans_vel_new*mass*new_velocity)/2; Wext = dt/2*(trans_vel*P_old + trans_vel_new*P_ext); energy[stepno][1] = energy[stepno-1][1] + Wint_fs[1][1]; energy[stepno][2] = energy[stepno-1][2] + Wint_fd[1][1]; energy[stepno][3] = T[1][1]; energy[stepno][4] = energy[stepno-1][4] + Wext[1][1]; } print "*** flag 1\n"; /* Tolerance is satisfied, update histories for this time step */ UpdateResponse(); P_old = P_ext; Fs_old = Fs; Fd_old = Fd; displ = new_displ; velocity = new_velocity; accel = new_accel; print "*** flag 2\n"; /* Store node displacement for this time step */ i = 1; for( node_no=1 ; node_no<=5 ; node_no=node_no+1 ) { dof = GetDof([node_no]); deck_displ[stepno][i] = displ[ dof[1][accel_dir] ][1]; i = i+1; } i = 1; for( node_no=7 ; node_no<=9 ; node_no=node_no+1 ) { dof = GetDof([node_no]); pier_displ[stepno][i] = displ[ dof[1][accel_dir] ][1]; i = i+1; } print "*** flag 3\n"; /* Store element shear force for this time step */ /* In the meantime, comment this section out ..... i = 1; for( elmt_no=5 ; elmt_no<=9 ; elmt_no=elmt_no+1 ) { stress = GetStress( [elmt_no], displ ); isolator_shear[stepno][i] = stress[1][accel_dir]; i=i+1; } */ print "*** flag 4\n"; /* stress = GetStress([bot_elmt2],displ); bottom_shear[stepno][1] = -stress[2][accel_dir]; stress = GetStress([bot_elmt3],displ); bottom_shear[stepno][2] = -stress[2][accel_dir]; stress = GetStress([bot_elmt4],displ); bottom_shear[stepno][3] = -stress[2][accel_dir]; */ print "*** flag 5\n"; if( flag == 1 ) { stepno = total_stepno; } print "deck_displ : "; print deck_displ[stepno][1],deck_displ[stepno][2],deck_displ[stepno][3], deck_displ[stepno][4],deck_displ[stepno][5],"\n"; print "pier_displ : "; print pier_displ[stepno][1], pier_displ[stepno][2], pier_displ[stepno][3],"\n"; print "isol_shear : "; print isolator_shear[stepno][1], isolator_shear[stepno][2], isolator_shear[stepno][3], isolator_shear[stepno][4], isolator_shear[stepno][5],"\n"; print "bottom_shear : "; print bottom_shear[stepno][1], bottom_shear[stepno][2], bottom_shear[stepno][3], "\n"; print "energy : "; print energy[stepno][1], energy[stepno][2], energy[stepno][3], energy[stepno][4],"\n"; } PrintMatrix(deck_displ,pier_displ); PrintMatrix(isolator_shear,bottom_shear); PrintMatrix(energy); quit;
Developed in March 1998 by Mark Austin
Copyright © 1998-1999, Department of Civil Engineering, University of Maryland