changed constitutive relation for every one time step

[cicilixiaoran]

Hello! Everyone, I want to simulate the changed elastic modulus and changed constitutive relation of concrete for every one time step. Here damage factor is the function of the concrete strain, and Ecxiu=(1-$dd)*$Ec, where Ecxiu is the elastic modulus after the concrete is damaged, dd is damage factor , Ec is initial elastic modulus. My code is as follows:

set IDconcUU 100;
set IDconcoverUU 1000;
set ColSecTagg 2000;
set numelement 3000;
set ok 0;
set controlTime [getTime];
while {$controlTime < $TmaxAnalysis && $ok == 0} {
set controlTime [getTime]
set ok [analyze 1 $DtAnalysis]
set bb [eleResponse 1 section 1 fiber -$ro 0.00 stressStrain]
set wordch [lindex $bb 1]
set word [expr abs($wordch)];
set dd [expr function($word)]; # dd is damage factor, and it is the function of the strain
set Ec [expr 4789*sqrt($fcu/$MPa)*$MPa];
set Ecxiu [expr (1-$dd)*$Ec];

incr IDconcUU ;
incr IDconcoverUU ;
incr ColSecTagg;
incr numelement;

# confined concrete
set Kfc 1.1;
set fc1C [expr $Kfc*$fc];
set eps1C [expr 2.*$fc1C/$Ecxiu];
set fc2C [expr 0.2*$fc1C];
set eps2C [expr 7*$eps1C];
# unconfined concrete
set fc1U $fc;
set eps1U [expr 2.*$fc1U/$Ecxiu];
set fc2U [expr 0.2*$fc1U];
set eps2U -0.03;

uniaxialMaterial Concrete01 $IDconcUU $fc1C $eps1C $fc2C $eps2C ;
uniaxialMaterial Concrete01 $IDconcoverUU $fc1U $eps1U $fc2U $eps2U ;

section fiberSec $ColSecTagg {
set rc [expr $ro-$coverSec];
patch circ $IDconcUU $nfCoreT $nfCoreR 0 0 $ri $rc 0 360;
patch circ $IDconcoverUU $nfCoverT $nfCoverR 0 0 $rc $ro 0 360;
set theta [expr 360.0/$numBarsSec];
layer circ $IDreinf $numBarsSec $barAreaSec 0 0 $rc $theta 360;
}

set numIntgrPts 5;
element nonlinearBeamColumn $numelement 1 2 $numIntgrPts $ColSecTagg $ColTransfTag;
}

Is the above code right or not? Thank you!

[vesna]

This will not do what you want. The only way of updating parameters during analysis would be to make a new material that will do what you want.

[cicilixiaoran]

Dear vesna,
Thank you very much.These days.I have tried to do it as you said,but i don't successfully do it.Do you mean that i should make a new material every one time step inside an analysis loop? Can you set an example for me ?Thank you very much.

Best regards

[vesna]

Initially I thought that you need to write a new material in C++ and add it to the OpenSees source code. Since you are using Concrete01 material you can do what you want by using "parameter" and "update parameter" commands.

Look at this post: http://opensees.berkeley.edu/community/ ... =2&t=46673. It will give you the idea how you can solve your problem.

The problem in this post is not exactly the same as yours but here you will find an example of how you can change parameters of your model at a certain step of analysis.

[cicilixiaoran]

Dear vesna,
Thank you very much. These days I carefully read the example you provided and wrote a program about cantilever column under cyclic load with changed elastic modulus for every one time step, here damage equation is the function of the concrete strain, and defined as: dd=0, when 0≤ε≤εf ; dd=((ε-εf)•εu) /((εu-εf)•ε), when εf≤ε≤εu ,where, εf is peak strain(selected as 0.0001), εu is limit strain(selected as 0.03). and Ecxiu=(1-$dd)*$Ec, where Ecxiu is the elastic modulus after the concrete is damaged, dd is damage factor , Ec is initial elastic modulus. My code is as follows:
wipe;
model BasicBuilder -ndm 2 -ndf 3;
set dataDir Data;
file mkdir $dataDir/;

