Moment-Curvature Analysis with ConfinedConcrete01

[marcoita]

Hi to everyone,

I am trying to do some moment-curvature analyses but I have a problem with the ConfinedConceret01 material. If I use, i.e., the Concrete01 material for the confined core the analysis is correct and I obtain what I expect. If I use, instead, the ConfinedConcrete01 for the confined core the analysis fails to converge and I don't understand why. The only thing I change between the two analyses is the material for the confined core.

Thank you for your help.

The following is the opensees file that doesn't converge.

#--------------------------PARAMETERS--------------------------

model BasicBuilder -ndm 2 -ndf 3

file mkdir Data;

#---geometric parameters--------
set al 0.3; # section height
set la 0.3; # section lenght
set cf 0.022; # concrete cover
#-------------------------------

# ---material properties-----
set fy 391300.; # [kPa]
set fc 24900.; # [kPa]
#-------------------------

#----fiber-------------
set ncf 2; # n° fiber in cover
set nynuc 25; # n° fiber along al
set nznuc 25; # n° fiber along la
set ny [expr $nynuc+(2*$ncf)];
set nz [expr $nznuc+(2*$ncf)];
#--------------------------------

#---bars diameter----------------
set f_16 0.016; # 16 mm
set A_16 [expr 3.14159*(pow($f_16,2)/4)];
set nb_up 0; # n° of superior rebars
set nb_dw 0; # n° of inferior rebars
#-----------------------------------------------

#---bars diameter----------------
set f_18 0.018; # 18 mm
set A_18 [expr 3.14159*(pow($f_18,2)/4)];
set nb_up1 3; # n° of superior rebars
set nb_dw1 3; # n° of inferior rebars

#-----hoop------------------------------------------
set f_st 0.008; # 8 mm
set s 0.1; # hooop spacing 10 cm
set nbracciax 2;
set nbracciay 2;

#---material identification number---------
set CP 7; # cover concrete (Kent-Park - Concrete01)
set Acc 3; # rebar steel
set BGL 14; # confined concrete (ConfinedConcrete01)
#----------------------------


#---Coordinate----
set z1 [expr (-$la/2)+$cf+$f_st/2];
set z2 $z1;
set z12 $z1;
set z11 $z1;

set z4 [expr (-$la/2)];
set z3 $z4;
set z5 [expr -$z4];
set z6 $z5;

set z8 [expr (-$z1)];
set z9 $z8;
set z10 $z8;
set z7 $z8;

set y1 [expr $al/2];
set y4 $y1;
set y8 $y1;
set y5 $y1;

set y3 [expr -$al/2];
set y2 $y3;
set y6 $y3;
set y7 $y3;

set y9 [expr $y1-$cf-$f_18/2-$f_st];
set y12 $y9;
set y10 [expr $y7+$cf+$f_18/2+$f_st];
set y11 $y10;

set y13 [expr $y1-$cf-$f_st/2]
set y15 [expr $y7+$cf+$f_st/2]

set y14 $y13
set y16 $y15

set z13 $z1
set z14 $z8

set z14 $z8
set z15 $z14

set z16 $z13

#---------------------MATERIALS-----------------------

#----STEEL_01:BILINEARE---------------------------------------------------------------3
set Fy $fy; # [kPa]
set Ea 206000000.; # [kPa]
set epsy [expr $Fy/$Ea];
set b 0.;
uniaxialMaterial Steel01 3 $Fy $Ea $b;
#--------------------------------------------------------------------------------------

#----CONFINEDCONCRETE_01:----------------------------BGL---------------------------
set L1 [expr $al-2.*$cf-$f_st]

# uniaxialMaterial ConfinedConcrete01 $tag $secType $fpc $Ec -epscu $epscu $nu $L1 $phis $S $fyh $Es0 $haRatio $mu $phiLon -stRatio $stRatio
uniaxialMaterial ConfinedConcrete01 14 S1 $fc 30000000. -epscu -0.03 -varub $L1 $f_st $s $fy 210000000. 0.0 1000.0 $f_18 -stRatio 0.85
#-----------------------------------------------------------------------------------------------------

#----CONCRETE_01:KENT-SCOTT-PARK---------------------UNCONFINED CONCRETE----------------------------------7
# ---Confined core dimension----
set hc [expr $al-2*$cf]
set bc [expr $la-2*$cf]
#------------------------

