Einführung in PFC

Simulation der Partikeldynamik mittels der Diskrete Elemente Methode:
Einführung in PFC
Programmversionen
 Particle Flow Code in 2 Dimensions (PFC2D, Version 4.0)  begrenzt auf 1.000
Partikel
 Particle Flow Code in 3 Dimensions (PFC3D, Version 4.0)  begrenzt auf 5.000
Partikel
6 Bedienungshandbücher in Englisch
Programmebenen
 Kommando-Modus mit PFC-Kommandos in Befehlszeile (Abb. 1)
Abb. 1: PFC – Arbeitsbereich und Plot
 Batch-Modus mit PFC-Kommandos in Editordatei
Beispiel:
new wall id 4 face (0,10,0) (10,10,0) (10,10,10) (0,10,10) generiert plot wall=black ;und angezeigt plot show ;eine Wand wird  Fish-Programmiersprache mit Fish- und PFC-Kommandos
1
Beispiel:
new define make_wall_and_calculate hh=22 ;Zuweisungsoperation abc=hh*3+5 ;Rechnung command wall id 4 face (0,10,0) (10,10,0) (10,10,10) (0,10,10) ;PFC‐Kommando endcommand end make_wall_and_calculate plot wall=black plot show  C++-Programmiersprache (optional, verfügbar in Vollversion)
Nützliche Tastaturbefehle:
F3
Wiederholung der letzten Eingabezeile
Esc
abbrechen des Programms
Bei aktivem Plot:
Strg + C
Richtungstasten
M
Strg + R
X
Y
Z
Ploteinstellungen
Bewegung des Plots
heranzoomen
automatische Zentrierung des Plots mit Rotation (0,0,0)
Rotation um die x-Achse
Rotation um die y-Achse
Rotation um die z-Achse
Verringerung der Rechenzeit:
Partikelgröße
keine kleinen Partikel (monodisperse Partikel)
Steifigkeit
keine hohen Steifigkeiten
Ausschalten der Plots
set pinterval xyz
Plotaktualisierung (aller xyz Zeitschritte) anpassen/verringern
Differential Density Scaling Erklärung weiter hinten im Dokument
Grundlegende PFC-Kommandos
Nützliche Befehle:
; - Semikolon
& - und
pause key
Kommentar
Fortsetzung der auf 80 Zeichen begrenzten Eingabezeile auf der
nächsten Zeile
unterbricht das Programm (Any Key to continue)
new
leert den Speicher, damit ein neues Programm gestartet werden kann
Alle Befehle weisen dieselbe Struktur auf. Sie bestehen aus einem primären Befehlswort
gefolgt von einem oder mehreren Keywords und einem numerischen Input.
COMMAND keyword value . . . <keyword value . . . >
2
Die Befehle können ausgeschrieben werden, die fettgedruckten Buchstaben müssen aber
mindestens angegeben werden.
Befehle, Keywords und numerische Werte können durch Leerzeichen, Komma(s), Runde
Klammern ( ) oder ein Gleichheitszeichen = getrennt werden. Groß- und Kleinbuchstaben
sind erlaubt.
Beispiel:
new wall id=1 kn=1e8 ks=1e8 friction=0.1 face= (0,0,0) (10,0,0) (10,10,0) (0,10,0) ;bottom wall wall id=2 kn=1e8 ks=1e8 friction=0.1 face= 0,0,0 0,10,0 0,10,10 0,0,10 ;left wall wall id=3 kn=1e8 ks=1e8 friction=0.1 face= 10 0 0 10 0 10 10 10 10 10 10 0 ;right wall wall id 4 kn 1e8 ks 1e8 friction 0.1 face 0 0 0 0 0 10 10 0 10 10 0 0 ;front wall WALL id=5 KN=1e8 KS=1e8 FRICTION=0.1 FACE= (0,10,0) (10,10,0) (10,10,10) (0,10,10) ;back wall w id=6 kn=1e8 ks=1e8 f=0.1 fa= (0,10,10) (10,10,10) (10,0,10) (0,0,10) ;upper wall pause key plot create schwarze_wand ;Plot definieren plot add wall black ;Wände zum Plot hinzufügen plot show ;Plot anzeigen Befehlswort
Keywords
wall
id
kn
ks
friction
face
…
Numerische Eingabewerte
Wand
Wand
Identifikationsnummer
Normalsteifigkeit
Schersteifigkeit
Reibungskoeffizient
bei einer begrenzten (finiten) Wand folgen daraufhin die
Koordination
1e8, 0, 10
aktive Seite  Korkenzieherregel / Rechte-Faust-Regel
3
GENERATE und BALL Befehl
Ball-Befehl
Generate-Befehl
Überlappung
Tabelle 1: Beschreibung des Ball- und Generate-Befehls
Ball-Befehl
Generate-Befehl
ball id=26 radius=0.5 x=5 y=2 z=1
generate id=1,25 radius=0.2,1
ball id=27 radius=1.0 x=4 y=1.5 z=1.1 x=0,10 y=-5,-0.5 z=0,10
genaue x-y-z-Platzierung der Partikel
zufällige Anordnung der Partikel
Schlüsselwort radius bestimmt den
die Partikelradien sind zufällig
Partikelradius
gleichverteilt oder durch das
Schlüsselwort gauss normalverteilt
Partikel können sich überlappen,
Partikel überlappen sich nicht
abhängig von ihren Koordinaten
durch Überlappungen können
keine repulsiven Kräfte
repulsive Kräfte auftreten
Beispiel:
new generate id=1,25 radius=0.2,1 x=0,10 y=0.5,5 z=0,5 ;Anzahl und Radius mehrerer Partikel, die in einem quaderförmigen Bereich zufällig verteilt generiert werden ball id=26 radius=0.5 x=‐3 y=1.5 z=1 ;einzelnes Partikel mit genauen Koordinaten generieren ball id=27 radius=1.0 x=‐2 y=1.5 z=1 plot create schwarze_wand ;Plot definieren plot add wall black ;Wände anzeigen plot add ball range id 1 25 blue ;Partikel mit id 1 bis 25 anzeigen plot add ball range id 26 27 red ;Partikel 26 und 27 anzeigen plot add axes black ;Achsen anzeigen plot set background white ;Hintergrund weiß plot show ;anzeigen des Plots 4
property
change
initialize
Definition der Eigenschaften existierender Partikel
Synonym zu property
Synonym zu property
Stoffeigenschaften für die Partikeln und ihre Verbindungen
property density=2650 kn=1e8 ks=1e8
Dichte, Normal- und Schersteifigkeit
property radius multiply=1.51
Vergrößerung des Radius
Die Anzahl, die Koordinaten, die Radien, die Dichten und die Steifigkeiten der Partikel, die
Koordinaten und die Steifigkeiten der Wände und die Erdbeschleunigung werden festgelegt.
Daraus ergeben sich die wirkenden Gewichts- und Abstoßkräfte, die Beschleunigungen und
Spannungen in den Überlappungen.
Beispiel:
new plot create schwarze_wand plot add wall black plot add ball range id 1 25 blue plot add ball range id 26 27 red plot add axes black plot set background white plot show generate id=1,25 radius=0.2,1 x=0,10 y=0.5,5 z=0,5 ball id=26 radius=0.5 x=‐3 y=1.5 z=1 ball id=27 radius=1.0 x=‐2 y=1.5 z=1 property density=2650 kn=1e8 ks=1e8 ;zählt für alle existierenden Partikel pause key range name=bälle_1_25 id=1,25 ;Definition des Bereichs der Partikel mit id 1 bis 25 property radius multiply=1.51 range=bälle_1_25 ;Zuweisung der Vergrößerung des Radius für den definierten Bereich pause key set gravity=(0,0,‐9.81) ;Gravitationsbeschleunigung cycle 6000 ;Zeitschritte Tabelle 2: Keywords in den Befehlen
Stoffdaten
Randbedingungen:
 wall-Befehl für Wände
 wall-Befehl für Wände
 property-Befehl für Partikel  property-Befehl für Partikel
