Problem when using central difference method

[buddhasarah]

Hello, I tried to analyze an elastic cantilever beam by using central difference method but got very weird results (deformation increases to unreasonably large level). But I can get good results if newmark method was applied using the same code. THANK YOU VERY MUCH in advance if you can help me.

# ---- source code ------------------------------------------------------------------

wipe; # clear memory of all past model definitions
file mkdir Data; # create data directory
model BasicBuilder -ndm 2 -ndf 3; # Define the model builder, ndm=#dimension, ndf=#dofs


# define GEOMETRY -------------------------------------------------------------
set LCol 3000; # column length
set Weight 6000000000.; # superstructure weight

# calculated parameters
set PCol $Weight; # nodal dead-load weight per column
set g 9806; # g.
set Mass [expr $PCol/$g]; # nodal mass


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

# Single point constraints -- Boundary Conditions
fix 1 1 1 1; # node DX DY RZ

# nodal masses:
mass 2 $Mass $Mass 1e-9; # node#, Mx My Mz, Mass=Weight/g, neglect rotational inertia at nodes

# Define ELEMENTS & SECTIONS -------------------------------------------------------------

# MATERIAL parameters

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

element elasticBeamColumn 1 1 2 250000 30000 62500000000 $ColTransfTag;

recorder Drift -file Drift.txt -time -iNode 1 -jNode 2 -dof 1 -perpDirn 2 ; # lateral drift

# define GRAVITY -------------------------------------------------------------
pattern Plain 1 Linear {
#load 2 0 -$PCol 0
load 2 [expr 0.1*$PCol] 0 0
}
system ProfileSPD
integrator LoadControl 1.0 1 1.0 10.0
test EnergyIncr 1.0e-12 10 0
algorithm Newton
numberer RCM
constraints Plain
analysis Static
# Perform the analysis
analyze 10



# ---------------------------------------
# Create and Perform the dynamic analysis
# ---------------------------------------

# Remove the static analysis & reset the time to 0.0
wipeAnalysis
setTime 0.0

# Now remove the loads and let the beam vibrate
remove loadPattern 1

# Create the transient analysis

algorithm Linear
numberer RCM
constraints Plain
integrator CentralDifferenceNoDamping
system Diagonal
analysis Transient
analyze 10000 0.025

# ------ Part of results -----------------------------------------------------------------------------------
1 0.96
2 1.92
3 2.88
4 3.84
5 4.8
6 5.76
7 6.72
8 7.68
9 8.64
10 9.6
0.025 9.59796
0.05 9.07433
0.075 8.1178e+017
0.1 -1.26841e+036
0.125 1.98189e+054
0.15 -3.0967e+072
0.175 4.83859e+090
0.2 -7.56029e+108
0.225 1.1813e+127
0.25 -1.84577e+145
0.275 2.88402e+163
0.3 -4.50628e+181
0.325 7.04107e+199
0.35 -1.10017e+218
0.375 1.71901e+236
0.4 -2.68596e+254
0.425 4.1968e+272
0.45 -6.55751e+290
0.475 1.#INF
0.5 -1.#IND
0.525 -1.#IND
0.55 -1.#IND
0.575 -1.#IND
0.6 -1.#IND
0.625 -1.#IND
0.65 -1.#IND
0.675 -1.#IND
0.7 -1.#IND
0.725 -1.#IND
0.75 -1.#IND
0.775 -1.#IND
0.8 -1.#IND
0.825 -1.#IND
0.85 -1.#IND
0.875 -1.#IND
0.9 -1.#IND
0.925 -1.#IND
0.95 -1.#IND
0.975 -1.#IND

[fmk]

CentralDifference does have a stability condition, i.e. dT must be smaller than a certain value otherwise method is unstable and stupidly large response is predicted.

[buddhasarah]

Thank you very much Dr. Frank~
The problem has been solved.
Actually I first tried to reduce the time step to an extremely small degree. However, convergence cannot be achieved until I increase the nodal mass in rotational direction to a certain degree.
Although the convergence problem of the typical example presented above has been solved, the result of some models in which large number of DOF exist always disperses no matter how small time step is.
Can you give me some suggestions for this?

Thank you very much!

[fmk]

you do realize that with CentralDifference with no damping you only have some factor times the mass matrix on the left hand side of the equation and so if no man at a rotational dof (or even a very small one) the matrix equation is either cannot be solved or is ill-conditioned .. you might change the solver to see error messages.