# ----ro-------
set ro 0
#---------------------

#----fcdCC ---------
set K [expr 1+$ro*$fy/$fc]
puts K=$K

set fcd [expr $K*$fc]
set gammac 1.5
set fcdCP [expr $fcd/$gammac]
puts fcdCP=$fcdCP
#----------------------

#----fcuCC ---------
set fcuCP [expr 0.2*$fcdCP]
puts fcuCP=$fcuCP
#----------------------

#----epsc0CC---------
set epsc0 0.002
set epsc0CP [expr $K*$epsc0]
puts epsc0CP=$epsc0CP
#-----------------------------

#----epsuCC--------------
set fyZ [expr $fy/1000.]; # [MPa]
set fcZ [expr $fc/1000.]; # [MPa]
set h2 [expr $hc-$f_st]
set Z [expr 0.5/(((3. + 0.29*$fcZ)/(145.*$fcZ - 1000.))+((3./4.)*$ro*(pow(($h2/$s),0.5)))-$epsc0CP)]
puts Z=$Z

set epsu1 [expr 0.004+0.9*$ro*($fyZ/300.)]; # Scott (1982)
set epsu2 [expr 0.8/$Z+$epsc0CP]; # Mander

if {$epsu1>$epsu2} {
set epsuCP $epsu1
} else {
if {$epsu1 < $epsu2} {
set epsuCP $epsu2
}
}
puts epsuCP=$epsuCP
#-------------------------------------
uniaxialMaterial Concrete01 7 $fcdCP $epsc0CP $fcuCP $epsuCP;
#--------------------------------------------------------------------------------------------------------

#--------SECTION----------------------
set sezRett 3;

section fiberSec $sezRett {
patch quad $BGL $nynuc $nznuc $y13 $z13 $y14 $z14 $y15 $z15 $y16 $z16
patch quad $CP $ny $ncf $y5 $z5 $y6 $z6 $y7 $z7 $y8 $z8
patch quad $CP $ny $ncf $y1 $z1 $y2 $z2 $y3 $z3 $y4 $z4
patch quad $CP $ncf $nznuc $y8 $z8 $y14 $z14 $y13 $z13 $y1 $z1
patch quad $CP $ncf $nznuc $y15 $z15 $y7 $z7 $y2 $z2 $y16 $z16

layer straight $Acc $nb_up $A_16 $y9 $z9 $y12 $z12
layer straight $Acc $nb_up1 $A_18 $y9 $z9 $y12 $z12
layer straight $Acc $nb_dw $A_16 $y10 $z10 $y11 $z11
layer straight $Acc $nb_dw1 $A_18 $y10 $z10 $y11 $z11
}

#----------------RECORDER------------------------------------
# Create recorder
recorder Node -file Data/Moment_curvature.out -time -node 1002 -dof 3 disp; # output moment (col 1) & curvature (col 2)