set mm 1;
set kN 1;
set s 1;
set LunitTXT "mm";
set FunitTXT "kN";
set TunitTXT "sec";
set cm [expr 1e1*$mm];
set m [expr 1e2*$cm];
set mm2 [expr pow($mm,2.0)];
set N [expr $kN*1e-3];
set kNm [expr $kN*$m];
set Pa [expr $N/pow($m,2.0)];
set MPa [expr pow(10.,6)*$Pa];
set GPa [expr pow(10.,9)*$Pa];
set a [expr $m/pow($s,2.0)];
set g [expr 9.81*$a];
set kg [expr $N/$g];
set PI [expr 2.0*asin(1.0)];
set Ubig 1.0e10;
set Usmall [expr 1./$Ubig];

set LCol [expr 1.5*$m]; # column length
set Dcol [expr 0.15*$m];
set Rcol [expr $Dcol/2];
set ACol [expr $PI*pow($Rcol,2)];
set IzCol [expr 1./64*$PI*pow($Dcol,4)];

# nodal coordinates:
node 1 0 0; # node#, X, Y
node 2 0 $LCol

fix 1 1 1 1; # node DX DY RZ
fix 2 0 0 0; # node DX DY RZ

set IDctrlNode 2;
set IDctrlDOF 1;

set IDconcCore 1;
set IDconcCover 2;
set IDreinf 3;

set fcu [expr 30*$MPa];
set fc [expr -0.76*$fcu];
set Ec [expr 4789*sqrt($fcu/$MPa)*$MPa];

set fco [expr 0.76*pow(($Ec/4789/$MPa),2)*$MPa];
set xco [expr 2.*$fco/$Ec];

set n 0.1;
set axialLoad [expr $n*$fco*$PI*pow($Rcol,2)];

# confined concrete
set fyh [expr 239*$MPa]; # yield strength of stirrups
set rs 0.013; # ratio of volumne of hoop reinforcement
set Kfc [expr 1+($rs*$fyh)/$fco];
set fc1C [expr -$Kfc*$fco];
set eps1C [expr 2.*$fc1C/$Ec];
set fc2C [expr 0.2*$fc1C];
set eps2C [expr 5*$eps1C];
# unconfined concrete
set fc1U [expr -$fco];
set eps1U [expr -$xco];
set fc2U [expr 0.2*$fc1U];
set eps2U [expr 5*$eps1U];

set Fy [expr 366*$MPa];
set Es [expr 196*$GPa];
set Bs 0.01;
set R0 18;
set cR1 0.925;
set cR2 0.15;

uniaxialMaterial Concrete01 $IDconcCore $fc1C $eps1C $fc2C $eps2C;
uniaxialMaterial Concrete01 $IDconcCover $fc1U $eps1U $fc2U $eps2U;
uniaxialMaterial Steel02 $IDreinf $Fy $Es $Bs $R0 $cR1 $cR2;

# section GEOMETRY -------------------------------------------------------------

set coverSec [expr 25*$mm]; # Column cover
set Rcore [expr ($Rcol-$coverSec)];
set numBarsSec 6; # number of longitudinal bars
set Dbar [expr 10*$mm];
set barAreaSec [expr $PI*pow($Dbar,2)/4.0]; # area of bars
set ColSecTag 1;
set rsc [expr (6*$barAreaSec)/($PI*pow($Rcol,2))];

set ri 0.0; # inner radius
set ro [expr $Dcol/2]; # outer radius
set nfCoreR 8; # number of radial divisions in the core
set nfCoreT 8; # number of theta divisions in the core
set nfCoverR 4; # number of radial divisions in the cover
set nfCoverT 8; # number of theta divisions in the cover

section fiberSec $ColSecTag {
set rc [expr $ro-$coverSec];
patch circ $IDconcCore $nfCoreT $nfCoreR 0 0 $ri $rc 0 360;
patch circ $IDconcCover $nfCoverT $nfCoverR 0 0 $rc $ro 0 360;
set theta [expr 360.0/$numBarsSec];
layer circ $IDreinf $numBarsSec $barAreaSec 0 0 $rc $theta 360;
}

