ContactMaterial2D and divide-by-zero errorr

[barkhordari]

I try to model soil and structure. I don't know how joint soil nodes to structure nodes. I use ContactMaterial2D but opensees shown below error:

:operator/<double fact> - divide-by-zero errorr coming.

I would be grateful if you could help me to solve it.

[barkhordari]

wipe


#---SOIL GEOMETRY
set nssh 3; ##### number of same section in height
set h4 4;
set numst 9; ## number of story
set bay 3.0;
set underf 3.0;
set underh 3.5 ;
set len 200.0 ;
set h1 10.0;
set h2 10.0;
set h3 15.0;
set bs 18.0 ;
set fw 2.0;
set df [expr ($underf*$underh)] ;
set dim1 2.0 ;
set dim2 1.0 ;
set dim3 0.5 ;
set rho 1.5;
set Vs 350;
set nu 0;

set G [expr $rho*$Vs*$Vs]
set E [expr 2*$G*(1+$nu)]
set bulk [expr $E/(3*(1-2*$nu))]

# soil cohesion (kPa)
set cohesion 95.0
# peak shear strain
set gammaPeak 0.1
# soil friction angle
set phi 0.0
# reference pressure
set refPress 80.0
# pressure dependency coefficient
set pressCoeff 0.0





model BasicBuilder -ndm 2 -ndf 2



#### CREATE LAGRANGE MULTIPLIER NODES FOR BEAM CONTACT ELEMENTS
#-----------------------------------------------------------------------------------------

for {set k 1} {$k <=((($df*2)/$dim3))} {incr k 1} {
node [expr 900000+$k] 0.00 0.00
}



#build node layer1
set count 0;
for {set j 1} {$j<=[expr ($h1/$dim1)+1]} {incr j 1} {
set count [expr $j-1]
for {set i 1} {$i<=[expr ($len/$dim1)+1]} {incr i 1} {
node [expr (int(100000+$i+($j-1)*(($len/$dim1)+1)))] [expr ((-$len/2)+($i-1)*$dim1)] [expr (0+$count*$dim1)]

}
};


#build node layer2

set count 0;
for {set j 1} {$j<=[expr ($h2/$dim2)+1]} {incr j 1} {
set count [expr $j-1]
for {set i 1} {$i<=[expr ($len/$dim2)+1]} {incr i 1} {
node [expr (int(200000+$i+($j-1)*(($len/$dim2)+1)))] [expr ((-$len/2)+($i-1)*$dim2)] [expr ($h1+$count*$dim2)]

}
};


#build node layer3

set count 0;
for {set j 1} {$j<=[expr ($h3/$dim3)+1]} {incr j 1} {
set count [expr $j-1]
for {set i 1} {$i<=[expr ((($len/2)-($bs/2))/$dim3)+1]} {incr i 1} {
node [expr int(300000+$i+($j-1)*(((($len/2)-($bs/2))/$dim3)+1))] [expr ((-$len/2)+($i-1)*$dim3)] [expr ($h1+$h2+$count*$dim3)]

}
};


#build node layer5

set count 0;
for {set j 1} {$j<=[expr ($h3/$dim3)+1]} {incr j 1} {
set count [expr $j-1]
for {set i 1} {$i<=[expr ((($len/2)-($bs/2))/$dim3)+1]} {incr i 1} {
node [expr int(500000+$i+($j-1)*(((($len/2)-($bs/2))/$dim3)+1))] [expr (($bs/2)+($i-1)*$dim3)] [expr ($h1+$h2+$count*$dim3)]

}
};



#build node layer4

set count 0;
for {set j 1} {$j<=[expr (($h3-$df-$fw)/$dim3)+1]} {incr j 1} {
set count [expr $j-1]
for {set i 1} {$i<=[expr ($bs/$dim3)+1]} {incr i 1} {
node [expr (int(400000+$i+($j-1)*(($bs/$dim3)+1)))] [expr ((-$bs/2)+($i-1)*$dim3)] [expr ($h1+$h2+$count*$dim3)]

}
};



#build node layer6

set count 0;
for {set j 1} {$j<=[expr ($fw/$dim3)+1]} {incr j 1} {
set count [expr $j-1]
for {set i 1} {$i<=[expr ($bs/$dim3)+1]} {incr i 1} {
node [expr (int(600000+$i+($j-1)*(($bs/$dim3)+1)))] [expr ((-$bs/2)+($i-1)*$dim3)] [expr ($h1+$h2+$h3-$df-$fw+$count*$dim3)]
}
};




