/ Forside / Karriere / Uddannelse / Højere uddannelser / Nyhedsindlæg
Login
Glemt dit kodeord?
Brugernavn

Kodeord


Reklame
Top 10 brugere
Højere uddannelser
#NavnPoint
Nordsted1 1588
erling_l 1224
ans 1150
dova 895
gert_h 800
molokyle 661
berpox 610
creamygirl 610
3773 570
10  jomfruane 570
Finite element kontra...
Fra : Nick Janson


Dato : 28-03-08 23:16

Jeg har hidtil, med nogenlunde success, lavet for-sjov simuleringer af
varmediffusion, luftstrøminger og deformation af soft bodies med et system
af partikler og fjedre. Det virker egentlig ok og med en fornuftig
integrator kan det være ganske stabilt.
Nu er jeg så begyndt at se på finite element-metoden og minsanten om ikke
jeg er løbet ind i problemer. Der er ganske vidst et utal af introduktioner
af metoden på nettet, men ingen er tilsyneladende introducerende nok. Jeg
vil mene jeg har fat i grundprincipperne, men jeg mangler en del i at kunne
lave min egen løsning, og jeg savner i den grad en step-by-step gennemgang
af et konkret eksempel med så lidt som tre elementer eller der omkring. Den
slags ser jeg ingen steder.

Med fjedre-partikel-modellen ser man på værdierne for partiklerne,
eksempelvis positionen, og finder så fjedrekrafter og integrerer frem i
tiden. Det giver i sig selv en animation hvis man viser alle tidsskridt. Med
finite element virker det som om man finder alle interne og eksterne krafter
(som for partiklerne), opstiller det som Ax=b hvor A er stivhedsmatrixen, x
er partiklernes forskydning og b er de krafter man ved systemet har, og
derefter løser man for x for at finde ligevægtstilstanden. Hvordan skal man
animere dette hvis man går fra en tilstand direkte til ligevægten uden
mellemskridt? Er det overhovedet korrekt hvad jeg siger om FEM/FEA her,
eller har jeg misset mere end jeg troede?

En anden ting jeg er lidt i tvivl om er hvorfor fjeder-partikel-modellen er
FEM underlegen, for det må del velsagtens være siden FEM dominerer som den
gør.


 
 
Nick Janson (29-03-2008)
Kommentar
Fra : Nick Janson


Dato : 29-03-08 13:15

Måske jeg bare skal starte med første usikkerhed som er helt banal.
Hvis man laver en gængs finite elemente analyse på et problem, så finder man
ligevægtstilstanden og det er det? Intet med integration over tid, vel?

Helt simpelt 1D-problem er en fjeder med fjederkonstanten k=100N/m Den har
hvilelængden 4 og nuværende længde 4. I hver ende er en node n0 og n1. n0 er
fikseret i position 0 og n1 er fri i position 4m.
Så påvirkes n1 med en kraft på -200N som søger at trykke den mod n0.

n0------------n1 <=-200N

Opgaven er herefter at finde ud af hvordan dette ender.
Man siger (og det er trivielt for dette eksempel)

-200N = delta n1 * k =>
delta n1 = -200/k =>
delta n1 = -2

Og dette er løsningen for her er der ligevægt?

Det jeg er usikker på er hvordan det nu er løsningen, for hvis en ekstern
kraft påvirkede n1 som beskrevet, så ville n1 flytte sig, men den vil også
få en hastighed som vil gøre at den vil svinge over ligevægtspositionen og
så svinge tilbage og hvis ikke dæmpet, så vil den blive ved med det. Hvad er
det så man finder? Det er tilsyneladende kun ligevægtstilstanden, men hvis
jeg skulle visualisere effekten af denne eksterne kraft så ville jeg jo også
gerne se svingningen, som er helt fraværende.


Ulrik Smed (29-03-2008)
Kommentar
Fra : Ulrik Smed


Dato : 29-03-08 21:46

