# Input file for Matrix Solution to Material Softening Bar

```/*
*  =======================================================================
*  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;
```

Developed in March 1998 by Mark Austin