#-----------------------ANALYSIS-----------------------------
proc MomentCurvature2D { sezRett axialLoad maxK {numIncr 100} } {
##################################################
# A procedure for performing section analysis (only does
# moment-curvature, but can be easily modified to do any mode
# of section reponse.)
#
# MHS
# October 2000
# modified to 2D and to improve convergence by Silvia Mazzoni, 2006
#
# Arguments
# secTag -- tag identifying section to be analyzed
# axialLoad -- axial load applied to section (negative is compression)
# maxK -- maximum curvature reached during analysis
# numIncr -- number of increments used to reach maxK (default 100)
#
# Sets up a recorder which writes moment-curvature results to file
# section$secTag.out ... the moment is in column 1, and curvature in column 2

# Define two nodes at (0,0)
node 1001 0.0 0.0
node 1002 0.0 0.0

# Fix all degrees of freedom except axial and bending
fix 1001 1 1 1
fix 1002 0 1 0

# Define element
# tag ndI ndJ secTag
element zeroLengthSection 2001 1001 1002 $sezRett


# Define constant axial load
pattern Plain 3001 "Constant" {
load 1002 $axialLoad 0.0 0.0
}

# Define analysis parameters
integrator LoadControl 0 1 0 0
system SparseGeneral -piv; # Overkill, but may need the pivoting!
test EnergyIncr 1.0e-9 10
numberer Plain
constraints Plain
algorithm Newton
analysis Static

# Do one analysis for constant axial load
analyze 1

# Define reference moment
pattern Plain 3002 "Linear" {
load 1002 0.0 0.0 1.0
}

# Compute curvature increment
set dK [expr $maxK/$numIncr]

# Use displacement control at node 1002 for section analysis, dof 3
integrator DisplacementControl 1002 3 $dK 1 $dK $dK

# Do the section analysis
set ok [analyze $numIncr]

# ----------------------------------------------if convergence failure-------------------------
set IDctrlNode 1002
set IDctrlDOF 3
set Dmax $maxK
set Dincr $dK
set TolStatic 1.e-9;
set testTypeStatic EnergyIncr
set maxNumIterStatic 6
set algorithmTypeStatic Newton
if {$ok != 0} {
# if analysis fails, we try some other stuff, performance is slower inside this loop
set Dstep 0.0;
set ok 0
while {$Dstep <= 1.0 && $ok == 0} {
set controlDisp [nodeDisp $IDctrlNode $IDctrlDOF ]
set Dstep [expr $controlDisp/$Dmax]
set ok [analyze 1]; # this will return zero if no convergence problems were encountered
if {$ok != 0} {; # reduce step size if still fails to converge
set Nk 4; # reduce step size
set DincrReduced [expr $Dincr/$Nk];
integrator DisplacementControl $IDctrlNode $IDctrlDOF $DincrReduced
for {set ik 1} {$ik <=$Nk} {incr ik 1} {
set ok [analyze 1]; # this will return zero if no convergence problems were encountered
if {$ok != 0} {
# if analysis fails, we try some other stuff
# performance is slower inside this loop global maxNumIterStatic; # max no. of iterations performed before "failure to converge" is ret'd
puts "Trying Newton with Initial Tangent .."
test NormDispIncr $TolStatic 2000 0
algorithm Newton -initial
set ok [analyze 1]
test $testTypeStatic $TolStatic $maxNumIterStatic 0
algorithm $algorithmTypeStatic
}
if {$ok != 0} {
puts "Trying Broyden .."
algorithm Broyden 8
set ok [analyze 1 ]
algorithm $algorithmTypeStatic
}
if {$ok != 0} {
puts "Trying NewtonWithLineSearch .."
algorithm NewtonLineSearch 0.8
set ok [analyze 1]
algorithm $algorithmTypeStatic
}
if {$ok != 0} {; # stop if still fails to converge
puts [format $fmt1 "PROBLEM" $IDctrlNode $IDctrlDOF [nodeDisp $IDctrlNode $IDctrlDOF] $LunitTXT]
return -1
}; # end if
}; # end for
integrator DisplacementControl $IDctrlNode $IDctrlDOF $Dincr; # bring back to original increment
}; # end if
}; # end while loop
}; # end if ok !0
# -----------------------------------------------------------------------------------------------------
global LunitTXT; # load time-unit text
if { [info exists LunitTXT] != 1} {set LunitTXT "Length"}; # set blank if it has not been defined previously.

set fmt1 "%s Pushover analysis: CtrlNode %.3i, dof %.1i, Curv=%.4f /%s"; # format for screen/file output of DONE/PROBLEM analysis
if {$ok != 0 } {
puts [format $fmt1 "PROBLEM" $IDctrlNode $IDctrlDOF [nodeDisp $IDctrlNode $IDctrlDOF] $LunitTXT]
} else {
puts [format $fmt1 "Analisi terminata!" $IDctrlNode $IDctrlDOF [nodeDisp $IDctrlNode $IDctrlDOF] $LunitTXT]
}

}



#--------------------------LOADS----------------------------
# [F]=[kN]
# set AXIAL LOAD --------------------------------------------------------
set axialLoad -400; # + Tension, - Compression

# set maximum Curvature:
set Ku 0.08;
set numIncr 300; # Number of analysis increments to maximum curvature (default=100)

# Call the section analysis procedure
MomentCurvature2D $sezRett $axialLoad $Ku $numIncr

[marcoita]

I moved up the part of the analysis from "global LunitTXT..."; in this way I can solve the error regarding "fmt1" but the analysis does not converge yet and the message "WARNING SuperLU::solve(void)- Error 1 returned in factorization dgstrf()" is still there!