Nick Janson wrote:
> Måske jeg bare skal starte med første usikkerhed som er helt banal.
> Hvis man laver en gængs finite elemente analyse på et problem, så
> finder man ligevægtstilstanden og det er det? Intet med integration
> over tid, vel?
> Helt simpelt 1D-problem er en fjeder med fjederkonstanten k=100N/m
> Den har hvilelængden 4 og nuværende længde 4. I hver ende er en node
> n0 og n1. n0 er fikseret i position 0 og n1 er fri i position 4m.
> Så påvirkes n1 med en kraft på -200N som søger at trykke den mod n0.
>
> n0------------n1 <=-200N
>
> Opgaven er herefter at finde ud af hvordan dette ender.
> Man siger (og det er trivielt for dette eksempel)
>
> -200N = delta n1 * k =>
> delta n1 = -200/k =>
> delta n1 = -2
>
> Og dette er løsningen for her er der ligevægt?
>
> Det jeg er usikker på er hvordan det nu er løsningen, for hvis en
> ekstern kraft påvirkede n1 som beskrevet, så ville n1 flytte sig, men
> den vil også få en hastighed som vil gøre at den vil svinge over
> ligevægtspositionen og så svinge tilbage og hvis ikke dæmpet, så vil
> den blive ved med det. Hvad er det så man finder? Det er
> tilsyneladende kun ligevægtstilstanden, men hvis jeg skulle
> visualisere effekten af denne eksterne kraft så ville jeg jo også
> gerne se svingningen, som er helt fraværende.

Jeg sidder netop selv og pjatter lidt med noget bølgesimulering efter samme
model med lodder og fjedre. Jeg vil også godt se bølgerne bevæge sig, og
prøver at lave en simulator der kan vise hvordan en impuls ser ud
forskellige steder på en 2D flade (måske senere et 3D rum). Man kan tegne en
scene i Paint og lægge ind, med mure og dæmpning og lydgivere og mikrofoner.
Så jeg ved ikke om man vil kunne bruge denne FEM metode til det, men det
ville sikkert kunne hæve hastigheden væsentligt hvis man kan.

Jeg opdagede lige en sjov ting med en 1D række lodder og fjedre. Når man
sender en puls ind i enden, vandrer der en bølge igennem rækken. Efter
bølgen er der noget 'ripple' som klinger langsomt ud. Det ligner en lavpas
filtrering. Jeg prøvede at sample udsvinget af et bestemt lod lidt inde i
rækken over tid, og køre en FFT på samplen. Det viste et næsten ideelt
lavpasfilter, som når man FFT'er en sinc kurve. Overgangsfrekvensen bliver
dybere jo længere man går hen i rækken og sampler.

--
Ulrik Smed
Aarhus, Denmark



Nick Janson (29-03-2008)
Kommentar
Fra : Nick Janson


Dato : 29-03-08 22:23

> Jeg opdagede lige en sjov ting med en 1D række lodder og fjedre. Når man
> sender en puls ind i enden, vandrer der en bølge igennem rækken. Efter
> bølgen er der noget 'ripple' som klinger langsomt ud. Det ligner en lavpas
> filtrering. Jeg prøvede at sample udsvinget af et bestemt lod lidt inde i
> rækken over tid, og køre en FFT på samplen. Det viste et næsten ideelt
> lavpasfilter, som når man FFT'er en sinc kurve. Overgangsfrekvensen bliver
> dybere jo længere man går hen i rækken og sampler.

Afhængigt af dine masser og fjederkonstanter i den type simulering, der vil
du opleve et low pass filter. Det er fornuftigt nok, og realistisk nok. Hvis
du har en højfrekvens svingning af et punkt, så vil dette punkt gennem
fjedrene påvirke nabopunkterne, men hvis punktet selv svinger tilbage igen
før naboen er kommet op i fart, så vil naboen aldrig få nogen stor
amplitude. Det er mig bekendt sådan dæmpninger fungerer i praksis. Svage
fjedre eller store masser giver begge lave acelerationer F/m = (dx*k)/m går
mod nul for k gående mod nul eller m gående mod uendelig.
Du kan eventuelt lave din bitmaptegning så den indikerer netop masser eller
fjedrekonstanter for dit gitter og dermed indsætte blokke af forskellige
materialer. Så kan du se nogle sjove effekter.
En anden morsom ting er at lade en række punkter svinge som en helhed.
Eksempelvis en række af 100 nabopunkter der svinger sammen og forestiller en
højttalermembran (i en 2d verden er det en linie). Se hvad der sker med
bølgerne ved forskellige frekvenser. Højfrekvens er meget retningsbestemt,
som i virkeligheden, mens lavfrekvens breder sig mere som en kuglebølge.

Hvilken integrator bruger du forøvrigt? Får år siden da jeg først legede med
det, da brugte jeg explicit euler, hvilket naturligvis var simpelt men
elendigt. Senere brugte jeg RK4 (runge kutta 4). til sidst med variabel
skridtstørelse. Det gør at man kan øge hastigheden betragteligt.


Nick Janson (29-03-2008)
Kommentar
Fra : Nick Janson


Dato : 29-03-08 22:27

> En anden morsom ting er at lade en række punkter svinge som en helhed.
> Eksempelvis en række af 100 nabopunkter der svinger sammen og forestiller
> en højttalermembran (i en 2d verden er det en linie). Se hvad der sker med
> bølgerne ved forskellige frekvenser. Højfrekvens er meget retningsbestemt,
> som i virkeligheden, mens lavfrekvens breder sig mere som en kuglebølge.