geomTransf Linear 1;

element nonlinearBeamColumn 1 1 2 5 $ColSecTag 1;

parameter 1 element 1 section $ColSecTag Ec;

recorder Node -file $dataDir/DFree.out -time -node 2 -dof 1 2 3 disp;
recorder Node -file $dataDir/RBase.out -time -node 1 -dof 1 2 3 reaction;
recorder Element -file $dataDir/DefoColSec3.out -time -ele 1 section 3 fiber -$ro 0.00 stressStrain;

pattern Plain 3001 "Constant" {
load 2 0.0 -$axialLoad 0.0
}

set Tol 1.0e-8;
constraints Plain;
numberer Plain;
system BandGeneral;
test NormDispIncr $Tol 6 ;
algorithm Newton;
integrator LoadControl 0 1 0 0;
analysis Static;
analyze 1;

set IDctrlNode 2;
set IDctrlDOF 1;

set qvfuweiyi [expr 5*$mm];

set iDmax "1 2 3 4 5 ";
set Dincr [expr 0.1*$qvfuweiyi];
set Fact $qvfuweiyi;
set CycleType Full;
set Ncycles 1;

set Hload [expr $axialLoad];
pattern Plain 200 Linear {;
load 2 $Hload 0.0 0.0 0.0 0.0 0.0
}
constraints Plain
numberer RCM
system BandGeneral
test EnergyIncr 1.e-8 6 0;
algorithm Newton;
integrator DisplacementControl $IDctrlNode $IDctrlDOF $Dincr
analysis Static

proc GeneratePeaks {Dmax {DincrStatic 0.01} {CycleType "Full"} {Fact 1} } {;
file mkdir data
set outFileID [open data/tmpDsteps.tcl w]
set Disp 0.
puts $outFileID "set iDstep { ";puts $outFileID $Disp;puts $outFileID $Disp;
set Dmax [expr $Dmax*$Fact];
if {$Dmax<0} {;
set dx [expr -$DincrStatic]
} else {
set dx $DincrStatic;
}
set NstepsPeak [expr int(abs($Dmax)/$DincrStatic)]
for {set i 1} {$i <= $NstepsPeak} {incr i 1} {;
set Disp [expr $Disp + $dx]
puts $outFileID $Disp;
}
if {$CycleType !="Push"} {
for {set i 1} {$i <= $NstepsPeak} {incr i 1} {;
set Disp [expr $Disp - $dx]
puts $outFileID $Disp;
}
if {$CycleType !="HalfCycle"} {
for {set i 1} {$i <= $NstepsPeak} {incr i 1} {;
set Disp [expr $Disp - $dx]
puts $outFileID $Disp;
}
for {set i 1} {$i <= $NstepsPeak} {incr i 1} {;
set Disp [expr $Disp + $dx]
puts $outFileID $Disp;
}
}
}
puts $outFileID " }";
close $outFileID
source data/tmpDsteps.tcl;
return $iDstep
}

