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.
*
* The element behavior is modeled with FIBER_2D finite elements.
*
* Written By : Wane-Jang Lin July 1997
* =======================================================================
*/
/* Define problem specific parameters */
NDimension = 2;
NDofPerNode = 3;
MaxNodesPerElement = 2;
GaussIntegPts = 2;
/* Define material and section properties for element 1 */
L1 = 2 m;
h1 = 30 cm; b1 = 10 cm;
E1 = 30000 N/m^2; Et1 = -0.1*E1; fy1 = 1000 N/m^2;
/* Define material and section properties for element 2 */
L2 = 1.5 m;
h2 = 20 cm; b2 = 10 cm;
E2 = 20000 N/m^2; Et2 = E2; fy2 = 2500 N/m^2;
/* Generate grid of nodes for finite element model */
total_node = 3;
total_elmt = 2;
StartMesh();
AddNode(1, [ 0 m, 0 m]);
AddNode(2, [ L1, 0 m]);
AddNode(3, [ L1+L2, 0 m]);
/* Attach elements to grid of nodes */
AddElmt( 1, [1, 2], "elmt_attr1");
AddElmt( 2, [2, 3], "elmt_attr2");
/* Define element, section, and material properties */
ElementAttr("elmt_attr1") { type = "FIBER_2D";
section = "section1";
material = "material1";
fiber = "fib_name1"; }
ElementAttr("elmt_attr2") { type = "FIBER_2D";
section = "section2";
material = "material2";
fiber = "fib_name2"; }
SectionAttr("section1") { width = b1; depth = h1; area = b1*h1; }
SectionAttr("section2") { width = b2; depth = h2; area = b2*h2; }
MaterialAttr("material1") { poisson = 0.25;
E = E1; Et = Et1; yield = fy1; }
MaterialAttr("material2") { poisson = 0.25;
E = E2; Et = Et2; yield = fy2; }
/* Define fiber attributes for element 1 */
no_fiber = 4; no_fiber_type = 1;
b = b1; h = h1; dh = h/no_fiber;
ks = E1; kt = Et1; fy = fy1;
fcoord = Matrix([ 1, no_fiber ]);
farea = Matrix([ 1, no_fiber ]);
fmap = Matrix([ 1, no_fiber ]);
for( i=1 ; i <= no_fiber ; i=i+1 ) {
fcoord[1][i] = h/2 - dh/2 - (i-1)*dh;
farea[1][i] = b*h/no_fiber;
fmap[1][i] = 1;
}
fattr = Matrix([ 3, 1 ]);
for( j=1 ; j <= no_fiber_type ; j=j+1 ) {
fattr[1][j] = ks;
fattr[2][j] = kt;
fattr[3][j] = fy;
}
FiberAttr( no_fiber, "fib_name1" ) { FiberMaterialAttr = fattr;
FiberCoordinate = fcoord;
FiberArea = farea;
FiberMaterialMap = fmap; }
/* Define fiber attributes for element 2 */
no_fiber = 4; no_fiber_type = 1;
b = b2; h = h2; dh = h/no_fiber;
ks = E2; kt = Et2; fy = fy2;
fcoord = Matrix([ 1, no_fiber ]);
farea = Matrix([ 1, no_fiber ]);
fmap = Matrix([ 1, no_fiber ]);
for( i=1 ; i <= no_fiber ; i=i+1 ) {
fcoord[1][i] = h/2 - dh/2 - (i-1)*dh;
farea[1][i] = b*h/no_fiber;
fmap[1][i] = 1;
}
fattr = Matrix([ 3, 1 ]);
for( j=1 ; j <= no_fiber_type ; j=j+1 ) {
fattr[1][j] = ks;
fattr[2][j] = kt;
fattr[3][j] = fy;
}
FiberAttr( no_fiber, "fib_name2" ) { FiberMaterialAttr = fattr;
FiberCoordinate = fcoord;
FiberArea = farea;
FiberMaterialMap = fmap; }
/* Setup Boundary Conditions */
FixNode(1, [1,1,1]);
FixNode(2, [0,1,1]);
FixNode(3, [0,1,1]);
EndMesh();
/* Get the initial stiffness, force, and system displacement */
Ks = Stiff();
NodeLoad( total_node, [ 0 N, 0 N, 0 N*m ] );
dP = ExternalLoad();
displ = Solve( Ks, dP );
/* Get d.o.f. in x- direction at nodes 2 and 3 */
load_dir = 1;
dof_matrix = GetDof([2]);
dof1 = dof_matrix[1][load_dir];
dof_matrix = GetDof([3]);
dof2 = dof_matrix[1][load_dir];
/* Allocate memory for response matrix [total_step+1][4] */
print " \n";
print "* Column [1] -- resistant force *\n";
print "* Column [2] -- total elongation *\n";
print "* Column [3] -- element elongation 1 *\n";
print "* Column [4] -- element elongation 2 *\n";
total_step = 45;
response = ColumnUnits( Zero([total_step+1,4]), [N,cm,cm,cm] );
response[1][1] = dP[dof2][1];
response[1][2] = displ[dof2][1];
response[1][3] = displ[dof1][1];
response[1][4] = displ[dof2][1] - displ[dof1][1];
/* Apply unit load and compute displacement increment */
NodeLoad( total_node, [ 1 N, 0 N, 0 N*m ] );
dP = ExternalLoad();
dp = Solve( Ks, dP );
NodeLoad( total_node, [ -1 N, 0 N, 0 N*m ] );
dP = ExternalLoad();
/* Incremental displacement added */
index = 0;
for( step=1 ; step <= total_step ; step = step + 1 ) {
/* Increment displacement by "dp" */
displ = displ + dp;
/* State determination for all elements */
ElmtStateDet( dp );
/* Assemble structure resistant force */
PR = InternalLoad( displ );
PrintMatrix(PR);
/* Update element information */
UpdateResponse();
/* Store displacements in "response" matrix */
response[step+1][1] = PR[dof2][1];
response[step+1][2] = displ[dof2][1];
response[step+1][3] = displ[dof1][1];
response[step+1][4] = displ[dof2][1] - displ[dof1][1];
/* Calculate new delta displacement after yielding */
if(index == 1 ) {
NodeLoad( total_node, [ -1 N, 0 N, 0 N*m ] );
dP = ExternalLoad();
dp = Solve( Ks, dP );
index = 0;
}
/* Adjust delta displacement at yielding step */
if( abs(PR[dof1][1]) > 0.01 N ) {
Ks = Stiff(); /* assemble new structure stiffness */
index = 1;
step = step - 1;
dp[dof1][1] = 0 m;
dp[dof2][1] = PR[dof1][1]/Ks[2][2];
}
}
PrintMatrix(response);
print "\n";
print "===========================\n";
print "Nonlinear Analysis Complete\n";
print "===========================\n";
quit;
Developed in March 1998 by Mark Austin
Last Modified March 5, 1998
Copyright © 1998, Mark Austin,
Department of Civil Engineering,
University of Maryland