WICHTIG: Der Betrieb von goMatlab.de wird privat finanziert fortgesetzt. - Mehr Infos...

Mein MATLAB Forum - goMatlab.de

Mein MATLAB Forum

 
Gast > Registrieren       Autologin?   

Partner:




Forum
      Option
[Erweitert]
  • Diese Seite per Mail weiterempfehlen
     


Gehe zu:  
Neues Thema eröffnen Neue Antwort erstellen

Nicht-linearer 1-Massen-Schwinger

 

Harald
Forum-Meister

Forum-Meister


Beiträge: 24.448
Anmeldedatum: 26.03.09
Wohnort: Nähe München
Version: ab 2017b
     Beitrag Verfasst am: 24.04.2017, 19:54     Titel:
  Antworten mit Zitat      
Hallo,

dazu wird es am besten sein, wenn du die Lösung der x-Werte wieder in die Funktion einsetzt. Etwa so:

Code:
allK = zeros(size(y1, 1), 1);
for I = 1:size(y1, 1)
allK = Kfun(y1(I,:), k);
end


Grüße,
Harald
Private Nachricht senden Benutzer-Profile anzeigen


Codename
Themenstarter

Forum-Anfänger

Forum-Anfänger


Beiträge: 27
Anmeldedatum: 16.07.15
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 27.04.2017, 15:50     Titel:
  Antworten mit Zitat      
Hallo Harald,

zunächst vielen Dank für deine wertvolle Hilfe.

Ich möchte nun mein System erweitern und die Masse zwischen 2 Stabfedern einspannen, die jeweils nur Druckkräfte auf die Masse ausüben. Die zusätzlichen Kräfte durch die Stabfedern sind Fl(x) und Fr(x)
Code:

[t1,y1] = ode45(@(t,x) [x(2); (-x(2)*c-k*x(1)-Fr(x)+Fl(x)+F(t))/m] ,t ,[0 0]);


Zur Berechnung dieser Kräfte habe ich ein Werkstoffmodell zur Verfügung, welches mir in Abhängigkeit von der Dehnung eps der Stabfedern eine Spannung sigma ausgibt.
Der Spannung-Dehnungs-Verlauf dieser Stabfedern beschreibt dabei eine Hysterese.
Über den Zusammenhang Kraft = Spannung*Fläche möchte ich dann Fr und Fl berechnen.

Die Spannung der jeweiligen Stabfeder wird durch ein zeitinkrementelles Prädiktor-Korrektor-Verfahren berechnet. Bei diesem Verfahren wird die Spannung zum Zeitpunkt tn aus einem Prädiktor und einem Korrektor-Schritt berechnet.

Die Herausforderung liegt nun in der Umsetzung dieses Vorgehens. Wie verknüpfe ich das Lösungsverfahren der DGL (ode45) im Programm "solver" mit dem Prädiktor-Korrektor-Verfahren im Programm "Werkstoffmodell"?

Werkstoffmodell.m
 Beschreibung:

Download
 Dateiname:  Werkstoffmodell.m
 Dateigröße:  20.69 KB
 Heruntergeladen:  323 mal
Solver.m
 Beschreibung:

Download
 Dateiname:  Solver.m
 Dateigröße:  3.05 KB
 Heruntergeladen:  301 mal
Main_test.m
 Beschreibung:

Download
 Dateiname:  Main_test.m
 Dateigröße:  2.78 KB
 Heruntergeladen:  305 mal
Private Nachricht senden Benutzer-Profile anzeigen
 
Harald
Forum-Meister

Forum-Meister


Beiträge: 24.448
Anmeldedatum: 26.03.09
Wohnort: Nähe München
Version: ab 2017b
     Beitrag Verfasst am: 27.04.2017, 16:44     Titel:
  Antworten mit Zitat      
Hallo,

die Klasse Werkstoffmodell ist sehr umfangreich. Ich habe keine Ahnung, wie sie genutzt werden soll. Daher kann ich erst recht nicht sagen, wie sie mit ode45 verknüpft werden kann.
Solltest du den Autor der Klasse kennen, spreche doch mal ihn darauf an.