set yingbf 0.0001 # peak strain
set yingbu 0.03 # limit strain
foreach Dmax $iDmax {
set iDstep [GeneratePeaks $Dmax $Dincr $CycleType $Fact];
for {set i 1} {$i <= $Ncycles} {incr i 1} {
set zeroD 0
set D0 0.0
foreach Dstep $iDstep {
set D1 $Dstep
set Dincr [expr $D1 - $D0]
integrator DisplacementControl $IDctrlNode $IDctrlDOF $Dincr
analysis Static
set bb [eleResponse 1 section 3 fiber -$ro 0.00 stressStrain]

set wordch [lindex $bb 1]
set word [expr abs($wordch)];
if {$word<$yingbf} {; # avoid the divide by zero
set dd 0
} else {
set dd [expr (($word-$yingbf)*$yingbu)/(($yingbu-$yingbf)*$word)];

}
set Ecxiu [expr (1-$dd)*$Ec];
updateParameter 1 $Ecxiu;

set outFileIDdamage [open data/damage.out "a"]
set outFileIDmodulus [open data/modulus.out "a"]
puts $outFileIDdamage $dd;
close $outFileIDdamage
puts $outFileIDmodulus $Ecxiu;
close $outFileIDmodulus
set ok [analyze 1]

if {$ok != 0} {
if {$ok != 0} {
puts "Trying Newton with Initial Tangent .."
test NormDispIncr 1.e-8 2000 0
algorithm Newton -initial
set ok [analyze 1]
test EnergyIncr 1.e-8 6 0
algorithm Newton
}
if {$ok != 0} {
puts "Trying Broyden .."
algorithm Broyden 8
set ok [analyze 1 ]
algorithm Newton
}
if {$ok != 0} {
puts "Trying NewtonWithLineSearch .."
algorithm NewtonLineSearch 0.8
set ok [analyze 1]
algorithm Newton
}
if {$ok != 0} {
return -1
}; # end if
}; # end if
set D0 $D1;
}; # end Dstep
}; # end i
}; # end of iDmaxCycl
Is the above code right or not? Thank you very much!

Best regards

[vesna]

Why don't you do a simple pushover analysis as a test to see if your coding it right.

[cicilixiaoran]

Thank you, I do a pushover analysis with changed elastic modulus for every one time step, the code is as follows:


wipe;
model BasicBuilder -ndm 2 -ndf 3;
set dataDir Data;
file mkdir $dataDir/;

set mm 1;
set kN 1;
set s 1;
set LunitTXT "mm";
set FunitTXT "kN";
set TunitTXT "sec";
set cm [expr 1e1*$mm];
set m [expr 1e2*$cm];
set mm2 [expr pow($mm,2.0)];
set N [expr $kN*1e-3];
set kNm [expr $kN*$m];
set Pa [expr $N/pow($m,2.0)];
set MPa [expr pow(10.,6)*$Pa];
set GPa [expr pow(10.,9)*$Pa];
set a [expr $m/pow($s,2.0)];
set g [expr 9.81*$a];
set kg [expr $N/$g];
set PI [expr 2.0*asin(1.0)];
set Ubig 1.0e10;
set Usmall [expr 1./$Ubig];

set LCol [expr 1.5*$m];
set Dcol [expr 0.15*$m];
set Rcol [expr $Dcol/2];
set ACol [expr $PI*pow($Rcol,2)];
set IzCol [expr 1./64*$PI*pow($Dcol,4)];



node 1 0 0;
node 2 0 $LCol

fix 1 1 1 1;
fix 2 0 0 0;

set IDctrlNode 2;
set IDctrlDOF 1;

set IDconcCore 1;
set IDconcCover 2;
set IDreinf 3;

set fcu [expr 30*$MPa];
set fc [expr -0.76*$fcu];
set Ec [expr 4789*sqrt($fcu/$MPa)*$MPa];

set fco [expr 0.76*pow(($Ec/4789/$MPa),2)*$MPa];
set xco [expr 2.*$fco/$Ec];

set n 0.1;
set axialLoad [expr $n*$fco*$PI*pow($Rcol,2)];

# confined concrete
set fyh [expr 239*$MPa]; # yield strength of stirrups
set rs 0.013; # ratio of volumne of hoop reinforcement
set Kfc [expr 1+($rs*$fyh)/$fco];
set fc1C [expr -$Kfc*$fco];
set eps1C [expr 2.*$fc1C/$Ec];
set fc2C [expr 0.2*$fc1C];
set eps2C [expr 5*$eps1C];
# unconfined concrete
set fc1U [expr -$fco];
set eps1U [expr -$xco];
set fc2U [expr 0.2*$fc1U];
set eps2U [expr 5*$eps1U];

set Fy [expr 366*$MPa];
set Es [expr 196*$GPa];
set Bs 0.01;
set R0 18;
set cR1 0.925;
set cR2 0.15;