Identifikations radius
 x-, y-, z-Koordinaten
nummer:
 kn – Normalsteifigkeit
 xdisplacement, ydisplacement,
 ball – Partikel
zdisplacement – x-, y-, z-Weg
 ks – Schersteifigkeit
 clump – Clump
 xvelocity, yvelocity, zvelocity –
 density – Dichte
x-, y-, z-Geschwindigkeit
 history –
 friction – Reibung
Aufzeichnen von  damping – viskose
 xspin, yspin, zspin – x-, y-, zMesswerten
Rotationsgeschwindigkeit
Dämpfung
 measure –
 gravity – Erdbeschleunigung
Kontaktbindung
Messkreis
 n_bond – Normalsteifigkeit  xforce, yforce, zforce – x-, y-,
 wall – Wand
5
programminterne
Parameter
 s_bond – Schersteifigkeit
Parallelbindung
 pb_kn – Normalsteifigkeit
 pb_ks – Schersteifigkeit
 pb_nstrength –
Normalfestigkeit
 pb_sstrengt –
Scherfestigkeit
 pb_radius – Radius
 etc.
Ausgabefenster
 position –
Position
 size – Größe
 Farbe
 etc.


z-Kraft
xmoment, ymoment, zmoment
– x-, y-, z-Drehmoment
etc.
Tabelle 3: Definition der Eingabe- und Ausgabewerte
Eingabewert
Länge
m
Dichte
kg/m3
Gravitationsbeschleunigung
m/s2
Steifigkeit der Wände und Partikel N/m
Steifigkeit der Parallelbindung
Pa/m
N/(m2m)
Festigkeit von Kontaktbindungen
N
Festigkeit von Parallelbindungen
N/m2
Ausgabewert Kraft
N
Spannung
Pa
Zerlegung der Kontaktkraft in eine Normal- und eine Scherkomponente
F
Fn
F
kn
Fs
kn
Kontaktfläche
ks
F – Kontaktkraft
Fn – Normalkraft
Fs – Scherkraft
kn – Normalsteifigkeit
ks – Schersteifigkeit
6
Kontaktkräfte dargestellt als freie ungedämpfte Federschwingung
Freier ungedämpfter
Federschwinger
Masse bringt Zug- und
Druckkräfte auf die Feder auf
Lineares Kontaktgesetz
Masse bringt nur eine Druckkraft auf die Feder auf.
F
F
k
Un
kn
x
m
ks
Us
m
m
x – Weg/Auslenkung
k – Federkonstante
m – Masse
Un, Us – Überlappung in Normal- und Scherrichtung
kn, ks – Steifigkeit in Normal- und Scherrichtung
Wirkende Federkraft
F  k  x
Druck
F
Normalkraft
Scherkraft
Fn = -k n  U n
Fn
Fs = -k s  U s
Fs
ks
xm
x
-xm
Un
kn
Us
Druck
Zug
Reihen- und Parallelschaltung
zweier Federn
Kontakt zweier Partikel entspricht der Reihenschaltung
zweier Federn, welche zusammengepresst werden / sich
überlappen.
Berechnung der Normalsteifigkeit:
k1
k2
k1
m
k2
k1
k ges  k 1  k 2
m
k2
1
1
1
 
