/* * ======================================================================= * Compute force-displacement relationship for a material softening bar. * * -- One element remains elastic. * -- A second element has material softening behavior with Et = -0.1*E. * * We use the matrix language to implement the algorithm, and then later * on, compare to the 2D fiber finite element implementation. * * Written By : Wane-Jang Lin July 1997 * ======================================================================= */ /* Plastic spring material, Bi-linear */ L1 = 2 m; h1 = 30 cm; b1 = 10 cm; A1 = b1*h1; E1 = 30000 N/m^2; Et1 = -0.1*E1; fy1 = 1000 N/m^2; Ks1 = E1*A1/L1; Kt1 = Et1*A1/L1; Fy1 = fy1*A1; es1 = Fy1/Ks1; /* Elastic spring material */ L2 = 1.5 m; h2 = 20 cm; b2 = 10 cm; A2 = h2*b2; E2 = 20000 N/m^2; Et2 = E2; fy2 = 2500 N/m^2; Ks2 = E2*b2*h2/L2; /* Temporary matrix to store element tangent Ks at each step */ tangent = [Ks1 ; Ks2]; /* Initial condition, unstressed , matrix(2x1)=[element1, element2] */ p = [0 cm; 0 cm]; /* structure displacement */ PRe = [0 N; 0 N]; /* element resisting force */ PR = [0 N; 0 N]; /* structure resisting force */ Q_saved = [0 N; 0 N]; q_saved = [0 cm; 0 cm]; /* Transformation matrix Lele, d_q = Lele*d_p[ele] for each element */ /* Transform structural displacements d_p to element deformations d_q */ L = [-1 , 1]; /* Force interpolation matrices b(x), D(x) = b(x)*Q */ /* relate section force D(x) with element force Q */ bx = 1; /* Assemble initial structure tangent stiffness matrix BigK */ BigK = [ Ks1+Ks2, -Ks2; -Ks2, Ks2 ]; PrintMatrix ( BigK ); /* Initial displacement increment */ d_P = [ 0 N ; 1 N]; d_p = Solve( BigK, d_P ); PrintMatrix ( d_p ); /* Allocate response matrix[total_step+1][4] */ print "\n"; print "* response [][1] = resistant force *\n"; print "* response [][2] = total elongation *\n"; print "* response [][3] = element elongation 1 *\n"; print "* response [][4] = element elongation 2 *\n"; total_step = 45; response = ColumnUnits(Zero([total_step+1,4]),[N,cm,cm,cm]); response[1][1] = 0 N; response[1][2] = 0 cm; response[1][3] = 0 cm; response[1][4] = 0 cm; /* Increase displacement, structure determination */ index = 0; for ( k = 1; k <= total_step ; k = k+1 ) { p = p + d_p; /* State determination for each element */ for( ele = 1 ; ele <= 2; ele = ele+1 ) { /* extract element displacements from global displacements */ if( ele==1 ) { d_pe = [ 0 m; d_p[1][1] ]; } if( ele==2 ) { d_pe = [ d_p[1][1]; d_p[2][1] ]; } Q = Q_saved[ele][1]; /* retrieve values from (j-1) */ q = q_saved[ele][1]; Dx = bx*Q; dx = bx*q; K = tangent[ele][1]; kx = tangent[ele][1]; fx = 1/kx; rx = 0 cm; DUx = 1E+7 N; d_q = QuanCast(L*d_pe); /* element deformation increment */ q = q + d_q; /* convergence of forces for element "ele" */ while( abs(DUx) > 0.00001 N ) { d_Q = K*d_q; /* element force increment */ Q = Q + d_Q; /* update element force */ /* determine the section force increments */ /* repeat for all integration points of the element */ d_Dx = bx*d_Q; /* section force increment */ d_dx = rx + fx*d_Dx; /* section deformation increment */ Dx = Dx + d_Dx; /* update section force */ dx = dx + d_dx; /* update section deformation vector */ /* get new section tangent flexibility f(x) from new d(x) */ /* and section resisting force DR(x) from material properties */ if( ele == 1 ) { if( abs(dx) <= es1 ) { kx = Ks1; fx = 1/kx; DRx = kx*dx; } if( dx > es1 ) { kx = Kt1; fx = 1/kx; DRx = Fy1 + kx*(dx-es1); } if( dx < -es1 ) { kx = Kt1; fx = 1/kx; DRx = -Fy - kx*(dx+es1); } } if( ele == 2 ) { kx = Ks2; DRx = kx*dx; } DUx = Dx - DRx; /* section unbalanced force */ rx = fx*DUx; /* section residual deformation */ /* finish for all integration points of the element */ /* update the element flexibility and stiffness matrices */ F = bx*fx*bx; K = 1/F; /* check for element convergence */ s = bx*rx; /* element residual deformation */ d_q = -s; } /* j */ Q_saved[ele][1] = Q; q_saved[ele][1] = q; tangent[ele][1] = K; PRe[ele][1] = Q; /* element resisting force */ } /* Assemble structure resistant force */ PR[1][1] = PRe[1][1] - PRe[2][1]; PR[2][1] = PRe[2][1]; /* Store output in response matrix */ response[k+1][1] = PR[2][1]; response[k+1][2] = p[2][1]; response[k+1][3] = p[1][1]; response[k+1][4] = p[2][1] - p[1][1]; /* Calculate new delta displacement after yielding */ if( index == 1 ) { d_P[1][1] = 0 N; d_P[2][1] = -1 N; d_p = Solve( BigK, d_P ); index = 0; } /* Adjust delta displacement at yielding step */ if( abs(PR[1][1]) > 0.01 N ) { /* assemble new structure stiffness */ BigK[1][1] = tangent[1][1] + tangent[2][1]; BigK[1][2] = -tangent[2][1]; BigK[2][1] = -tangent[2][1]; BigK[2][2] = tangent[2][1]; index = 1; k = k-1; d_p[1][1] = 0 cm; d_p[2][1] = PR[1][1]/tangent[2][1]; } } PrintMatrix(response); /* End of Input */ print "\n"; print "===========================\n"; print "Nonlinear Analysis Complete\n"; print "===========================\n"; quit;