Du kan se et eksempel på dette her
http://www.greenleaf.dk/tag/projects/wave/index.html
Jeg ved ikke hvilken integrationsmetode der bruges, men det kører da
nogenlunde hurtigt og ser umiddelbart korrekt ud.
Der er en exe man kan downloade og prøve.


Ulrik Smed (29-03-2008)
Kommentar
Fra : Ulrik Smed


Dato : 29-03-08 23:57

Nick Janson wrote:
>> En anden morsom ting er at lade en række punkter svinge som en
>> helhed. Eksempelvis en række af 100 nabopunkter der svinger sammen
>> og forestiller en højttalermembran (i en 2d verden er det en linie).
>> Se hvad der sker med bølgerne ved forskellige frekvenser.
>> Højfrekvens er meget retningsbestemt, som i virkeligheden, mens
>> lavfrekvens breder sig mere som en kuglebølge.
>
> Du kan se et eksempel på dette her
> http://www.greenleaf.dk/tag/projects/wave/index.html
> Jeg ved ikke hvilken integrationsmetode der bruges, men det kører da
> nogenlunde hurtigt og ser umiddelbart korrekt ud.
> Der er en exe man kan downloade og prøve.

Ser fedt ud! Ja, jeg har prøvet at lave en membran med en række
punkter der følges ad. Og tegne horn og refleksporte og alt muligt, og se
resonanser og retningsvirkning og reflektioner osv., det er meget
fascinerende. Men jeg undrede mig over den perfekte lavpasvirkning der er i
'luften' alene, i 1D versionen. Man kan også se 'ripples' i 2D, når man slår
an med en skarp puls.