uniaxialMaterial Concrete01 $IDconcCore $fc1C $eps1C $fc2C $eps2C;
uniaxialMaterial Concrete01 $IDconcCover $fc1U $eps1U $fc2U $eps2U;
uniaxialMaterial Steel02 $IDreinf $Fy $Es $Bs $R0 $cR1 $cR2;

# section GEOMETRY -------------------------------------------------------------
set coverSec [expr 25*$mm]; # Column cover
set Rcore [expr ($Rcol-$coverSec)];
set numBarsSec 6; # number of longitudinal bars
set Dbar [expr 10*$mm];
set barAreaSec [expr $PI*pow($Dbar,2)/4.0]; # area of bars
set ColSecTag 1;
set rsc [expr (6*$barAreaSec)/($PI*pow($Rcol,2))];

set ri 0.0; # inner radius
set ro [expr $Dcol/2]; # outer radius
set nfCoreR 8; # number of radial divisions in the core
set nfCoreT 8; # number of theta divisions in the core
set nfCoverR 4; # number of radial divisions in the cover
set nfCoverT 8; # number of theta divisions in the cover

section fiberSec $ColSecTag {
set rc [expr $ro-$coverSec];
patch circ $IDconcCore $nfCoreT $nfCoreR 0 0 $ri $rc 0 360;
patch circ $IDconcCover $nfCoverT $nfCoverR 0 0 $rc $ro 0 360;
set theta [expr 360.0/$numBarsSec];
layer circ $IDreinf $numBarsSec $barAreaSec 0 0 $rc $theta 360;
}

geomTransf Linear 1;

element nonlinearBeamColumn 1 1 2 5 $ColSecTag 1;

parameter 1 element 1 section $ColSecTag Ec;

recorder Node -file $dataDir/DFree.out -time -node 2 -dof 1 2 3 disp;
recorder Node -file $dataDir/RBase.out -time -node 1 -dof 1 2 3 reaction;
recorder Element -file $dataDir/DefoColSec3.out -time -ele 1 section 3 fiber -$ro 0.00 stressStrain;

pattern Plain 3001 "Constant" {
load 2 0.0 -$axialLoad 0.0
}

constraints Plain;
numberer Plain;
system BandGeneral;
test NormDispIncr 1.0e-8 6 ;
algorithm Newton;
integrator LoadControl 0 1 0 0;
analysis Static;
analyze 1;

set Dmax [expr 0.03*$LCol];
set Dincr [expr 0.001*$LCol];

set Hload [expr $axialLoad];
pattern Plain 3002 Linear {;
load 2 $Hload 0.0 0.0 0.0 0.0 0.0
}
constraints Plain
numberer RCM
system BandGeneral
test EnergyIncr 1.e-8 6 0;
algorithm Newton;
integrator DisplacementControl $IDctrlNode $IDctrlDOF $Dincr
analysis Static

set yingbf 0.0001
set yingbu 0.03

set ok 0;
set controlDisp 0.0;
set D0 0.0;
set Dstep [expr ($controlDisp-$D0)/($Dmax-$D0)]
while {$Dstep < 1.0 && $ok == 0} {
set controlDisp [nodeDisp $IDctrlNode $IDctrlDOF ]
set Dstep [expr ($controlDisp-$D0)/($Dmax-$D0)]
set bb [eleResponse 1 section 3 fiber -$ro 0.00 stressStrain]

set wordch [lindex $bb 1]
set word [expr abs($wordch)];

if {$word<$yingbf} {;
set dd 0
} else {
set dd [expr (($word-$yingbf)*$yingbu)/(($yingbu-$yingbf)*$word)];

}

set Ecxiu [expr (1-$dd)*$Ec];
updateParameter 1 $Ecxiu;

set outFileIDdamage [open data/damage.out "a"]
set outFileIDmodulus [open data/modulus.out "a"]
puts $outFileIDdamage $dd;
close $outFileIDdamage
puts $outFileIDmodulus $Ecxiu;
close $outFileIDmodulus
set ok [analyze 1]

if {$ok != 0} {
puts "Trying Newton with Initial Tangent .."
test NormDispIncr 1.e-8 2000 0
algorithm Newton -initial
set ok [analyze 1 ]
test EnergyIncr 1.e-8 6 0
algorithm Newton
}
if {$ok != 0} {
puts "Trying Broyden .."
algorithm Broyden 8
set ok [analyze 1 ]
algorithm Newton
}
if {$ok != 0} {
puts "Trying NewtonWithLineSearch .."
algorithm NewtonLineSearch .8
set ok [analyze 1 ]
algorithm Newton
}
};