#************************************************************************************
# 6. DEFINE SOIL MATERIALS

nDMaterial PressureIndependMultiYield 555 2 $rho $G $bulk $cohesion $gammaPeak \
$phi $refPress $pressCoeff -16\
1.00e-6 1.000 2.00e-6 1.000 5.00e-6 0.996 \
1.00e-5 0.984 2.00e-5 0.975 5.00e-5 0.922 \
1.00e-4 0.850 2.00e-4 0.734 5.00e-4 0.532 \
1.00e-3 0.367 2.00e-3 0.224 5.00e-3 0.139 \
1.00e-2 0.085 2.00e-2 0.051 5.00e-2 0.027 \
1.00e-1 0.021




nDMaterial ElasticIsotropic 666 2e15 0.0


#################################*******************************************

# . DEFINE SOIL ELEMENTS

set wgtY [expr -9.81*$rho]
set wgtX 0.0

#layer 1 element
for {set j 1} {$j<=[expr $h1/$dim1]} {incr j 1} {
for {set i 1} {$i<=[expr ($len/$dim1)]} {incr i 1} {
set nI [expr int(100000+$i+($j-1)*(($len/$dim1)+1))]
set nJ [expr int(100000+1+$i+($j-1)*(($len/$dim1)+1))]
set nK [expr int(100000+($len/$dim1)+2+$i+($j-1)*(($len/$dim1)+1))]
set nL [expr int(100000+($len/$dim1)+1+$i+($j-1)*(($len/$dim1)+1))]

element quad [expr int($i+100000+($j-1)*($len/$dim1))] $nI $nJ $nK $nL 1.0 "PlaneStrain" 555 0.0 0.0 $wgtX $wgtY

}
};



#layer 2 element
for {set j 1} {$j<=[expr $h2/$dim2]} {incr j 1} {
for {set i 1} {$i<=[expr ($len/$dim2)]} {incr i 1} {

set nI [expr int(200000+$i+($j-1)*(($len/$dim2)+1))]
set nJ [expr int(200000+1+$i+($j-1)*(($len/$dim2)+1))]
set nK [expr int(200000+($len/$dim2)+2+$i+($j-1)*(($len/$dim2)+1))]
set nL [expr int(200000+($len/$dim2)+1+$i+($j-1)*(($len/$dim2)+1))]

element quad [expr int($i+200000+($j-1)*($len/$dim2))] $nI $nJ $nK $nL 1.0 "PlaneStrain" 555 0.0 0.0 $wgtX $wgtY


}
};


#layer 3 element
for {set j 1} {$j<=[expr $h3/$dim3]} {incr j 1} {
for {set i 1} {$i<=[expr ((($len/2)-($bs/2))/$dim3)]} {incr i 1} {

set nI [expr int(300000+$i+($j-1)*(((($len/2)-($bs/2))/$dim3)+1))]
set nJ [expr int(300000+1+$i+($j-1)*(((($len/2)-($bs/2))/$dim3)+1))]
set nK [expr int(300000+((($len/2)-($bs/2))/$dim3)+2+$i+($j-1)*(((($len/2)-($bs/2))/$dim3)+1))]
set nL [expr int(300000+((($len/2)-($bs/2))/$dim3)+1+$i+($j-1)*(((($len/2)-($bs/2))/$dim3)+1))]

element quad [expr int($i+300000+($j-1)*((($len/2)-($bs/2))/$dim3))] $nI $nJ $nK $nL 1.0 "PlaneStrain" 555 0.0 0.0 $wgtX $wgtY

}
};



#layer 5 element
for {set j 1} {$j<=[expr $h3/$dim3]} {incr j 1} {
for {set i 1} {$i<=[expr ((($len/2)-($bs/2))/$dim3)]} {incr i 1} {

set nI [expr int(500000+$i+($j-1)*(((($len/2)-($bs/2))/$dim3)+1))]
set nJ [expr int(500000+1+$i+($j-1)*(((($len/2)-($bs/2))/$dim3)+1))]
set nK [expr int(500000+((($len/2)-($bs/2))/$dim3)+2+$i+($j-1)*(((($len/2)-($bs/2))/$dim3)+1))]
set nL [expr int(500000+((($len/2)-($bs/2))/$dim3)+1+$i+($j-1)*(((($len/2)-($bs/2))/$dim3)+1))]

element quad [expr int($i+500000+($j-1)*((($len/2)-($bs/2))/$dim3))] $nI $nJ $nK $nL 1.0 "PlaneStrain" 555 0.0 0.0 $wgtX $wgtY


}
};