Jeg ved ikke hvad min integrator hedder, men mit program (C#) har 3 arrays,
et positions (eller pressure) array, og et array for x-hastighed og et for
y-hastighed. Jeg beregner først hastighederne ved at tage den ''gamle'
hastighed og lægger forskellen mellem punktet på hver side til.
vx[x, y] += (p[x - 1, y] - p[x, y]) / wavespeed;
vy[x, y] += (p[x, y - 1] - p[x, y]) / wavespeed;

Så tilføjer jeg scenedata ved at påvirke hastighederne. Dæmpning sker ved at
dividere hastighederne med en konstant de steder jeg har tegnet dæmpning.
vx[x, y] /= 1.07F;

Og så beregner jeg positionerne ved at kigge på de 4 omkringliggende
hastigheder.
p[x, y] += (vx[x, y] - vx[x + 1, y] + vy[x, y] - vy[x, y + 1]) / wavespeed;

"wavespeed" er er 2 som udgangspunkt. Hvis jeg gør den større, udbreder
bølgerne sig langsommere. Og hvis den bliver ret meget under 2 går det hele
i sving.

Jeg har lavet absorberende kanter på fladen ved at sætte hastigheden lig med
den ene nabo-position. Her for højre side.
vx[xsize, y] = p[xsize - 1, y];

Men det virker ikke helt perfekt, der er noget reflektion tilbage, især ved
høje frekvenser. Det problem ved jeg ikke om kan løses, men det kunne være
fedt at slippe for den 'falske efterklang' fra kanterne.

Jeg har tidligere lavet det samme i VB6, skulle lige hilse og sige at det
kører en DEL bedre i C#!

--
Ulrik Smed
Aarhus, Denmark



Nick Janson (30-03-2008)
Kommentar
Fra : Nick Janson


Dato : 30-03-08 11:34

> Men jeg undrede mig over den perfekte lavpasvirkning der er i 'luften'
> alene, i 1D versionen. Man kan også se 'ripples' i 2D, når man slår an med
> en skarp puls.

Et fjedersystem som beskrevet vil altid være low pass. Spørgsmålet er så
hvor skarp grænsen er og hvor den ligger.
Du har et vidst tidsskridt h og sampler dit system deri. Hvis du har
frekvenser med periodelængder under h/2 så forsvinder de.

> Jeg ved ikke hvad min integrator hedder, men mit program (C#) har 3
> arrays, et positions (eller pressure) array, og et array for x-hastighed
> og et for y-hastighed.

X og y hastighed? Jeg ville have troet du lod partiklerne svinge op og ned i
plnet og ikke i planet. Da jeg legede dem denne slags, da svingede
partiklerne kun langs planets normalvektor. Det er også sådan det gøres på
det link jeg sendte sidst. Hvad er x-y-hastigheden?

Hvis man ønsker at se på trykket over tid, så vil man med helt simpel
integration (explicit euler) sige tryk(t+h)=tryk(t)+h*tryk'(t), hvor
tryk'(t) er hastigheden af trykændringen til tiden t. du vil også sige
tryk'(t+h) = tryk'(t) + h*tryk''(t). Trykændringen ændrer sig altså også.
Med fjedersystemer, som jeg gjorde det, så er trykket partiklens højde over
eller under planet. Ændringen i tryk er så dens hastighed langs planets
normalvektor og ændringen i ændring af tryk er partiklens acceleration som
er forskel i højde mellem partikel og nabopartikel ganget en konstant
(fjederkonstanten) divideret med partiklens masse.

Hvis man forestiller sig et gitter, i et plan, af partikler som er forbundet
til deres naboer med fjedre og du trækker en partikel op fra planet, så vil
den hive sine naboer med og de hiver deres etc. Bemærk at partikler kun
bevæger sig normalt til planet og ikke i planet.

Din metode virker lidt adskildt den matematiske beskrivelse, hvor man jo
opskriver en tilstandsvektor og så stepper i tiden. Jeg tror du vil få
vanskeligt ved at udvide den med bedre integrationsmetoden. Overvej evt at
beskrive partiklerne som en vektor med position og hastighed (begge dele i
1D) samt masse. Beskriv få fjedrene som skalarer der er fjederkonstanterne.
Du kan så implementere forskellige integrationsmetoder helt efter bogen.

det simple eksempel er stadig explicit euler (findes der en mere simpel?)
hvor du siger
system(t+h) = system(t)+h*system'(t)
Hvis dit system er vektoren
system = [pos0 vel0 pos1 vel1 pos2 vel2 ... posn veln]^T
så er den afledte
system' =[vel0 f0/m0 vel1 f1/m1 vel2 f2/m2 ... veln fn/mn]^T
f0 er alle krafterne på partikel 0 som afhænger af dens og dens naboers
positioner samt fjederkonstanterne for nabofjedrene.
Kraften på p0 fra p1 (som er en nabo) er så
f01 = (p0.pos-p1.pos-hvilelængde)*fjederkonstant
Du summerer alle krafter for alle partikler og har derefter den afledte
systemvektor og kan steppe.

Det er, så vidt jeg kan se, en korrekt beskrivelse af systemet.

> "wavespeed" er er 2 som udgangspunkt. Hvis jeg gør den større, udbreder
> bølgerne sig langsommere. Og hvis den bliver ret meget under 2 går det
> hele i sving.

Svingningerne er netop et udtryk for numerisk ustabilitet som simple
integrationsmetoder generelt vil opleve.

> Jeg har lavet absorberende kanter på fladen ved at sætte hastigheden lig
> med den ene nabo-position. Her for højre side.
> vx[xsize, y] = p[xsize - 1, y];

At sætte kantens hastighed til nabopartikelns hastighed svarer til at sætte
systemets (trykkets) første afledte lig værdien inde i domænet. Det er
Neumann grænsebetingelser
http://en.wikipedia.org/wiki/Neumann_boundary_condition og betyder i
princippet at du simulerer på et system hvor partiklerne udenfor dit område
alle har samme hastighed som din sidste partikel. For mig at se burde det
virke, mne måske jeg bare ikke tænker nok over det, for du siger jeg det
ikke gør,

> Men det virker ikke helt perfekt, der er noget reflektion tilbage, især
> ved høje frekvenser. Det problem ved jeg ikke om kan løses, men det kunne
> være fedt at slippe for den 'falske efterklang' fra kanterne.

Hvis jeg skal slynge et gæt ud, relateret til mine egne partikler, så er
problemet måske at nok sætter du hastigheden til det samme som sidste
simulerede partikel, men det tager jo tid for den derefter at komme også
flytte sig. Du har en partikel med positionen 10 og hastigheden 5. Naboen
(kanten) får så også hastigheden 5 men har stadig en eller anden position
(hvad er den egentlig) som eksempelvis er 0. I næste skridt vil kanten nok
have samme fart som naboen, men ikke samme position og dermed påvirker den
stadig...?

> Jeg har tidligere lavet det samme i VB6, skulle lige hilse og sige at det
> kører en DEL bedre i C#!

Kender det. Jeg kører også en del bedre i C#


Ulrik Smed (30-03-2008)
Kommentar
Fra : Ulrik Smed


Dato : 30-03-08 20:34

Nick Janson wrote:
>> Men jeg undrede mig over den perfekte lavpasvirkning der er i
>> 'luften' alene, i 1D versionen. Man kan også se 'ripples' i 2D, når
>> man slår an med en skarp puls.
>
> Et fjedersystem som beskrevet vil altid være low pass. Spørgsmålet er
> så hvor skarp grænsen er og hvor den ligger.
> Du har et vidst tidsskridt h og sampler dit system deri. Hvis du har
> frekvenser med periodelængder under h/2 så forsvinder de.

Det virker som om grænsen er ideelt skarp, og at grænsefrekvensen flytter
sig nedad jo længere man kommer fra kilden. Hvordan forholder det sig med
rigtig lyd i luft? Er der en tilsvarende grænsefrekvens der er afhængig af
afstanden?

> X og y hastighed? Jeg ville have troet du lod partiklerne svinge op
> og ned i plnet og ikke i planet. Da jeg legede dem denne slags, da
> svingede partiklerne kun langs planets normalvektor. Det er også
> sådan det gøres på det link jeg sendte sidst. Hvad er x-y-hastigheden?

Hm, ja det er egentlig noget vrøvl at kalde dem x og y hastighed, kan jeg
godt se. Jeg har jo kun én position til hver partikel, så de svinger kun i
én retning. Jeg lægger jo begge hastighederne til samme position. Nu har jeg
lavet en funktion der kun bruger 2 arrays. Det virker også, men jeg har lidt
problemer med at mine vægge 'opsuger' noget energi fra bølgen, så man ikke
kan tegne et smalt rør og sende en bølge igennem. Det var vist egentlig
derfor jeg oprindelig lavede x og y 'hastigheder', for så kunne jeg lave
lodrette og vandrette vægge der virkede rigtigt. Men det burde nu også være
muligt med den anden model, det må der lige leges lidt mere med.

> Hvis man ønsker at se på trykket over tid, så vil man med helt simpel
> integration (explicit euler) sige tryk(t+h)=tryk(t)+h*tryk'(t), hvor
> tryk'(t) er hastigheden af trykændringen til tiden t. du vil også sige
> tryk'(t+h) = tryk'(t) + h*tryk''(t). Trykændringen ændrer sig altså
> også. Med fjedersystemer, som jeg gjorde det, så er trykket
> partiklens højde over eller under planet. Ændringen i tryk er så dens
> hastighed langs planets normalvektor og ændringen i ændring af tryk
> er partiklens acceleration som er forskel i højde mellem partikel og
> nabopartikel ganget en konstant (fjederkonstanten) divideret med
> partiklens masse.
> Hvis man forestiller sig et gitter, i et plan, af partikler som er
> forbundet til deres naboer med fjedre og du trækker en partikel op
> fra planet, så vil den hive sine naboer med og de hiver deres etc.
> Bemærk at partikler kun bevæger sig normalt til planet og ikke i
> planet.

Ja, det må være det jeg har lavet nu. Hastigheden laver jeg ved at lægge de
4 nabopositioner sammen og dele med 4, og så trække positionen fra, for at
få en samlet fjederkraft fra de 4 nabopunkter.
v[x, y] += ((p[x - 1, y] + p[x + 1, y] + p[x, y - 1] + p[x, y + 1])/4 - p[x,
y]);

Og positionen lægger jeg blot hastigheden til.
p[x, y] += v[x, y];

> Din metode virker lidt adskildt den matematiske beskrivelse, hvor man
> jo opskriver en tilstandsvektor og så stepper i tiden. Jeg tror du
> vil få vanskeligt ved at udvide den med bedre integrationsmetoden.
> Overvej evt at beskrive partiklerne som en vektor med position og
> hastighed (begge dele i 1D) samt masse. Beskriv få fjedrene som
> skalarer der er fjederkonstanterne. Du kan så implementere
> forskellige integrationsmetoder helt efter bogen.
> det simple eksempel er stadig explicit euler (findes der en mere
> simpel?) hvor du siger
> system(t+h) = system(t)+h*system'(t)
> Hvis dit system er vektoren
> system = [pos0 vel0 pos1 vel1 pos2 vel2 ... posn veln]^T
> så er den afledte
> system' =[vel0 f0/m0 vel1 f1/m1 vel2 f2/m2 ... veln fn/mn]^T
> f0 er alle krafterne på partikel 0 som afhænger af dens og dens
> naboers positioner samt fjederkonstanterne for nabofjedrene.
> Kraften på p0 fra p1 (som er en nabo) er så
> f01 = (p0.pos-p1.pos-hvilelængde)*fjederkonstant
> Du summerer alle krafter for alle partikler og har derefter den
> afledte systemvektor og kan steppe.
>
> Det er, så vidt jeg kan se, en korrekt beskrivelse af systemet.

Lyder meget som det jeg har nu, bortset fra at jeg har sparet
fjederkonstanterne og partikelmasserne. Men de svarer vel til at jeg delte
med min 'wavespeed' før.

>> "wavespeed" er er 2 som udgangspunkt. Hvis jeg gør den større,
>> udbreder bølgerne sig langsommere. Og hvis den bliver ret meget
>> under 2 går det hele i sving.
>
> Svingningerne er netop et udtryk for numerisk ustabilitet som simple
> integrationsmetoder generelt vil opleve.

Hm, denne integrationsmetode, er det blot det at man ganger med en
fjederkonstant eller en partikelmasse, så man kun får overført en del af
fjederkraften eller hastigheden for hvert tidsskridt?

> Hvis jeg skal slynge et gæt ud, relateret til mine egne partikler, så
> er problemet måske at nok sætter du hastigheden til det samme som
> sidste simulerede partikel, men det tager jo tid for den derefter at
> komme også flytte sig. Du har en partikel med positionen 10 og
> hastigheden 5. Naboen (kanten) får så også hastigheden 5 men har
> stadig en eller anden position (hvad er den egentlig) som eksempelvis
> er 0. I næste skridt vil kanten nok have samme fart som naboen, men
> ikke samme position og dermed påvirker den stadig...?

Det lyder sandsynligt. Det virker som om det netop er forskellen mellem to
tidsskridt der reflekteres. Den forskel er mindre jo lavere frekvensen er,
så derfor reflekteres høje frekvenser mest. En stor, blød bølge absorberes
næsten perfekt.

Hm, jeg føler mig lidt handicappet omkring alle de matematiske begreber og
måder at skrive det på. Har ikke rigtig 'forstand' på matematik på det plan.
Kan bedre fatte nogle linier kode, haha.

>> Jeg har tidligere lavet det samme i VB6, skulle lige hilse og sige
>> at det kører en DEL bedre i C#!
>
> Kender det. Jeg kører også en del bedre i C#

Fedt, så kan vi jo udveksle kode. Jeg kører med den gratis C# 2005
Express.

--
Ulrik Smed
Aarhus, Denmark



Nick Janson (30-03-2008)
Kommentar
Fra : Nick Janson


Dato : 30-03-08 21:16

> Det virker også, men jeg har lidt problemer med at mine vægge 'opsuger'
> noget energi fra bølgen, så man ikke kan tegne et smalt rør og sende en
> bølge igennem. Det var vist egentlig derfor jeg oprindelig lavede x og y
> 'hastigheder', for så kunne jeg lave lodrette og vandrette vægge der
> virkede rigtigt. Men det burde nu også være muligt med den anden model,
> det må der lige leges lidt mere med.

Jeg lavede mine reflekterende vægge ved at fastholde kanternes position. Det
kan let gøres ved at sætte massen til uendelig. Hvis du vil lave absorbtion
kan du måske bare slette yderste fjeder, så sidste reele partikel ikke er
forbundet til andet end naboerne inde i domænet. Det burde svare til en
direlict grænsebetingelse hvor alle partikler til højre for domænet har
samme position som sidste partikel indenfor.

> Ja, det må være det jeg har lavet nu. Hastigheden laver jeg ved at lægge
> de 4 nabopositioner sammen og dele med 4, og så trække positionen fra, for
> at få en samlet fjederkraft fra de 4 nabopunkter.
> v[x, y] += ((p[x - 1, y] + p[x + 1, y] + p[x, y - 1] + p[x, y + 1])/4 -
> p[x, y]);
>
> Og positionen lægger jeg blot hastigheden til.
> p[x, y] += v[x, y];

Bare for en ordens skyld... du løser dette som en helhed, ikke? Ikke noget
med at løse for en partikel og så for næste, som nu har en nabo der er
ændret etc? Nævner det bare fordi jeg en gang (for længe længe siden)
forglemte mig selv og gjorde det sådan.


> Lyder meget som det jeg har nu, bortset fra at jeg har sparet
> fjederkonstanterne og partikelmasserne. Men de svarer vel til at jeg delte
> med min 'wavespeed' før.

Pointen med at beskrive det i mere gængse termer er egentlig at du så uden
det store besvær kan prøve med andre integrationsmetoder.

> Hm, denne integrationsmetode, er det blot det at man ganger med en
> fjederkonstant eller en partikelmasse, så man kun får overført en del af
> fjederkraften eller hastigheden for hvert tidsskridt?

Jo mindre skridt, jo længere tid tager det, og jo flere samples er der
undervejs, før en fjeder trækker sig sammen.
Hvis du har en partikel på positionen x=10m og massen m=1kg samt farten
v=0m/s forbundet til et punkt i 0 med en fjeder med fjederkonstanten
k=1000N/m og hvilelængden L=0, så vil du have en fjederkraft på -10m*1000N/m
= -10000N
Hvis du har tidsskridt h=1s så siger du
x(t+h) = x(t) + h* v(t)
v(t+h) = v(t) + h* a(t)

hvor

a(t) = -10000N/1kg = 10000m/s^2

det giver
x(t+h) = 10 + 1s*0m/s = 10m
v(t+h) = 0 + 1s*-10000m/s/s = -10000m/s

Positionen er stadig uændret så kraften er den samme og accelerationen er
uændret.

Næste skridt er
x = 10m + 1s*-10000m/s = -9990m
v = -10000m/s + 1s*-10000m/s/s = -20000m/s

Nu er positionen -9990m og kraften er 9990m*1000N/m = 9990000N

Accelerationen er nu ekstram til højre hvor den før var til venstre. Den er
også markant større. Når den så svinger tilbage næste gang så bliver den
endnu værre.
pointen er at man har antaget at accelerationen er konstant over hele
tidsskridtet hvilket er forkert. Den falder hurtigt og ændrer fortegn så
snart x=0 passeres. Dette er problemet med ustabilitet ved store tidsskridt.
Du må undskylde hvis jeg tærsker langhalm, på sjusket vis endda, på noget du
godt ved.


> Hm, jeg føler mig lidt handicappet omkring alle de matematiske begreber og
> måder at skrive det på. Har ikke rigtig 'forstand' på matematik på det
> plan. Kan bedre fatte nogle linier kode, haha.

ingen grund til at føle dig handikappet. Du har ikke mistet nogle ben som
ikke kan gro ud igen. Mit niveau er bestemt ikke imponerende og det er kun
indenfor de senere år jeg har kunet forstå emnet godt nok til at ikke hive
hår ud i tide og utide. Nu har jeg så andre ting jeg kan blive frustreret
over.
Det eneste væsentlige er lige nu bare at forstå princippet bag at beskrive
et system som en vektor og at forstå hvad de "afledte" af en funktion er.

så giver notation som f(t+h) = f(t)+h*f'(t) mening og hvis du så kan løse
ligningssystemer så vil selv umiddelbart underlig formulering som
f(t+h) = f(t)+h*f'(t+h) hvor næste positions værdi beregnes ud fra den
afledte i næste punkt, som du jo stadig ikke kender, give fin mening.

Du skal da også være velkommen til at stille sprørgsmål hvis du går i gang
med "oversættelsen" af din kode. Jeg har praktisk erfaring med den model og
der er nok andre her som bedre kan besvare hvad der er realistisk og hvad
der er artifacts.


Ulrik Smed (30-03-2008)
Kommentar
Fra : Ulrik Smed


Dato : 30-03-08 23:33

Nick Janson wrote:

> Jeg lavede mine reflekterende vægge ved at fastholde kanternes
> position. Det kan let gøres ved at sætte massen til uendelig.

Ja, det giver en reflektion, men den bliver negativ. Hvis jeg bruger det til
at lave et rør med, så svarer det til en membran der sidder fast i rørets
vægge.

> Hvis du
> vil lave absorbtion kan du måske bare slette yderste fjeder, så
> sidste reele partikel ikke er forbundet til andet end naboerne inde i
> domænet. Det burde svare til en direlict grænsebetingelse hvor alle
> partikler til højre for domænet har samme position som sidste
> partikel indenfor.

Det giver en positiv reflektion. Faktisk netop det jeg skal bruge til vægge.
Men når jeg kun har én array med hastigheder, i stedet for de x og y 'nogle'
jeg brugte før (er vel egentlig bare lagrede fjeder-spændinger), så kan jeg
ikke komme til at 'klippe en fjeder over', ved at nulstille
fjederspændingen. Det var det der var årsagen til at jeg havde x og y arrays
før. Problemet bunder nok i at jeg - for perfomancens skyld - først beregner
fjederspændinger uden at kigge på scene-data, derefter hiver i de spændinger
scene-dataene peger på (nulstiller dem for at få en væg), og så beregner
positioner, igen uden at kigge på scene-data.

> Bare for en ordens skyld... du løser dette som en helhed, ikke? Ikke
> noget med at løse for en partikel og så for næste, som nu har en nabo
> der er ændret etc? Nævner det bare fordi jeg en gang (for længe længe
> siden) forglemte mig selv og gjorde det sådan.

Hehe, jeps, har været der selv. Jojo, først lave alle hastighederne,
og derefter alle positionerne.

> Jo mindre skridt, jo længere tid tager det, og jo flere samples er der
> undervejs, før en fjeder trækker sig sammen.
> Hvis du har en partikel på positionen x=10m og massen m=1kg samt
> farten v=0m/s forbundet til et punkt i 0 med en fjeder med
> fjederkonstanten k=1000N/m og hvilelængden L=0, så vil du have en
> fjederkraft på -10m*1000N/m = -10000N
> Hvis du har tidsskridt h=1s så siger du
> x(t+h) = x(t) + h* v(t)
> v(t+h) = v(t) + h* a(t)
>
> hvor
>
> a(t) = -10000N/1kg = 10000m/s^2
>
> det giver
> x(t+h) = 10 + 1s*0m/s = 10m
> v(t+h) = 0 + 1s*-10000m/s/s = -10000m/s
>
> Positionen er stadig uændret så kraften er den samme og
> accelerationen er uændret.
>
> Næste skridt er
> x = 10m + 1s*-10000m/s = -9990m
> v = -10000m/s + 1s*-10000m/s/s = -20000m/s
>
> Nu er positionen -9990m og kraften er 9990m*1000N/m = 9990000N
>
> Accelerationen er nu ekstram til højre hvor den før var til venstre.
> Den er også markant større. Når den så svinger tilbage næste gang så
> bliver den endnu værre.
> pointen er at man har antaget at accelerationen er konstant over hele
> tidsskridtet hvilket er forkert. Den falder hurtigt og ændrer fortegn
> så snart x=0 passeres. Dette er problemet med ustabilitet ved store
> tidsskridt. Du må undskylde hvis jeg tærsker langhalm, på sjusket vis
> endda, på noget du godt ved.

Ja, det kan jeg godt forstå, det var det der skete når min 'wavespeed' blev
for lav. Det lidt sjove er at hvis man laver en 1D række, så bliver den
perfekt stabil og uden lavpasvirkning hvis man netop overfører hele energien
på et step, uden at integrere. Ligesom den her Flash-ting, når man sætter
"damping" til 0.
http://phet.colorado.edu/simulations/stringwave/stringWave.swf

> Du skal da også være velkommen til at stille sprørgsmål hvis du går i
> gang med "oversættelsen" af din kode. Jeg har praktisk erfaring med
> den model og der er nok andre her som bedre kan besvare hvad der er
> realistisk og hvad der er artifacts.

Takker for det, jeg leger lidt videre og spø'r hvis det er.

--
Ulrik Smed
Aarhus, Denmark



Carsten Svaneborg (31-03-2008)
Kommentar
Fra : Carsten Svaneborg


Dato : 31-03-08 00:29

Ulrik Smed wrote:
> Jeg sidder netop selv og pjatter lidt med noget bølgesimulering efter
> samme model med lodder og fjedre.

Det burde være muligt at løse sådan en model analysisk. Fjeder modeller
kan generelt opskrives som en matrix så E=0.5 x^T C x + 0.5 v^T M v

Hvor xi er den i'te komponent af vektoren af udsving af punkt i, v er
vektoren af hastigheder, og C er en symmetrisk matrix så Cij er halvdelen
af fjeder konstanten mellem punkt i og j, mens M er en diagonal matrice
med masserne af de enkelte punkter.

Det er sent, men diagonaliseres C matricen, så afkobles de enkelte
svingning modes som findes i gitteret, den potentielle energi bliver
så forsimplet til Epot = 0.5 sum_p k_p a_p^2. Hvor k_p er fjeder
konstanten for den p'te mode (egenværdien af C matricen), mens
a_p er amplituden af den mode.

Du bør så få en bevægelsesligning hvor hver mode, men tricket er at disse
nu er uafhængige. Så fra at have N koblede differential ligninger, har du
N uafhængige differential ligninger, som kan løses analytisk.

Alt arbejdet ligger så i at diagonaliserer matricen for hvem der er
forbundet med hvem...

--
Mvh. Carsten Svaneborg
http://gauss.ffii.org softwarepatent database

Nick Janson (30-03-2008)
Kommentar
Fra : Nick Janson


Dato : 30-03-08 11:53

Jeg tror jeg misforstod noget.
Jeg læste om FEM i forbindelse med simulation af "soft bodies" hvor man
bruger tetraeder som elementer og ser på elasticitet og regner løs.
Det er jo ikke FEM, det er blot een anvendelse af det. Det havde været som
hvis jeg troede finite difference havde med en speciel simuleringsopgave at
gøre og ikke bare var en måde at estimere det afledte. Hm. Ja, så forstår
jeg nok bedre at ingen har svaret.



Søg
Reklame
Statistik
Spørgsmål : 177459
Tips : 31964
Nyheder : 719565
Indlæg : 6408186
Brugere : 218881

Månedens bedste
Årets bedste
Sidste års bedste