Grüße,
Harald
Private Nachricht senden Benutzer-Profile anzeigen
 
Codename
Themenstarter

Forum-Anfänger

Forum-Anfänger


Beiträge: 27
Anmeldedatum: 16.07.15
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 28.04.2017, 11:25     Titel:
  Antworten mit Zitat      
Hallo Harald,

die Vorgehensweise des Programms soll wie folgt aussehen:

Zu Beginn wird die Funktion "Werkstoffmodell.fnSetPrestress" ausgeführt. Diese berechnet die Vorspannung und -dehnung in den Stabfedern zum Zeitpunkt t=0. Im Anschluss daran wird die Verlagerung x der Masse zum Zeitpunkt t=0 berechnet. In meinem Fall beträgt diese 0, da noch keine äußere Kräften wirken und die Startbedingung x(0)=0 ist.
Im nächsten Zeitschritt (t1) wirkt dann die Anregungskraft F(t1) auf das System, wodurch es ausgelenkt wird. Die Auslenkung der Masse führt dazu, dass die rechte Stabfeder weiter gestaucht und die linke entlastet wird. Dies führt dazu, dass die Kraft Fr(t1)>Fl(t1) ist.
Zur Berechnung von Fr(t1) und Fl(t1) wird zunächst die Dehnung eps_r(t1) und eps_l(t1) berechnet.

Die Funktion "Werkstoffmodell.fnCurrentState(eps_r_n,eps_l_n,t_n)" berechnet dann ein delta_eps und delta_t aus der Differenz der jeweiligen Werte zum aktuellen Zeitpunkt t1 und dem vorherigen Zeitpunkt t0. Diese Werte werden dann in die Funktionen "Werkstoffmodell.fnRunPredictorStep; Werkstoffmodell.fnEvaluateFlowCondition; Werkstoffmodell.fnUpdateTemperature;" eingesetzt, um die Spannung der Stabfedern sig_r und sig_l zum Zeitpunkt t1 zu berechnen.
Die Stabfederkräfte Fr und Fl können nun über F=sig*A berechnet werden und das Runge-Kutta Verfahren kann bis zum Zeitpunkt t1 die DGL lösen und die Verlagerung der Masse x(t1) ausgeben.
Im letzten Schritt wird die Funktion "Werkstoffmodell.UpdateStateVariables" aufgerufen, um die Zustandsgrößen abzuspeichern.

Nachdem die zuvor beschriebenen Funktionen ausgeführt wurden, wird der ganze Ablauf für den nächsten Zeitschritt (t2) durchgeführt, solange bis t_end erreicht wird.

Ich hoffe das Vorgehen ist so wie beschrieben nachvollziehbar. Über eine Hilfestellung zur Implementierung dieses Vorgehens würde ich mich sehr freuen.

Viele Grüße
Private Nachricht senden Benutzer-Profile anzeigen
 
Codename
Themenstarter

Forum-Anfänger

Forum-Anfänger


Beiträge: 27
Anmeldedatum: 16.07.15
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 01.05.2017, 12:45     Titel:
  Antworten mit Zitat      
Hallo,
Ich möchte das Thema noch einmal pushen, da ich bisher keine Lösung zu meinem Problem gefunden habe.
Ist das von mir beschriebene Vorgehen so überhaupt mit dem ode-solver (Runge-Kutta) möglich?

An der Umsetzung scheitert es aktuell daran, dass die Klasse "Werkstoffmodell" zu jedem Zeitpunkt aufgerufen werden muss, wohingegen die Funktion Ode45 einmalig für den gesamten Zeitbereich aufgerufen wird. Liege ich damit richtig?

Um dieses Problem zu entgehen, könnte ich alternativ zum ode-solver die Newmark-beta-Methode implementieren. Diese Methode wird wie die Klasse "Werkstoffmodell" zu jedem Zeitpunkt aufgerufen. Macht diese Methode in meinem Fall mehr Sinn?

Über eine Rückmeldung zu diesem Thema würde ich mich sehr freuen.
Private Nachricht senden Benutzer-Profile anzeigen
 
Harald
Forum-Meister

Forum-Meister


