*heading
 Alfano Crisfield non-symmetric multi-delamination analysis
 3-dimensional model in ABAQUS/Standard using COH3D8 elements
*preprint, model=yes, echo=yes, history=yes
*parameter
** User input data defining the interface properties
** Ultimate strenth in tensile and mode II:
 ultI   =  3.30E6
 ultII  =  7.0E6
** ultI   =  66.0E6
** ultII  =  160.0E6
** cohesive layer modulus (N/m**2):
** Emod = 8.5E9
 Emod = 6.42E13
** Density
** density1 = 1.e3
** density2 = 0.75e3
 density1 = 2.0388e6
 density2 = 1.5291e6
** Fracture toughness (N/m):
 GIc   = 0.33e3
 GIIc  = 0.80e3
** B-K parameter:
 eta = 2.284
** thickness
 thick = 0.01325e-3
** width in the plane strain direction
 width = 20.0e-3
*parameter
 thickLayer = 0.1325e-3
 t1 = 10.0*thickLayer
 t2 =  2.0*thickLayer
 t3 = 12.0*thickLayer
 length = 200.0e-3
 z1 = t1
 z2 = z1 + thick
 z3 = z2 + t2
 z4 = z3 + thick
 z5 = z4 + t3
*node
 1001,       0.,   0.
 1201, <length>,   0.
 2001,       0., <z1>
 2201, <length>, <z1>
 3001,       0., <z2>
 3201, <length>, <z2>
 4001,       0., <z3>
 4201, <length>, <z3>
 5001,       0., <z4>
 5201, <length>, <z4>
 6001,       0., <z5>
 6201, <length>, <z5>
**
 11001,       0.,   0., <width>
 11201, <length>,   0., <width>
 12001,       0., <z1>, <width>
 12201, <length>, <z1>, <width>
 13001,       0., <z2>, <width>
 13201, <length>, <z2>, <width>
 14001,       0., <z3>, <width>
 14201, <length>, <z3>, <width>
 15001,       0., <z4>, <width>
 15201, <length>, <z4>, <width>
 16001,       0., <z5>, <width>
 16201, <length>, <z5>, <width>
*ngen, nset=nfront
 1001, 1201
 2001, 2201
 3001, 3201
 4001, 4201
 5001, 5201
 6001, 6201
*ngen, nset=nback
 11001, 11201
 12001, 12201
 13001, 13201
 14001, 14201
 15001, 15201
 16001, 16201
*nset, nset=nall
 nfront,nback
*nset, nset=right
 1201,
 2201,
 3201,
 4201,
 5201,
 6201,
 11201,
 12201,
 13201,
 14201,
 15201,
 16201, 
*element, type=C3D8I
 1001, 2001,2002,12002,12001, 1001,1002,11002,11001,
 3001, 4001,4002,14002,14001, 3001,3002,13002,13001, 
 5001, 6001,6002,16002,16001, 5001,5002,15002,15001, 
*elgen,elset=bottom
 1001, 200, 1, 1
*elgen,elset=middle
 3001, 200, 1, 1
*elgen,elset=top
 5001, 200, 1, 1
*elset,elset=all_solid
 bottom, middle, top
*solid section, elset=all_solid, mat=mat1, orient=global
**section controls,name=ehc,hourglass=enhanced
** 20.e-3,
*orientation, name=global
 1,0,0, 0,1,0
*material, name=mat1
*elastic, type=engineering
**    E1,    E2,    E3, nu12, nu13, nu23,   G12,   G13,
**   G23
 115.0E9, 8.5E9, 8.5E9, 0.29, 0.29,  0.3, 4.5E9, 3.3E9
 4.5E9,
*density
<density1>
*nset,nset=bondbot,generate
2001,2061
2081,2201
12001,12061
12081,12201
*nset,nset=bondtop,generate
5041,5201
15041,15201
*surface,name=bot_top
bottom,s1
*surface,name=mid_bot
middle,s2
*surface,name=mid_top
middle,s1
*surface,name=top_bot
top,s2
*boundary
** right, 2,2
 dispbc-bot, 1,1
 dispbc-top, 1,1
 nfront,3,3,
 nback,3,3,
*nset, nset=dispbc-bot
 1001,11001
*nset, nset=dispbc-top
 6001,16001
*nset, nset=dispbc
 1001, 6001
*contact pair,interaction=lower,adjust=bondbot
bot_top,mid_bot
*contact pair,interaction=upper,adjust=bondtop
top_bot,mid_top
*surface interaction,name=lower
*cohesive behavior,eligibility=original contacts
<Emod>,<Emod>,<Emod>
*DAMAGE INITIATION,CRITERION=QUADS
<ultI>, <ultII>, <ultII>
*DAMAGE EVOLUTION,TYPE=ENERGY,MIXED MODE BEHAVIOR=POWER LAW,POWER=1.0
<GIc>, <GIIc>, <GIIc>
*DAMAGE STABILIZATION
 1.e-6
*surface interaction,name=upper
*cohesive behavior,eligibility=original contacts
<Emod>,<Emod>,<Emod>
*DAMAGE INITIATION,CRITERION=QUADS
<ultI>, <ultII>, <ultII>
*DAMAGE EVOLUTION,TYPE=ENERGY,MIXED MODE BEHAVIOR=POWER LAW,POWER=1.0
<GIc>, <GIIc>, <GIIc>
*DAMAGE STABILIZATION
 1.e-6
**bulk viscosity
**0.09,1.8
**
**step, nlgeom=yes, UNSYMM=YES, inc=1000
*step, nlgeom=yes, inc=1000, extrapolation=parabolic
*static, stabilize, factor=7.e5
 0.01, 1.0, 1.e-8
*boundary
 dispbc-bot, 2,2, -20.0e-3
 dispbc-top, 2,2,  20.0e-3
*nset,nset=qa_test_nset
6001
*output,field,freq=9999
*node output
u,
*element output
s,
*output,field,freq=9999
*node output,nset=qa_test_nset
 u,
************************
*output,history,freq=1
*node output, nset=qa_test_nset
 u2,rf2
**energy output, variable=preselect
*monitor, dof=2, node=6001
*controls, parameter=field
 , 1.0
*controls, parameter=line search
 6, 5.0, 0.20,0.25,0.15
*controls, analysis=discontinuous
*controls, parameter=time incrementation
 , , , , , , , 10,
*end step