Input file for Analysis of Highway Bridge Structure.
Input file for Analysis of Highway Bridge Structure.
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