Do you agree with that? Thank you very much！

[cicilixiaoran]

Dear vesna,Thank you very much.I do a pushover analysis with changed elastic modulus for every one time step, the code is as follows:

wipe;
model BasicBuilder -ndm 2 -ndf 3;
set dataDir Data;
file mkdir $dataDir/;

set mm 1;
set kN 1;
set s 1;
set LunitTXT "mm";
set FunitTXT "kN";
set TunitTXT "sec";
set cm [expr 1e1*$mm];
set m [expr 1e2*$cm];
set mm2 [expr pow($mm,2.0)];
set N [expr $kN*1e-3];
set kNm [expr $kN*$m];
set Pa [expr $N/pow($m,2.0)];
set MPa [expr pow(10.,6)*$Pa];
set GPa [expr pow(10.,9)*$Pa];
set a [expr $m/pow($s,2.0)];
set g [expr 9.81*$a];
set kg [expr $N/$g];
set PI [expr 2.0*asin(1.0)];
set Ubig 1.0e10;
set Usmall [expr 1./$Ubig];

set LCol [expr 1.5*$m];
set Dcol [expr 0.15*$m];
set Rcol [expr $Dcol/2];
set ACol [expr $PI*pow($Rcol,2)];
set IzCol [expr 1./64*$PI*pow($Dcol,4)];



node 1 0 0;
node 2 0 $LCol

fix 1 1 1 1;
fix 2 0 0 0;

set IDctrlNode 2;
set IDctrlDOF 1;

set IDconcCore 1;
set IDconcCover 2;
set IDreinf 3;

set fcu [expr 30*$MPa];
set fc [expr -0.76*$fcu];
set Ec [expr 4789*sqrt($fcu/$MPa)*$MPa];

set fco [expr 0.76*pow(($Ec/4789/$MPa),2)*$MPa];
set xco [expr 2.*$fco/$Ec];

set n 0.1;
set axialLoad [expr $n*$fco*$PI*pow($Rcol,2)];

# confined concrete
set fyh [expr 239*$MPa]; # yield strength of stirrups
set rs 0.013; # ratio of volumne of hoop reinforcement
set Kfc [expr 1+($rs*$fyh)/$fco];
set fc1C [expr -$Kfc*$fco];
set eps1C [expr 2.*$fc1C/$Ec];
set fc2C [expr 0.2*$fc1C];
set eps2C [expr 5*$eps1C];
# unconfined concrete
set fc1U [expr -$fco];
set eps1U [expr -$xco];
set fc2U [expr 0.2*$fc1U];
set eps2U [expr 5*$eps1U];

set Fy [expr 366*$MPa];
set Es [expr 196*$GPa];
set Bs 0.01;
set R0 18;
set cR1 0.925;
set cR2 0.15;

uniaxialMaterial Concrete01 $IDconcCore $fc1C $eps1C $fc2C $eps2C;
uniaxialMaterial Concrete01 $IDconcCover $fc1U $eps1U $fc2U $eps2U;
uniaxialMaterial Steel02 $IDreinf $Fy $Es $Bs $R0 $cR1 $cR2;

# section GEOMETRY -------------------------------------------------------------
set coverSec [expr 25*$mm]; # Column cover
set Rcore [expr ($Rcol-$coverSec)];
set numBarsSec 6; # number of longitudinal bars
set Dbar [expr 10*$mm];
set barAreaSec [expr $PI*pow($Dbar,2)/4.0]; # area of bars
set ColSecTag 1;
set rsc [expr (6*$barAreaSec)/($PI*pow($Rcol,2))];