k ges k1 k 2
1
1
1
 
k ges k1 k 2 Berechnung der Schersteifigkeit ist analog.
Steifigkeit (kn, ks) und Eigenkreisfrequenz ω0 für ein Kontinuum (Wände und Partikel)
Vergleich zwischen
7
mechanischem Schwinger
Ruhelage
schwingendem Kontinuum
Auslenkung
Auslenkung
Ruhelage
lo
k
m
A0
E
l1
m
FF
x(t)
lo(t)
FS
m
m
FT
FT
Federschwinger
Auslenkung:
x(t)
Hookesches Gesetz für
Federschwinger:
FF = -kx
FF - Federkraft
k - Federkonstante
x - Weg
Federsteifigkeit:
k = m  ω02

Trägheitskraft: FT = ma = mx
Bewegungsgleichung:
FT = FF
 = -kx
mx
k

x+ x=0
m
Eigenkreisfrequenz:
k
ω0 = 2πf 0 =
m
Stabschwinger
Auslenkung:
l( t )  l 0 ( t )  l 0
l1 ( t )  l 0
l0
 - Dehnung
Hookesches Gesetz für das Kontinuum
l
  E  E
l0
EA 0
FS = -E  A 0  ε = Δl
l0
E - Elastizitätsmodul
A0 - Querschnittsfläche
l0 - Ausgangslänge
Kontinuum Steifigkeit:
EA 0
k=
= m  ω02
l0
 = m  Δl = ml ε