Beiträge: 24.448
Anmeldedatum: 26.03.09
Wohnort: Nähe München
Version: ab 2017b
     Beitrag Verfasst am: 01.05.2017, 14:25     Titel:
  Antworten mit Zitat      
Hallo,

Zitat:
An der Umsetzung scheitert es aktuell daran, dass die Klasse "Werkstoffmodell" zu jedem Zeitpunkt aufgerufen werden muss, wohingegen die Funktion Ode45 einmalig für den gesamten Zeitbereich aufgerufen wird. Liege ich damit richtig?

Du kannst die Klasse ja in der "rechten Seite" der DGL aufrufen, und dann wird sie auch in jedem Zeitschritt aufgerufen. Problematisch kann sein, dass bei ode45 auch mal ein Schritt abgelehnt wird, wenn Fehlertoleranzen nicht eingehalten werden.

Erster Schritt dürfte sein, die "rechte Seite" mal komplett in eine Funktion auszulagern statt ein Anonymous Function Handle zu verwenden.

Grüße,
Harald
Private Nachricht senden Benutzer-Profile anzeigen
 
Codename
Themenstarter

Forum-Anfänger

Forum-Anfänger


Beiträge: 27
Anmeldedatum: 16.07.15
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 08.05.2017, 09:15     Titel:
  Antworten mit Zitat      
Hallo Harald,

ich habe die Funktion nun so geschrieben:

Code:
function [Fres,Fr,Fl,time,l_r,l_l,eps_n_r,eps_n_l,sig_n_r,sig_n_l] = Fresfun(x,t,Werkstoffmodell,Solver)

if x(1) >= 0
eps_n_r = (Solver.l_0 - Solver.l_v + x(1)) / Solver.l_0;
eps_n_l = (Solver.l_0 - Solver.l_v - x(1)) / Solver.l_0;
l_r = Solver.l_v - x(1);
l_l = Solver.l_v + x(1);
display ('x(1)>=0');

    if l_r >= Solver.l_0
    warning('obj.l_r >= obj.l_0');
    warning('Berechnung nicht korrekt!');
    elseif l_l >= Solver.l_0
    warning('l_l >= Solver.l_0');
    warning('Berechnung nicht korrekt!');    
    end
       
elseif x(1) < 0
eps_n_r = (Solver.l_0 - Solver.l_v + x(1)) / Solver.l_0;
eps_n_l = (Solver.l_0 - Solver.l_v - x(1)) / Solver.l_0;
l_r = Solver.l_v - x(1);
l_l = Solver.l_v + x(1);
display ('x(1)<0');

    if l_r >= Solver.l_0
    warning('obj.l_r >= obj.l_0');
    warning('Berechnung nicht korrekt!');
    elseif l_l >= Solver.l_0
    warning('l_l >= Solver.l_0');
    warning('Berechnung nicht korrekt!');    
    end
   
end
 
display('x'); display(x(1));
display('t'); display(t);
display('dt'); display(Solver.dt);

time = t;

%% Methoden des Werkstoffmodells
% tic
Werkstoffmodell.fnComputeCurrentState(eps_n_r,eps_n_l,t);
Werkstoffmodell.fnRunPredictorStep;
Werkstoffmodell.fnEvaluateFlowCondition;
Werkstoffmodell.fnUpdateTemperature;
% toc

sig_n_r = Werkstoffmodell.sig_n_r;
sig_n_l = Werkstoffmodell.sig_n_l;


Fr = Werkstoffmodell.sig_n_r*Solver.Area*1e6; %Probenkraft rechts [N]
Fl = Werkstoffmodell.sig_n_l*Solver.Area*1e6; %Probenkraft links [N]
Fres = Fl-Fr; %Probenkraft links - rechts [N]

display ('Fr');
display (Fr);
display ('Fl');
display (Fl);
display ('Fres');
display (Fres);

% Werkstoffmodell.fnUpdateStateVariables;
display ('--------------------------------');
 


Meine Fragen dazu:

1.) Wie du bereits erwähnt hast rechnet ode45 je nach Fehlertoleranz mit unterschiedlichen Zeitschritten dt. Kann ich diese Zeitschritte auch als konstant vorgeben?
Dieser Ansatz funktioniert nicht so ganz, da dt von ode45 verändert wird:
Code:
t=0:obj.dt:obj.t_end;   % time scale
[obj.t1,obj.y1] = ode45(@(t,x) [x(2); (-x(2)*obj.c-obj.k*x(1)+Fresfun(x,t,obj.Werkstoffmodell,obj)+obj.F(t))/obj.m] ,t [0 0]); % Neue DGL


2.) Wie kann ich die von ode45 verwendeten Zeitschritte dt abspeichern bzw. auslesen? Für meine Klasse "Werkstoffmodell" benötige ich die genaue Schrittweite.

Viele Grüße
Private Nachricht senden Benutzer-Profile anzeigen
 
Harald
Forum-Meister

Forum-Meister


Beiträge: 24.448
Anmeldedatum: 26.03.09
Wohnort: Nähe München
Version: ab 2017b
     Beitrag Verfasst am: 08.05.2017, 09:27     Titel:
  Antworten mit Zitat      
Hallo,

zu 1.: man könnte höchstens 'MinStep' und 'MaxStep' über odeset auf den gleichen Wert setzen.

zu 2.: du könntest dir die vorherigen Zeitschritte als persistent-Variable abspeichern.

Grüße,
Harald
Private Nachricht senden Benutzer-Profile anzeigen
 
Codename
Themenstarter

Forum-Anfänger

Forum-Anfänger


Beiträge: 27
Anmeldedatum: 16.07.15
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 08.05.2017, 14:22     Titel:
  Antworten mit Zitat      
Hallo Harald,

die Option "MinStep" kann ich nicht einstellen.
Wie hast du das gemeint?
Private Nachricht senden Benutzer-Profile anzeigen
 
Harald
Forum-Meister

Forum-Meister


Beiträge: 24.448
Anmeldedatum: 26.03.09
Wohnort: Nähe München
Version: ab 2017b
     Beitrag Verfasst am: 08.05.2017, 14:32     Titel:
  Antworten mit Zitat      
Hallo,

mein Fehler, es gibt nur ein MaxStep.
Das wäre aber so oder so nur eine Notlösung gewesen: es ist ja nicht sinnvoll, einen Löser mit variabler Schrittweite dazu zu zwingen, eine konstante Schrittweite beizubehalten.

Stattdessen sollte man dann ein Verfahren mit fester Schrittweite verwenden, z.B. Runge-Kutta 4. Ordnung, siehe z.B. hier
https://www.mathworks.com/matlabcen.....runge-kutta-4th-order-ode

Grüße,
Harald
Private Nachricht senden Benutzer-Profile anzeigen
 
Codename
Themenstarter

Forum-Anfänger

Forum-Anfänger


Beiträge: 27
Anmeldedatum: 16.07.15
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 17.05.2017, 16:57     Titel:
  Antworten mit Zitat      
Hallo Harald,

ich habe nun deinen Hinweis umgesetzt und verschiedene fixed step solver ausprobiert, die ich hier gefunden habe:
https://de.mathworks.com/matlabcent.....lver-in-matlab-8-0-r2012b


Während des Testens dieser Verfahren musste ich allerdings feststellen, dass ich mit diesen solvern nur sinnvolle Ergebnisse bekomme, wenn die Masse mind. 100 kg beträgt. Für meinen Anwendungsfall möchte ich allerdings eine Masse von 0.004 kg vorgeben. In diesem Fall divergiert allerdings die Simulation für alle fixed step solver und die Ergebnisse sind nicht mehr zu gebrauchen. Auch eine Reduzierung des Zeitschritts von 1/1000 s auf 1/100000 s konnte dieses Problem nicht beheben.

Daraufhin habe ich wieder mit variable-step solver herumprobiert und konnte dabei gute Ergebnisse bei der Verwendung der stiff-ode-solver (ode15s,ode23s,ode23tb) generieren. Die nonstiff-ode-solver (ode45,ode23,ode113) konnte ich ebenfalls nicht verwenden, da die Simulationsdauer sehr hoch war.

