Singularity problem - macro-element (unreinforced masonry)

[malaplante]

Hi everyone,

I'm trying to reproduce that macro-element :

O’Donnell, A., Kurama, Y., and Taflanidis, A. “An analytical modeling framework for seismic risk assessment of unreinforced masonry structures.”

but I have an instability in my model, but I don't know where (look at the end for the error messages). Here is my code :

wipe;
file mkdir Data;
model BasicBuilder -ndm 2 -ndf 3;

# Geometry
set Hwall [expr 2641.*$mm];
set Lwall [expr 3000.*$mm];
set ewall [expr 190.*$mm];
set hblock [expr 190.*$mm];
set egrout [expr 10.*$mm];
set efiber [expr $hblock+2.*$egrout]
set gamma [expr atan(($Lwall-2.*$efiber)/($Hwall-2.*$efiber))]
set Pdl [expr 855.*$kN];
set Awall [expr $Lwall*$ewall];
set Izwall [expr 1./12.*$ewall*pow($Lwall,3)];

# Material Tag
set ZLTSmatTag 1;
set DTEmatTag 2;
set FBCEmatTag 3;
set FBCEmasTag 11;
set ZLRSmatTag 4;

# ZERO-LENGTH TRANSLATIONAL SPRINGS ---------------------------------------------------------------
set Ksi [expr 1000.*$kip_in]
set Ksc [expr 0.09*$Ksi]

set c 0.;
set mu 0.65;
set Fsu [expr ($c+$mu*$Pdl/$Awall)*$Awall]
set Fsc [expr 0.5*$Fsu]

set s1p $Fsc
set e1p [expr $Fsc/$Ksi]
set s2p $Fsu
set e2p [expr $e1p+($Fsu-$Fsc)/$Ksc]
set s3p [expr 1.005*$Fsu]
set e3p [expr 10.*$e2p]

set s1n [expr -$s1p]
set e1n [expr -$e1p]
set s2n [expr -$s2p]
set e2n [expr -$e2p]
set s3n [expr -$s3p]
set e3n [expr -$e3p]

set pinchX 1.0;
set pinchY 1.0;
set damage1 0.0;
set damage2 0.0;
set beta 0.0;

uniaxialMaterial Hysteretic $ZLTSmatTag $s1p $e1p $s2p $e2p $s3p $e3p $s1n $e1n $s2n $e2n $s3n $e3n $pinchX $pinchY $damage1 $damage2 $beta

# DIAGONAL TRUSS ELEMENTS -------------------------------------------------------------------------
set G [expr 5700.*$MPa]
set Kdi [expr $G*$Awall/(2*pow(sin($gamma),2)*$Hwall)]
set Kdc [expr 0.5*$Kdi]

set b [expr $Hwall/$Lwall]
set ft [expr 1.*$MPa]
set tu [expr $ft/$b*pow($Pdl/($Awall*$ft)+1,0.5)]
set Vu [expr $tu*$Awall]

set Fdc [expr $Vu/(2*sin($gamma))]
set ftd [expr 2.3*$MPa]
set tu2 [expr $ftd/$b*pow($Pdl/($Awall*$ftd)+1,0.5)]
set Vu2 [expr $tu2*$Awall]
set Fdu [expr $Vu2/(2*sin($gamma))]

set s1p $Fdc
set e1p [expr $Fdc/$Kdi]
set s2p $Fdu
set e2p [expr $e1p+($Fdu-$Fdc)/$Kdc]
set s3p [expr 1.005*$Fdu]
set e3p [expr 10.*$e2p]

set s1n [expr -$s1p]
set e1n [expr -$e1p]
set s2n [expr -$s2p]
set e2n [expr -$e2p]
set s3n [expr -$s3p]
set e3n [expr -$e3p]

set pinchX 1.0;
set pinchY 1.0;
set damage1 0.0;
set damage2 0.0;
set beta 0.0;

uniaxialMaterial Hysteretic $DTEmatTag $s1p $e1p $s2p $e2p $s3p $e3p $s1n $e1n $s2n $e2n $s3n $e3n $pinchX $pinchY $damage1 $damage2 $beta

# FIBER BEAM-COLUMN ELEMENTS ----------------------------------------------------------------------
set fc1U [expr -17.5*$MPa]; # UNCONFINED concrete (todeschini parabolic model), maximum stress
set Ec [expr -15000.*$MPa]; # Concrete Elastic Modulus (the term in sqr root needs to be in psi
set eps1U -0.003; # strain at maximum strength of unconfined concrete
set fc2U [expr 0.2*$fc1U*$MPa]; # ultimate stress
set eps2U -0.0105; # strain at ultimate stress

uniaxialMaterial Concrete01 $FBCEmasTag $fc1U $eps1U $fc2U $eps2U;