Trägheitskraft: F = mx
T
0
Bewegungsgleichung:
FT = FS
ml0ε = -EA 0 ε
EA 0
ε +
ε=0
ml0
Eigenkreisfrequenz:
EA 0
0 
ml0
8
Kontaktaufgabe der Elastizitätstheorie, erstmals gelöst von H. Hertz 1881
L. D. Landau, E. M. Lifschitz: Lehrbuch der theoretischen Physik, Band VII,
Elastizitätstheorie
Akademie-Verlag Berlin 1989, S. 37-38
Beispiel: Kontakt zweier gleich großer elastischer Kugeln
F - wirkende Kraft
R - Radius der Kugel
d - Durchmesser = 2R
h - Überlappung
E - Elastizitätsmodul
 - Poissonzahl
k - Federsteifigkeit
F
h
F
1/ 3
 2  D2 

h  F 
 R 
F
F
k (E, , F, R )  
2 1/ 3
h
2/3 2  D 

F 
 R 
3 1 2
D 
2 E
F
1/3
 2  D2 
F 

 R 
k = a F d
2/3
1/3 1/3
 F  E2  d 
=
2 2
 9(1- ν ) 
 E 
a=
2 
 3(1- ν ) 
1/3
1/3
 F  E2 
=
 d1/3
2 2
 9(1- ν ) 
2/3
eine Konstante des Materials in N2/3/mm4/3
Würfel
Kugel
k = lE
k = a  F1/3d1/3
Basalt (E=90 kN/mm2, v=0.2)
80
Steifigkeit k in 10^5 N/m
F
k=
=
Δl
F  k h
2/3
70
60
50
40
30
20
9
10
5
0
1 2 3
4 5 6
7 8 9
10
Kugelradius r in mm
1
Kontaktkraft F
in N
9
Federsteifigkeit zweier gleicher elastischer Kugeln nach Hertz (1881) in
Abhängigk. von Material, Kontaktfläche und -kraft (Würfel im Vergleich)
1e+09
Federsteifigkeit in N/m
1e+08
1e+07
1e+06
1e+05
1e+04
1e+03
1e+02
0,001
0,01
0,1
1
10
100
Kugeldurchmesser/Würfelkantenlänge in mm
1000
Silikonkautschuk F=1mN
Polystyrol F=1mN
Kalkstein/Quarz F=1mN
Basalt F=1mN
Baustahl F=1mN
Silikonkautschuk F=1N
Polystyrol F=1N
Kalkstein/Quarz F=1N
Basalt F=1N
Baustahl F=1N
Silikonkautschuk F=1kN
Polystyrol F=1kN
Kalkstein/Quarz F=1kN
Basalt F=1kN
Baustahl F=1kN
Würfel-Silikonkautschuk
Würfel-Polystyrol
Würfel-Kalkstein/Quarz
Würfel-Basalt
Würfel-Baustahl
10
Randbedingungen
Berechnung in
Zeitschritten
gesetzt durch Wände/Partikel
- einwirkende Kräfte
- fixierte Positionen und Geschw.
Kontaktkräfte
Bewegungsgesetz
Kraft-Weg-Gesetz
angewendet auf Kontakte
Normalkraft:
Fn  k n  U n
Scherkraft : Fs  k s  U s
U.. Überlappung
angewendet auf Partikel
geradlinige Bewegung: Fi  m  x i
Rotationsbewegung:
i
Mi  J  
i..2 für 2D i..3 für 3D
Positionskorrektur von Partikeln,
Wänden und Kontakten
Einsatz der Software PFC3D der Fa. Itasca Co., Minnesota US
11
t=t0, xi=x0,i, vi=v0
Bewegungsgleichungen für jedes Primärpartikel i
Kräfte, Momente
 Kräftebilanz
Positionen,
Geschwindigkeiten
t=t+?t

p
l
 ( ij )  ( ij )
 ( ij )  ( ij )

d 2 ri
mi 2   (FK ,n FK ,s )   (FB,n  FB,s )  mi g
dt
j 1
j 1
 Momentenbilanz

dω i
Ii

dt
FK
FB
l und p
Ii
p
 (ij)  (i)
 (ij)  (ij)
(
r
F
)
M


 K K  B  M mg
l
j 1
j 1
Kontaktkraft Partikel-Partikel oder Partikel-Wand,
Federkraft in der Festkörperbrücke,
Zahl der Kontakte und Festkörperbrücken vom Partikel i
Flächenträgheitsmoment vom Partikel i
Kontaktmodell
Festkörperbrückenmodell
Kontaktnormalkräfte 
 ( ij )


FB(,ijn)  (k B ,ij ,n  Aij  s ij ,n )  n ij
FK ,n  (k ij ,n  s ija,n  ηij  s ij ,n )  n ij
Tangentiale Kontaktkräfte
 ( ij )