#layer 4 element
for {set j 1} {$j<=[expr ($h3-$df-$fw)/$dim3]} {incr j 1} {
for {set i 1} {$i<=[expr (($bs)/$dim3)]} {incr i 1} {
set nI [expr int(400000+$i+($j-1)*($bs/$dim3+1))]
set nJ [expr int(400000+1+$i+($j-1)*($bs/$dim3+1))]
set nK [expr int(400000+(($bs)/$dim3)+2+$i+($j-1)*($bs/$dim3+1))]
set nL [expr int(400000+(($bs)/$dim3)+1+$i+($j-1)*($bs/$dim3+1))]

element quad [expr int($i+400000+($j-1)*(($bs)/$dim3))] $nI $nJ $nK $nL 1.0 "PlaneStrain" 555 0.0 0.0 $wgtX $wgtY

}
};



#layer 6 element
for {set j 1} {$j<=[expr ($fw)/$dim3]} {incr j 1} {
for {set i 1} {$i<=[expr (($bs)/$dim3)]} {incr i 1} {
set nI [expr int(600000+$i+($j-1)*($bs/$dim3+1))]
set nJ [expr int(600000+1+$i+($j-1)*($bs/$dim3+1))]
set nK [expr int(600000+(($bs)/$dim3)+2+$i+($j-1)*($bs/$dim3+1))]
set nL [expr int(600000+(($bs)/$dim3)+1+$i+($j-1)*($bs/$dim3+1))]

element quad [expr int($i+600000+($j-1)*(($bs)/$dim3))] $nI $nJ $nK $nL 1.0 "PlaneStrain" 666 0.0 0.0 $wgtX $wgtY


}
};




#### define equal DOF for same node

##bound layer1 and layer2
set fpel [expr ($len/$dim1+1)*($h1/$dim1+1)-($len/$dim1)] ;
for {set k 1} {$k<=(($len/$dim1)+1)} {incr k 1} {
equalDOF [expr int(100000+$fpel-1+$k)] [expr int(200000+(2*$k)-1)] 1 2

};

##bound layer2 and layer3
set fpel [expr ($len/$dim2+1)*($h2/$dim2+1)-($len/$dim2)] ;
for {set k 1} {$k<=(((($len/2)-($bs/2))/$dim2)+1)} {incr k 1} {
equalDOF [expr int(200000+$fpel-1+$k)] [expr int(300000+(2*$k)-1)] 1 2

};

##bound layer2 and layer5
set fpel [expr ($len/$dim2+1)*($h2/$dim2+1)-(($len/2-$bs/2)/$dim2)];
for {set k 1} {$k<=(((($len/2)-($bs/2))/$dim2)+1)} {incr k 1} {
equalDOF [expr int(200000+$fpel-1+$k)] [expr int(500000+(2*$k)-1)] 1 2

};


##bound layer2 and layer4
set fpel [expr ($len/$dim2+1)*($h2/$dim2+1)-(($len-($len/2-$bs/2))/$dim2)];
for {set k 1} {$k<=($bs/$dim2+1)} {incr k 1} {
equalDOF [expr int(400000+(2*$k)-1)] [expr int(200000+$fpel-1+$k)] 1 2

};


##bound layer3 and layer4
set fpel [expr (($len/2-$bs/2)/$dim3)+1];
for {set k 1} {$k<=((($h3-$df-$fw)/$dim3)+1)} {incr k 1} {
equalDOF [expr int(400001+($k-1)*(($bs/$dim3)+1))] [expr int(300000+$fpel+($k-1)*((($len/2-$bs/2)/$dim3)+1))] 1 2
};



##bound layer4 and layer5
set fpel [expr $bs/$dim3];
for {set k 1} {$k<=((($h3-$df-$fw)/$dim3)+1)} {incr k 1} {
equalDOF [expr int(400001+$fpel+($k-1)*(($bs/$dim3)+1))] [expr int(500001+($k-1)*((($len/2-$bs/2)/$dim3)+1))] 1 2
};

##bound layer4 and layer6
set fpel [expr ($bs/$dim3+1)*(($h3-$df-$fw)/$dim3+1)-($bs/$dim3)];
for {set k 1} {$k<=($bs/$dim3+1)} {incr k 1} {
equalDOF [expr int(400000+$fpel+($k-1))] [expr int(600000+$k)] 1 2

};