set coverY [expr $Lwall/2.0]; # The distance from the section z-axis to the edge of the cover concrete -- outer edge of cover concrete
set coverZ [expr $ewall/2.0]; # The distance from the section y-axis to the edge of the cover concrete -- outer edge of cover concrete
set nfY 300;
set nfZ 1;
section Fiber $FBCEmatTag {;
patch rect $FBCEmasTag $nfY $nfZ -$coverY -$coverZ $coverY $coverZ;
}

# ZERO-LENGTH LINEAR-ELASTIC ROTATIONAL SPRINGS ---------------------------------------------------
uniaxialMaterial Elastic $ZLRSmatTag [expr 3.*$Ec*$Izwall/(0.5*$Hwall-$efiber)]

# Nodal coordinates --------------------------------------------------------------------------
# X Y
node 1 [expr -$Lwall/2+$efiber] $efiber
node 2 [expr $Lwall/2-$efiber] $efiber
node 3 [expr $Lwall/2-$efiber] [expr $Hwall-$efiber]
node 4 [expr -$Lwall/2+$efiber] [expr $Hwall-$efiber]
node 5 0 0
node 6 0 $efiber
node 7 0 [expr $Hwall-$efiber]
node 8 0 $Hwall

# Constraints ---------------------------------------------------------------------------
# node DX DY RZ
fix 5 1 1 1;
fix 1 0 0 0;
fix 2 0 0 0;
fix 3 0 0 0;
fix 4 0 0 0;
fix 6 0 0 0;
fix 7 0 0 0;
fix 8 0 0 0;


set TransfTag 1; # associate a tag to column transformation
geomTransf Linear $TransfTag;

# Element Tag --------------------------------------------------------------------------------
set ZLTSsecTag1 1
set ZLTSsecTag2 2
set DTEsecTag1 3
set DTEsecTag2 4
set FBCEsecTag1 5
set FBCEsecTag2 6
set ZLRSsecTag1 7

# RIGID LINKS -------------------------------------------------------------------------------------
equalDOF 6 1 1 2;
equalDOF 6 2 1 2;
equalDOF 7 3 1 2;
equalDOF 7 4 1 2;

# ZERO-LENGTH TRANSLATIONAL SPRINGS ---------------------------------------------------------------
element zeroLength $ZLTSsecTag2 7 8 -mat $ZLTSmatTag -dir 1

# DIAGONAL TRUSS ELEMENTS -------------------------------------------------------------------------
element twoNodeLink $DTEsecTag1 1 3 -mat $DTEmatTag -dir 1
element twoNodeLink $DTEsecTag2 2 4 -mat $DTEmatTag -dir 1

# FIBER BEAM-COLUMN ELEMENTS ----------------------------------------------------------------------
set numIntgrPts 3;
element forceBeamColumn $FBCEsecTag1 5 6 $numIntgrPts $FBCEmatTag $TransfTag
element forceBeamColumn $FBCEsecTag2 7 8 $numIntgrPts $FBCEmatTag $TransfTag

# ZERO-LENGTH LINEAR-ELASTIC ROTATIONAL SPRINGS ---------------------------------------------------
element zeroLength $ZLRSsecTag1 6 7 -mat $ZLRSmatTag -dir 6

pattern Plain 1 Linear {
load 8 0 -$Pdl 0
}

# Gravity-analysis parameters -- load-controlled static analysis
set Tol 1.0e-8; # convergence tolerance for test
constraints Lagrange; # how it handles boundary conditions
numberer Plain; # renumber dof's to minimize band-width (optimization), if you want to
system SparseGeneral ; # how to store and solve the system of equations in the analysis
test NormDispIncr $Tol 6 ; # determine if convergence has been achieved at the end of an iteration step
algorithm ModifiedNewton; # use Newton's solution algorithm: updates tangent stiffness at every iteration
set NstepGravity 10; # apply gravity in 10 steps
set DGravity [expr 1./$NstepGravity]; # first load increment;
integrator LoadControl $DGravity; # determine the next time step for an analysis
analysis Static; # define type of analysis static or transient
analyze $NstepGravity; # apply gravity

and I got as results :

WARNING SuperLU::solve(void)- Error 3 returned in factorization dgstrf()
WARNING ModifiedNewton::solveCurrentStep() -the LinearSysOfEqn failed in solve()

StaticAnalysis::analyze() - the Algorithm failed at iteration: 0 with domain at load factor -1
OpenSees analyze failed, returned: -3 error flag

[fmk]

as i have not a clue what it is you are trying to model .. i sugget fixing all nodes and gradualy releasing the fixities until you find the error.

[malaplante]

Thank for your hint. I have a question that might help me. I have brace in the macro-element with end release on each corner. I am using beam-column for the beam then I'm doing :

model BasicBuilder -ndm 2 -ndf 3;

then, do I have to restraint the rotation for each node of the truss elements who are connected with the node of beam with equalDOF ?