Da es sich in diesem Fall offensichtlich um ein numerisches Problem handelt, bin ich um jeder weitere Hilfe sehr dankbar.
Mir fehlt es momentan an der systematischen Herangehensweise an dieses Problem. Handelt es sich beispielsweise bei meinem System um eine steife DGL? Im Idealfall würde ich gerne wieder auf ein fixed-step-solver umschwenken. Welche Möglichkeiten (z.B. andere numerische Verfahren) bieten sich dafür an?

Grüße
Private Nachricht senden Benutzer-Profile anzeigen
 
Holger Gerach

Gast


Beiträge: ---
Anmeldedatum: ---
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 17.05.2017, 23:36     Titel:
  Antworten mit Zitat      
Hallo,

Ich hätte einen Extended Kalman Filter verwendet, den man leicht in einer m-function implementieren kann. Bei der gerigen Systemordnung dürfte das gut gehen. Alternative wäre ein Unscented Kalman Filter, falls die Schätzung nicht erwartungstreu genug ist.

Gruß Holger
 
Codename
Themenstarter

Forum-Anfänger

Forum-Anfänger


Beiträge: 27
Anmeldedatum: 16.07.15
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 18.05.2017, 10:13     Titel:
  Antworten mit Zitat      
Hallo Holger,

vielen Dank für den Hinweis.
Ich habe nun eine Funktion des Extended Kalman Filters unter http://de.mathworks.com/matlabcentr.....=5096219&tab=function gefunden.

Nun weiß ich allerdings nicht, wie ich diesen Filter in mein Programm implementieren kann. Habe ich es richtig verstanden, dass ich den Kalman-Filter verwende, nachdem ich meine DGL mit einem ode-solver gelöst habe?

Es wäre sehr hilfreich, wenn du mir eine Hilfestellung geben könntest.

Viele Grüße
Private Nachricht senden Benutzer-Profile anzeigen
 
Codename
Themenstarter

Forum-Anfänger

Forum-Anfänger


Beiträge: 27
Anmeldedatum: 16.07.15
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 20.05.2017, 15:40     Titel:
  Antworten mit Zitat      
Hallo zusammen,

hat noch jemand einen guten Hinweis zu meiner Problemstellung?

In der Literatur wird das von mir beschriebene System (nonlinear dynamics of a sdof oscillator with hysteresis) häufig durch einen "Newmark implicit algorithm" in Verbindung mit der "Newton-Raphson-Iteration" gelöst. Den Algorithmus habe ich mal in den Anhang gepackt.

Liefert dieses Verfahren überhaupt Vorteile gegenüber den in Matlab bereits implementiertetn ode-solvern (Runge-Kutta)?

Viele Grüße

Newmark_for_nonlinear_systems.pdf
 Beschreibung:

Download
 Dateiname:  Newmark_for_nonlinear_systems.pdf
 Dateigröße:  395.76 KB
 Heruntergeladen:  494 mal
Private Nachricht senden Benutzer-Profile anzeigen
 
Neues Thema eröffnen Neue Antwort erstellen

Gehe zu Seite Zurück  1, 2

Einstellungen und Berechtigungen
Beiträge der letzten Zeit anzeigen:

Du kannst Beiträge in dieses Forum schreiben.
Du kannst auf Beiträge in diesem Forum antworten.
Du kannst deine Beiträge in diesem Forum nicht bearbeiten.
Du kannst deine Beiträge in diesem Forum nicht löschen.
Du kannst an Umfragen in diesem Forum nicht mitmachen.
Du kannst Dateien in diesem Forum posten
Du kannst Dateien in diesem Forum herunterladen
.





 Impressum  | Nutzungsbedingungen  | Datenschutz | FAQ | goMatlab RSS Button RSS

Hosted by:


Copyright © 2007 - 2024 goMatlab.de | Dies ist keine offizielle Website der Firma The Mathworks

MATLAB, Simulink, Stateflow, Handle Graphics, Real-Time Workshop, SimBiology, SimHydraulics, SimEvents, and xPC TargetBox are registered trademarks and The MathWorks, the L-shaped membrane logo, and Embedded MATLAB are trademarks of The MathWorks, Inc.