##bound layer3 and layer6
set fpel [expr (($len/2-$bs/2)/$dim3+1)*(($h3-$df-$fw)/$dim3+1)];
for {set k 1} {$k<=($fw/$dim3+1)} {incr k 1} {
equalDOF [expr int(300000+$fpel+($k-1)*(($len/2-$bs/2)/$dim3+1))] [expr int(600001+($k-1)*($bs/$dim3+1))] 1 2

};


##bound layer5 and layer6
set fpel [expr (($len/2-$bs/2)/$dim3+1)*(($h3-$df-$fw)/$dim3+1)-(($len/2-$bs/2)/$dim3)];
for {set k 1} {$k<=($fw/$dim3+1)} {incr k 1} {
equalDOF [expr int(500000+$fpel+($k-1)*((($len/2)-($bs/2))/$dim3+1))] [expr int(600001+($bs/$dim3)+($k-1)*($bs/$dim3+1))] 1 2
};







puts "Finished creating all -ndf 2 nodes..."

#-----------------------------------------------------------------------------------------


#--------------------------************************------------------------------------
# CREATE structure NODES AND FIXITIES
#--------------------------************************----------------------------------------
puts "CREATE structure NODES AND FIXITIES..."


model BasicBuilder -ndm 2 -ndf 3

source Wsection.tcl; # procedure to define fiber W section
source BoxSection.tcl;
#####under groud node

for {set j 1} {$j<=($bs/$bay)+1} {incr j 1} {
for {set i 1} {$i<=($df/$dim3)+1} {incr i 1} {
node [expr int(700000+10000*$j+$i)] [expr (-$bs)/2+($j-1)*($bay)] [expr $h1+$h2+$h3-$df+($i-1)*$dim3]
}
};


##bound layer7 and layer6

set fpel [expr ($bs/$dim3+1)*(($fw)/$dim3+1)-($bs/$dim3)];
for {set k 1} {$k<=($bs/$bay+1)} {incr k 1} {
equalDOF [expr int(700001+$k*10000)] [expr int(600000+$fpel+($k-1)*($bay/$dim3))] 1 2

};




##### up the groud node

for {set j 1} {$j<=($bs/$bay)+1} {incr j 1} {
for {set i [expr int(($df/$dim3)+2)]} {$i<= [expr int(($df/$dim3)+1+$numst)] } {incr i 1} {
node [expr int(700000+10000*$j+$i)] [expr (-$bs)/2+($j-1)*($bay)] [expr $h1+$h2+$h3+($i-($df/$dim3+1))*$h4]

}
};





# define MATERIAL properties ----------------------------------------
set uFy 24e6
set uE 2e10
set ub 0.001
set uR0 15
set ucR1 0.925
set ucR2 0.15

uniaxialMaterial Steel02 1 $uFy $uE $ub $uR0 $ucR1 $ucR2


# ELEMENT properties
# Structural-Steel properties
# column sections: box250x25
set D 0.25
set B 0.25
set T 0.025
set nfTBl 16
set nfTBt 4
set nfSl 16
set nfSt 4
BoxSection 1 1 $D $B $T $nfTBl $nfTBt $nfSl $nfSt


# column sections: box200x20
set D 0.2
set B 0.2
set T 0.02
set nfTBl 16
set nfTBt 4
set nfSl 16
set nfSt 4
BoxSection 2 1 $D $B $T $nfTBl $nfTBt $nfSl $nfSt

# column sections: box200x15
set D 0.2
set B 0.2
set T 0.015
set nfTBl 16
set nfTBt 4
set nfSl 16
set nfSt 4
BoxSection 3 1 $D $B $T $nfTBl $nfTBt $nfSl $nfSt



# beam sections: W300x8-F250x15
set d 0.315 ; # depth
set bf 0.25 ; # flange width
set tf 0.015 ; # flange thickness
set tw 0.008 ; # web thickness
set nfdw 16; # number of fibers along dw
set nftw 4; # number of fibers along tw
set nfbf 16; # number of fibers along bf
set nftf 4; # number of fibers along tf
Wsection 4 1 $d $bf $tf $tw $nfdw $nftw $nfbf $nftf





geomTransf Linear 1

########### element substructures

## 1.column

for {set j 1} {$j<=($bs/$bay)+1} {incr j 1} {
for {set i 1} {$i<=[expr int($df/$dim3)]} {incr i 1} {
element nonlinearBeamColumn [expr int(710000+($j-1)*($df/$dim3)+$i)] [expr 700000+$j*10000+$i] [expr 700000+$j*10000+$i+1] 5 1 1
}
};