FK ,s  min (k ij ,s  s ija,s  ηij  s ij ,s ), μ ij  FK( ij,n)  t ij FB(,ijs)  (k B ,ij ,s  Aij  s ij ,s )  t ij


(ij)
Bruchkriterium
 FB,(ij)n M B
(ij)
(ij)
 ij 
 (ij) RB   max
oder
Aij
IB
kij
kij,B
sij
ij
ij
Kontaktsteifigkeit,
flächenbezogene Brückensteifigkeit,
Überlappung,
Dämpfungsparameter,
Gleitreibungskoeffizient,
 und 
max und max
MB
Aij
IB
 ij 
FB,(ij)s
Aij
(ij)
  max
normale und tangentiale Spannungen
in der Brückenquerfläche,
Bruchfestigkeit und Scherfestigkeit,
Moment in der Brückenquerfläche,
Querschnittsfläche der Festkörperbrücke,
Flächenträgheitsmoment (Aij)
12
Interactions between balls
Contact bond
Fn (tension)
Contact and slip
Normal force
Fn
bond breaks
(tension)
slip
Un (overlap)
n
k
Shear force
Fs
k
Fsmax
Fs
Fsc
slip when
Un>0 AND Fs>Fsmax
slip
contact
ks
Fnc
n
contact
Parallel bond
similar to contact bond with parallel bond
typical pb_kn and pb_ks
[stress/displacement]
Un (overlap)
bond breaks
Fsmax
Us
Entities
Existence
ball-ball and ball-wall
This model is always active, unless a
contact bond is present — in which
case, the contact bond model behavior
supersedes the slip model behavior.
Parameters
friction coefficient µ at the contact
[dimensionless]
PFC command
PROPERTY friction
ks
Us
ball-ball
The existence of a contact bond
precludes the possibility of slip, i.e., the
magnitude of the shear contact force is
not adjusted to remain less than the
allowable maximum of Fsmax
normal contact bond strength Fnc
shear contact bond strength Fsc
[force]
PROPERTY n_bond=x
PROPERTY s_bond=x
ball-ball
The two bond models (contact and parallel)
can occur simultaneously; thus, in the
absence of a contact bond, the slip model is
active in conjunction with the parallel-bond
model.
normal and shear strength [stress]
normal and shear stiffness
[stress/displacement]
bond radius [length]
PROPERTY pb_nstrength=x
PROPERTY pb_sstrength=x
PROPERTY pb_kn=x pb_ks=x
PROPERTY pb_radius=x
13
Installation
Deletion
Contact shape
point
Force/torque
developed at the
contact
Interpretation
Procedure
can only transmit a compression force,
i.e., no tensile force and torque at the
contact.
friction
The criterion of no-normal strength is
enforced by checking whether the
overlap is less than or equal to zero. If
it is, then both the normal and shear
contact forces are set to zero. Then the
contact is checked for slip conditions
by calculating the maximum allowable
shear contact force Fsmax = µ |Fn|
If |Fs| > Fsmax, then slip is allowed to
occur (during the next calculation
cycle) by setting the magnitude of Fs
equal to Fsmax
Contact and parallel bonds are installed at each real contact (with nonzero overlap) and
virtual contact (between 2 balls with separation less than 1e-6 times the mean radius of the
2 balls).
Both types of bonds are deleted if their strengths are exceeded; alternatively, the contact
bonds can be deleted by setting either the n_bond or s_bond values to zero, and the parallel
bonds can be deleted by setting either the pb_nstrength or pb_sstrength values to zero.
If the magnitude of the tensile normal contact force equals or exceeds the normal contact
bond strength, the bond breaks, and both the normal and shear contact forces are set to
zero. If the magnitude of the shear contact force equals or exceeds the shear contact bond
strength, the bond breaks, but the contact forces are not altered, provided that the shear
force does not exceed the friction limit, and provided that the normal force is compressive.
point
either a circular (SET disk off) or rectangular
cross-section (on) lying between the particles
can transmit both a force and torque
can only transmit a force, i.e., no
torque.
adhesion
cementatious material between the two balls
Instead, the magnitude of the shear
contact force is limited by the shear
contact bond strength. Contact bonds
also allow tensile forces to develop at
a contact. These forces arise from the
application of Eq. (1.10) when Un <
0 (i.e., there is no overlap). In this
case, the contact bond acts to bind the
balls together. The magnitude of the
tensile normal contact force is limited
by the normal contact bond strength.
14
Calculation of the cycle t - translational motion example of two attached balls
After ball creation following properties are defined:
- positions of the balls A and B x [i A ] , x [i B] ,
- radii R [ A ] and R [ B] ,
- density,
- normal stiffness K n (secant modulus) and shear stiffness k s (tangent modulus).
i .. 1, 2, 3-index of a vector (x, y, z-direction)
1 Given values at cycle start
x i( t t / 2 ) , x i( t ) , Fi( t ) , m, g i
2 Calculation of the geometry at cycle start
d is the distance between the ball centers:
d  x [i B]  x [i A ]
overlap U n defined as relative contact displacement in the normal direction:
U n  R [ A ]  R [ B]  d
the location of the contact point x [i C ] is given by
x [i C ]  x [i A ]  R [ A ]  12 U n
3 Calculation of the forces at cycle start
contact force vector Fi can be resolved into normal and shear components:
Fi  Fin  Fis
the normal contact force is calculated by
(1)
(2)
(3)
(4)
16
Fin  K n U n
(5)
The new shear contact force is found by summing the old shear force existing at the start of
the time step with the shear force-increment Fis . If the slip model is active the shear force
s
will be corrected below Fmax
when it is above:
Fis  Fis  Fis  F n
(6)
The shear force-increment F will be found by the following equations (7-10):
contact velocity Vi (relative motion at the contact)
Vi  x [i C ] [ B]  x [i C ] [ A ]
(7)
Vi  Vin  Vis the contact velocity can be resolved into normal and shear components
U si  Vis t
shear displacement-increment over a time step of t
(8)
(9)
s
i
Fis  k s U si
shear force-increment
(10)
The contribution of the final contact force (4) to the resultant force on the two balls in contact
is given:
Fi[ A ]  Fi[ A ]  Fi
(11)
[ B]
[ B]
Fi  Fi  Fi
(12)
[A]
[ B]
Fi and Fi are the resultant force, the sum of all externally applied forces acting on the two
particles.
4 Calculation at the middle (velocity) and end (position) of the time interval
x i is calculated at the mid-intervals t  nt / 2 , while x i , x i and Fi are computed at the
primary intervals of t  nt .
Fi  m(x i  g i )
(13)
1 ( t  t / 2)
x i( t ) 
(14)
x i
 x i( t t / 2 )