set ri 0.0; # inner radius
set ro [expr $Dcol/2]; # outer radius
set nfCoreR 8; # number of radial divisions in the core
set nfCoreT 8; # number of theta divisions in the core
set nfCoverR 4; # number of radial divisions in the cover
set nfCoverT 8; # number of theta divisions in the cover

section fiberSec $ColSecTag {
set rc [expr $ro-$coverSec];
patch circ $IDconcCore $nfCoreT $nfCoreR 0 0 $ri $rc 0 360;
patch circ $IDconcCover $nfCoverT $nfCoverR 0 0 $rc $ro 0 360;
set theta [expr 360.0/$numBarsSec];
layer circ $IDreinf $numBarsSec $barAreaSec 0 0 $rc $theta 360;
}

geomTransf Linear 1;

element nonlinearBeamColumn 1 1 2 5 $ColSecTag 1;

parameter 1 element 1 section $ColSecTag Ec;

recorder Node -file $dataDir/DFree.out -time -node 2 -dof 1 2 3 disp;
recorder Node -file $dataDir/RBase.out -time -node 1 -dof 1 2 3 reaction;
recorder Element -file $dataDir/DefoColSec3.out -time -ele 1 section 3 fiber -$ro 0.00 stressStrain;

pattern Plain 3001 "Constant" {
load 2 0.0 -$axialLoad 0.0
}

constraints Plain;
numberer Plain;
system BandGeneral;
test NormDispIncr 1.0e-8 6 ;
algorithm Newton;
integrator LoadControl 0 1 0 0;
analysis Static;
analyze 1;

set Dmax [expr 0.03*$LCol];
set Dincr [expr 0.001*$LCol];

set Hload [expr $axialLoad];
pattern Plain 3002 Linear {;
load 2 $Hload 0.0 0.0 0.0 0.0 0.0
}
constraints Plain
numberer RCM
system BandGeneral
test EnergyIncr 1.e-8 6 0;
algorithm Newton;
integrator DisplacementControl $IDctrlNode $IDctrlDOF $Dincr
analysis Static

set yingbf 0.0001
set yingbu 0.03

set ok 0;
set controlDisp 0.0;
set D0 0.0;
set Dstep [expr ($controlDisp-$D0)/($Dmax-$D0)]
while {$Dstep < 1.0 && $ok == 0} {
set controlDisp [nodeDisp $IDctrlNode $IDctrlDOF ]
set Dstep [expr ($controlDisp-$D0)/($Dmax-$D0)]
set bb [eleResponse 1 section 3 fiber -$ro 0.00 stressStrain]

set wordch [lindex $bb 1]
set word [expr abs($wordch)];

if {$word<$yingbf} {;
set dd 0
} else {
set dd [expr (($word-$yingbf)*$yingbu)/(($yingbu-$yingbf)*$word)];

}

set Ecxiu [expr (1-$dd)*$Ec];
updateParameter 1 $Ecxiu;

set outFileIDdamage [open data/damage.out "a"]
set outFileIDmodulus [open data/modulus.out "a"]
puts $outFileIDdamage $dd;
close $outFileIDdamage
puts $outFileIDmodulus $Ecxiu;
close $outFileIDmodulus
set ok [analyze 1]

if {$ok != 0} {
puts "Trying Newton with Initial Tangent .."
test NormDispIncr 1.e-8 2000 0
algorithm Newton -initial
set ok [analyze 1 ]
test EnergyIncr 1.e-8 6 0
algorithm Newton
}
if {$ok != 0} {
puts "Trying Broyden .."
algorithm Broyden 8
set ok [analyze 1 ]
algorithm Newton
}
if {$ok != 0} {
puts "Trying NewtonWithLineSearch .."
algorithm NewtonLineSearch .8
set ok [analyze 1 ]
algorithm Newton
}
};

Do you agree with that? Thank you very much！