## 2.beam
for {set j 1} {$j<=$bs/$bay} {incr j 1} {

for {set i 1} {$i<=$underf} {incr i 1} {
set node1 [expr int(700000+$j*10000+($i*$underh)/$dim3+1)]
set node2 [expr int(700000+($j+1)*10000+($i*$underh)/$dim3+1)]
element nonlinearBeamColumn [expr int(720000+($j-1)*$underf+$i)] $node1 $node2 5 4 1
}
};



########### element for superstructure

#### 1.build column

for {set j 1} {$j<=($bs/$bay+1)} {incr j 1} {
for {set i 1} {$i<=($numst)} {incr i 1} {

if {$i<=$nssh} {
element nonlinearBeamColumn [expr int(810000+$i+$numst*($j-1))] [expr int(700000+10000*$j+$df/$dim3+$i)] [expr int(700000+10000*$j+$df/$dim3+1+$i)] 5 1 1
}
if {$nssh<$i && $i<=2*$nssh} {
element nonlinearBeamColumn [expr int(810000+$i+$numst*($j-1))] [expr int(700000+10000*$j+$df/$dim3+$i)] [expr int(700000+10000*$j+$df/$dim3+1+$i)] 5 2 1
}

if {2*$nssh<$i && $i<=3*$nssh} {
element nonlinearBeamColumn [expr int(810000+$i+$numst*($j-1))] [expr int(700000+10000*$j+$df/$dim3+$i)] [expr int(700000+10000*$j+$df/$dim3+1+$i)] 5 3 1
}

}
};


## 2.build beam

for {set j 1} {$j<=($bs/$bay)} {incr j 1} {
for {set i 1} {$i<=($numst)} {incr i 1} {

if {$i<=$nssh} {
set node1 [expr int(700000+10000*$j+$df/$dim3+$i+1)]
set node2 [expr int(700000+10000*($j+1)+$df/$dim3+$i+1)]
element nonlinearBeamColumn [expr int(820000+$i+$numst*($j-1))] $node1 $node2 5 4 1
}
if {$nssh<$i && $i<=2*$nssh} {
set node1 [expr int(700000+10000*$j+$df/$dim3+$i+1)]
set node2 [expr int(700000+10000*($j+1)+$df/$dim3+$i+1)]
element nonlinearBeamColumn [expr int(820000+$i+$numst*($j-1))] $node1 $node2 5 4 1
}

if {2*$nssh<$i && $i<=3*$nssh} {
set node1 [expr int(700000+10000*$j+$df/$dim3+$i+1)]
set node2 [expr int(700000+10000*($j+1)+$df/$dim3+$i+1)]
element nonlinearBeamColumn [expr int(820000+$i+$numst*($j-1))] $node1 $node2 5 4 1
}
}
};



#-----------------------------------------------------------------------------------------
# CREATE CONTACT MATERIAL FOR BEAM CONTACT ELEMENTS
#-----------------------------------------------------------------------------------------

# two-dimensional contact material
nDMaterial ContactMaterial2D 777 0.1 1000.0 0.0 0.0

puts "Finished creating all contact materials..."

#-----------------------------------------------------------------------------------------
# CREATE BEAM CONTACT ELEMENTS
#-----------------------------------------------------------------------------------------


# set gap and force tolerances for beam contact elements
set gapTol 1.0e-8
set forceTol 1.0e-8

# define beam contact elements


set fpel22 [expr (($len/2-$bs/2)/$dim3+1)*(($h3-$df)/$dim3+1)];
set fpel33 [expr (($len/2-$bs/2)/$dim3+1)*(($h3-$df)/$dim3)+1];



for {set i 1} {$i<=($df/$dim3)} {incr i 1} {

set iNode1 [expr 710000+$i]
set jNode1 [expr 710000+$i+1]
set sNode1 [expr int(300000+$fpel22+($i-1)*((($len/2-$bs/2)/$dim3+1)))]
element BeamContact2D [expr $i+830000] $iNode1 $jNode1 $sNode1 [expr 900000+$i] 777 0.125 $gapTol $forceTol

set iNode2 [expr int(700000+10000*($bs/$bay+1)+$i)]
set jNode2 [expr int(700000+10000*($bs/$bay+1)+$i+1)]
set sNode2 [expr int(500000+$fpel33+($i-1)*((($len/2-$bs/2)/$dim3+1)))]
element BeamContact2D [expr $i+840000] $iNode2 $jNode2 $sNode2 [expr int(900000+$df/$dim3+$i)] 777 0.125 $gapTol $forceTol
};







puts "***************>>>DONE"