t
Insert (13) into (14) and solving for the velocities at time ( t  t / 2) result in



 F( t )
x i( t t / 2 )  x i( t t / 2)   i  g i t

 m
Finally, the velocities in (15) are used to update the position of the particle center as
x i( t  t )  x i( t )  x i( t  t / 2) t
(15)
(16)
17
Time steps
The computed solution produced by equations (1-16) will remain stable only if the timestep
does not exceed a critical timestep. A simplified procedure is implemented in PFC2D to
estimate the critical timestep at the start of each cycle. The actual timestep used in any cycle
is taken as a fraction of this estimated critical value. This fraction can be specified using the
SET safety_fac command.
t crit  m / k
(17)
where k is the stiffness of each spring and m the mass of each ball.
Example: t crit  1kg /(108 N / m)  10 4 s
Freie ungedämpfte Federschwingung
FF
FT
k, 0
m
m
0
statische Ruhelage
 kx  mx
k
x0
m
x   02 x  0
Es ist:
0  k / m
FF ... Federkraft
FT ... Trägheitskraft
x
Umlenkpunkt
x 
T  2 m / k
18
Differential Density Scaling
Inertial masses are multiplied by a factor, different for each ball, so that a unit timestep is
obtained. Gravitational masses are unaffected.
Normal unscaled operation
SET dt auto (default)
New different timesteps are calculated
before each cycle step, hence the
calculation time for many particles is
long.
Provides a true dynamic solution that is
valid for both steady- and non-steadystate conditions.
Fi is the resultant force, the sum of all
externally applied forces acting on the
particle:
Fi  m(x i  g i )
The finite-difference expressions for
the particle velocities:

 F( t )
x i( t  t / 2 )  x i( t t / 2)   i  g i t

 m
Differential density scaling
SET dt dscale
Only one timestep ist calculated for a particle
arrangement befor cycling many times and keep
uniform until particle arrangement changes. Hence
the general calculation time is reduced.
If one is interested only in obtaining the steady-state
solution, in which all particle accelerations are zero,
corresponding to either static equilibrium or steady
flow, differential density scaling may be used to
reduce the total number of cycles required to reach
the steady-state condition.
Only the final steady-state solution is valid; transient
states do not represent the true dynamic behaviour
of the system. For path-dependent phenomena that
involve development of mechanisms, this may result
in an unrealistic steady-state solution.
Reduce the total number of cycles required to reach
the steady-state condition.
The inertial mass of each particle is modified at the
start of each cycle such that the stability criterion of
(17) is satisfied with a timestep of unity:
Fi  (m i )x i  (m g )g i
where mi and mg are the inertial and gravitational
masses, respectively; mg is always equal to the
actual mass of the particle and, in the absence of
differential density scaling, mi = mg.
 F( t )  m g  
x i( t  t / 2 )  x i( t t / 2)   i i   i g i t
m m  
19
Mechanical Damping for dissipation of kinetic energy- translational motion example
Local Damping
DAMP (default) local 0.2
keyword 'default' for
newly created balls
(default local =0.7)
Combined Damping
Viscous Damping
DAMP (default) viscous normal 0.2
DAMP (default) viscous shear 0.2
keyword 'default' for newly created
contacts
(default viscous normal/shear=0)
acts on each ball
acts on each contact as a dashpot in
parallel with the existing contact model
A damping-force term is added to the unbalanced A damping force, Dn, Ds is added to the
contact force:
force in the equation of motion:
damp
F( i )  F(i )  mx ( i )
Fin  K n U n  D n
F(i) includes the contribution from the gravity force. Fis  Fis  Fis  D s
 Fi
F(di )   F(i ) sign ( v ( i ) )
Dn  cn Vn
Ds  cs V s
F(li ) 

2
The damping constants are not
d


Fi vi F
specified directly; instead, the critical
 dF( i ) 
  sign ( v ( i ) )
sign 
- + - Mdamping ratios in the normal and shear
 dt 
+


- - + M damping
directions (βn and βs ) are specified:
+ - + Mc n   n c crit
c s  s c scrit
n
+
+ + - M damping
c crit
n  2 m n  2 mk n
c scrit  2ms  2 mk s
ωn and ωs are the natural frequencies of
the undamped system, kn and ks are
the contact tangent stiffnesses
- only accelerating
- used for significant rigid- Viscous damping is characterized by
motion is damped,
body motion of a system in the critical damping ratio β. When β =
- damping is frequency addition to oscillatory
1, the system is said to be critically
independent
motion to be dissipated
damped, meaning that the response
- inappropriate for
- dissipate energy at a
decays to zero at the most rapid rate.
particles in free flight
slower rate compared to
Also, β = 1 represents the transition
under gravity or for
local damping
from an oscillatory response, when β <
impact of particles,
1, to an exponentially decaying
- inappropriate for
response when β > 1. When β < 1, the
systems in which large
system is said to be underdamped, or
groups of particles are
lightly damped, and when β > 1, the
driven by specifiedsystem is said to be overdamped, or
velocity boundary
heavily damped.
conditions.
20
Viscous Damping
 kx  mx  bx
x  2x   02 x  0
mx  bx  kx  0
x
xA
xmax
x A e  t
x(t=0)
/D
0
t
x Ae t cos( D t   )
x A  x max 
x ( t  0)
cos 
Die gebräuchlichsten Dämpfungsmaße sind:
b
Dämpfungskonstante des viskosen Dämpfers
b
Abklingkonstante in der Exponentialfunktion
 
 0 D
2m
b
b
D

dimensionsloses Dämpfungsmaß (Lehrsches Dämpfungsmaß)
bkrit 2 mk
2D

logarithmisches Dämpfungsmaß
1  D2
Schwingfall
0  
D<1
periodisch ablaufende gedämpfte
Schwingung mit Kreisfrequenz
aperiodischer Grenzfall
0  
D=1
aperiodischer Grenzfall:
möglichst in kurzer Zeit in
Ruhelage
 D   02   2
D  0
Allgemeine reelle Lösung des Schwingfalls:
x (t )  e t (C1 sin  D t  C 2 cos  D t )
x (t )  x A e t cos( D t   )
x
x
t
Kriechfall
0  
D>1
keine periodische Schwing.,
kein Durchgang durch die tAchse
D  0
x
t
t
21

Download Report