Gelöste Aufgaben/Kig1: Unterschied zwischen den Versionen

Aus numpedia
Zur Navigation springen Zur Suche springen
Keine Bearbeitungszusammenfassung
 
(6 dazwischenliegende Versionen desselben Benutzers werden nicht angezeigt)
Zeile 13: Zeile 13:


<onlyinclude>
<onlyinclude>
[[Datei:Kig1-01.png|200x200px|left|mini|Lageplan]]
[[Datei:Kig1-01.png|250x250px|left|mini|Lageplan|alternativtext=]]
Damit der Schalter zuverlässig funktioniert soll die Kraft ''F'' an der Führung so gewählt werden, dass die vorgegebene Kontaktkraft ''K'' eingestellt wird.Gesucht ist die analytische Lösung für ein Euler-Bernoulli-Modell der elastischen Struktur.
Damit der Schalter zuverlässig funktioniert soll die Kraft ''F'' an der Führung so gewählt werden, dass die vorgegebene Kontaktkraft ''K'' eingestellt wird.
Gesucht ist die analytische Lösung für ein Euler-Bernoulli-Modell der elastischen Struktur.
</onlyinclude>
</onlyinclude>
Wir lösen dazu das Gleichungssystem zur Bestimmung der Integrationskonstanten der Differentialgleichung des Euler-Bernoulli-Balkens und weiterer Unbekannten und geben die grafischen Lösungen für ''w, w', M'' und ''Q'' an.
Wir lösen dazu das Gleichungssystem zur Bestimmung der Integrationskonstanten der Differentialgleichung des Euler-Bernoulli-Balkens und weiterer Unbekannten und geben die grafischen Lösungen für ''w, w', M'' und ''Q'' an.
Zeile 28: Zeile 29:
Aus dem Bild sehen wir
Aus dem Bild sehen wir


<math>\begin{array}{ccl}
::<math>\begin{array}{ccl}
w_1(a) = w_B\\
w_1(a) = w_B\\
w_2(\ell-a) = w_C = d\\
w_2(\ell-a) = w_C = d\\
Zeile 35: Zeile 36:
und wir arbeiten im folgenden mit
und wir arbeiten im folgenden mit


<math>\begin{array}{ccl}
::<math>\begin{array}{ccl}
\alpha&=&a/\ell\\
\alpha&=&a/\ell\\
\beta &=& 1-\alpha
\beta &=& 1-\alpha
Zeile 66: Zeile 67:


Zusätzlich zum "klassischen" Randwertproblem haben wir hier Zwangsbedingung durch die starren Stäbe und die Kontaktkraft.
Zusätzlich zum "klassischen" Randwertproblem haben wir hier Zwangsbedingung durch die starren Stäbe und die Kontaktkraft.
==tmp==
In dieser Lösung arbeiten wir mit dimensionslosen Koordinaten für die unabhängige Koordinate ''x'' und die abhängige Koordinate ''w(x)''.


<!-------------------------------------------------------------------------------->
<!-------------------------------------------------------------------------------->
{{MyCodeBlock|title=Header
{{MyCodeBlock|title=Header
|text=Text
|text=
In dieser Lösung arbeiten wir mit dimensionslosen Koordinaten für die unabhängige Koordinate ''x'' und die abhängige Koordinate ''w(x)''.
|code=
|code=
<syntaxhighlight lang="lisp" line start=1>
<syntaxhighlight lang="lisp" line start=1>
1+1
/*******************************************************/
/* MAXIMA script                                      */
/* version: wxMaxima 18.10.1                           */
/* author: Andreas Baumgart                            */
/* last updated: 2019-12-07                            */
/* ref: TMC, Labor 1                                   */
/* description: analytical solution for beam with      */
/*              with parallel rod guidance            */
/*******************************************************/
</syntaxhighlight>
</syntaxhighlight>
}}
}}
==tmp==
<!-------------------------------------------------------------------------------->
<!-------------------------------------------------------------------------------->
{{MyCodeBlock|title=Declarations
{{MyCodeBlock|title=Declarations
|text=
|text=
Zeile 98: Zeile 99:
|code=
|code=
<syntaxhighlight lang="lisp" line start=1>
<syntaxhighlight lang="lisp" line start=1>
1+1
declare("ℓ", alphabetic);
/* system parameters                                  */
params: [k[B]    = gamma*EI/ℓ^3,
        q[0]    = m*g/ℓ,
        K        = kappa*Q[ref],
        d        = delta*W[ref],
        gamma    = 5,
        kappa    = 3,
        delta    = 1];
 
geometry : [beta  = 1-alpha,
            alpha = 1/3];
 
/* reference deflection selected:                        */
/* beam with doublesided bearing under constant line load */
dimless : [M[ref]=Q[ref]*ℓ, Q[ref]= 3*EI/ℓ^3*W[ref],
          F = Q[ref]*f, B[z] = Q[ref]*b,
          W[ref] = 5*q[0]*ℓ^4/(EI*384)];
 
dimless: subst(params, dimless);
params : subst(dimless, params);
</syntaxhighlight>
</syntaxhighlight>
}}
}}


==tmp==
<!-------------------------------------------------------------------------------->
<!-------------------------------------------------------------------------------->


Zeile 145: Zeile 166:
|code=
|code=
<syntaxhighlight lang="lisp" line start=1>
<syntaxhighlight lang="lisp" line start=1>
1+1
/******************************************************/
/* Boundary Value Problem Formulation                */
/* part I: field                                      */
dgl : EI*W[ref]*diff(w(xi),xi,4)/ℓ^4 = q[0];
dgl: solve(subst(dimless,subst(params,dgl)),'diff(w(xi),xi,4))[1];
 
params : append(params, [mu = rhs(dgl)]);
dgl: lhs(dgl) = mu;
 
/* generic solution */
displ : solve(integrate(integrate(integrate(integrate(dgl,xi),xi),xi),xi),w(xi));
sections: [[i=1, %c4=C[1,0], %c3=C[1,1], %c2=C[1,2], %c1=C[1,3]],
          [i=2, %c4=C[2,0], %c3=C[2,1], %c2=C[2,2], %c1=C[2,3]]];
 
/* section I */
define(  w[1](xi),  ratsimp(subst(sections[1],W[ref]*subst(displ,w(xi)))));
define(Phi[1](xi),    diff(w[1](xi),xi  )/ℓ  );
define(  M[1](xi), -EI*diff(w[1](xi),xi,2)/ℓ^2);
define(  Q[1](xi), -EI*diff(w[1](xi),xi,3)/ℓ^3);
 
/* section II */
define(  w[2](xi),  ratsimp(subst(sections[2],subst(displ,W[ref]*w(xi)))));
define(Phi[2](xi),    diff(w[2](xi),xi  )/ℓ  );
define(  M[2](xi), -EI*diff(w[2](xi),xi,2)/ℓ^2);
define(  Q[2](xi), -EI*diff(w[2](xi),xi,3)/ℓ^3);
</syntaxhighlight>
</syntaxhighlight>
}}
}}
==tmp==


<!-------------------------------------------------------------------------------->
<!-------------------------------------------------------------------------------->
{{MyCodeBlock|title=Boundary Conditions
{{MyCodeBlock|title=Boundary Conditions
|text=
|text=
Zeile 229: Zeile 269:
|code=
|code=
<syntaxhighlight lang="lisp" line start=1>
<syntaxhighlight lang="lisp" line start=1>
1+1
/******************************************************/
/* Boundary Value Problem Formulation                */
/* part II: boundary conditions                      */
node[A]: [ w[1](0) = 0,
          M[1](0) = 0];
node[B]: [ w[1](alpha) = w[2](0),
          Phi[1](alpha) = Phi[2](0),     
          -Q[1](alpha) + Q[2](0) + B[z] = 0,
          -M[1](alpha) + M[2](0) = 0];
node[C]: [ w[2](beta) = d,
          -Q[2](beta) - K = 0,
          -M[2](beta) = 0];
node[D]: [F - k[B]*w[2](0)-B[z] = 0];
 
BCs : append(node[A],node[B],node[C],node[D]);
scale: [W[ref], -EI*W[ref]/(ℓ^2), 4*W[ref]/(24), 4*W[ref]/(24*ℓ),
        EI*W[ref]/(ℓ^3), EI*W[ref]/(ℓ^2), 4*W[ref]/24, EI*W[ref]/(ℓ^3),
        EI*W[ref]/(ℓ^2), 5*m*g/128];
BCs : makelist(expand(BCs[i]/scale[i]),i,1,length(BCs));
BCs : subst(dimless,subst(params,subst(dimless,subst(params,BCs))));
</syntaxhighlight>
</syntaxhighlight>
}}
}}


==tmp==
<!-------------------------------------------------------------------------------->
<!-------------------------------------------------------------------------------->
{{MyCodeBlock|title=Prepare for Solver
{{MyCodeBlock|title=Prepare for Solver
Zeile 284: Zeile 341:
{\tilde{B}_z}
{\tilde{B}_z}
\end{pmatrix} =  
\end{pmatrix} =  
\begin{pmatrix}
\begin{pmatrix}
0\\
0\\
Zeile 299: Zeile 354:
|code=
|code=
<syntaxhighlight lang="lisp" line start=1>
<syntaxhighlight lang="lisp" line start=1>
1+1
/* integration constants = unknowns */
ICs : [C[1,0],C[1,1],C[1,2],C[1,3],C[2,0],C[2,1],C[2,2],C[2,3],f,b];
ACM: augcoefmatrix(BCs,ICs)$
/* system matrix and rhs */
AA :  submatrix(ACM,10)$
bb : - col(ACM,10)$
/* print OLE */
print(subst(params,AA),"*",transpose(ICs),"=",subst(params,bb))$
</syntaxhighlight>
</syntaxhighlight>
}}
}}
==tmp==


<!-------------------------------------------------------------------------------->
<!-------------------------------------------------------------------------------->
{{MyCodeBlock|title=Solving
{{MyCodeBlock|title=Solving
|text=
|text=
Zeile 344: Zeile 402:
\end{array} \right)</math>.
\end{array} \right)</math>.


Die Schaltkraft ''F'' ist negativ! Die Gewichtskraft des Balkens ist für die gewählte Konstellation zu groß, um die Kontaktkraft ohne ''F'' einzustellen. Wir schauen uns an, wie die erforderliche Kraft F über α verläuft:
[[Datei:Kig1-21.png|mini|Ratio F/Q]]Die Schaltkraft ''F'' ist negativ! Die Gewichtskraft des Balkens ist für die gewählte Konstellation zu groß, um die Kontaktkraft ohne ''F'' einzustellen. Wir schauen uns an, wie die erforderliche Kraft F über α verläuft:
[[Datei:Kig1-21.png|mini|Ratio F/Q]]<br clear="all"/>
<br clear="all"/>
|code=
|code=
<syntaxhighlight lang="lisp" line start=1>
<syntaxhighlight lang="lisp" line start=1>
1+1
/* ... solve .....*/
sol : ratsimp(solve(BCs, ICs))[1]$
 
/* check for required force */
part: expand(subst(geometry[1],subst(sol,f)));
plot2d(part, [alpha,0,1], [y,-50,0], [xlabel, "α ->"], [ylabel, "F/Q[ref] ->"])$
</syntaxhighlight>
</syntaxhighlight>
}}
}}


==tmp==
<!-------------------------------------------------------------------------------->
<!-------------------------------------------------------------------------------->
{{MyCodeBlock|title=Post-Processing
{{MyCodeBlock|title=Post-Processing
|text=
|text=
Zeile 363: Zeile 424:


==== ... für ''ϕ(x)'': ====
==== ... für ''ϕ(x)'': ====
[[Datei:Kig1-23.png|mini|phi(x)]]<br clear="all"/>
[[Datei:Kig1-23.png|mini|''ϕ(x)''|ohne]]<br clear="all"/>


==== ... für ''M(x)'': ====
==== ... für ''M(x)'': ====
[[Datei:Kig1-24.png|mini|M(x)]]<br clear="all"/>
[[Datei:Kig1-24.png|mini|''M(x)''|ohne]]<br clear="all"/>


==== ... für ''Q(x)'': ====
==== ... für ''Q(x)'': ====
[[Datei:Kig1-25.png|mini|Q(x)]]<br clear="all"/>
[[Datei:Kig1-25.png|mini|''Q(x)''|ohne]]<br clear="all"/>


====... für die Lager-Reaktionskräfte:====
====... für die Lager-Reaktionskräfte:====
Zeile 377: Zeile 438:
|code=
|code=
<syntaxhighlight lang="lisp" line start=1>
<syntaxhighlight lang="lisp" line start=1>
1+1
/******************************************************/
/* Post - Processing                                  */
/* part II: boundary conditions                      */
 
reactForces: [A[z] = Q[1](0)];
expand(subst(geometry,subst(params,subst(dimless,subst(params,subst(sol, subst(sol, reactForces)))))));
/* plot displacements */
fcts: [[ w [1](xi), w [2](xi)],
      [Phi[1](xi),Phi[2](xi)],
      [ M [1](xi), M [2](xi)],
      [ Q [1](xi), Q [2](xi)]];
facts: [1/W[ref], 1/(W[ref]/ℓ), 1/M[ref], 1/Q[ref]];
textlabels : ["w(x)/W[ref]→", "w'(x)/(W[ref]/ℓ)→", "M(x)/M[ref]→", "Q(x)/Q[ref]→"];
for i: 1 thru 4 do(
  f : expand(subst(geometry,subst(dimless,subst(params,facts[i]*[subst(sol, fcts[i][1]),
                                                                subst(sol, fcts[i][2])])))),
  preamble: if i<=2 then "set yrange [] reverse" else "set yrange []",
  plot2d([[parametric,      t, subst(t,xi,f[1]), [t,0,subst(geometry,alpha)]],
          [parametric, subst(geometry,alpha)+t, subst(t,xi,f[2]), [t,0,subst(geometry,beta)]]],
                            [legend, "sec. I", "sec. II"],
                            [gnuplot_preamble, preamble],
                            [xlabel, "x/ℓ ->"],
                            [ylabel, textlabels[i]]))$
</syntaxhighlight>
</syntaxhighlight>
}}
}}

Aktuelle Version vom 14. September 2022, 07:33 Uhr


Aufgabenstellung

Das skizzierte System ist ein elektrischer Schalter, der in C einen Stromkreis schließen soll. Es besteht aus einem durchgehenden elastischen Stab ABC (Euler-Bernoulli-Balken: Masse m, Biegesteifigkeit EI, Länge ), der in A gelenkig gelagert ist, in B mit einer Führung aus drei starren, masselosen Stäben gelenkig verbunden ist und dem Kontakt in C. Die Führung ist durch eine Feder (Steifigkeit k) mit der Umgebung verbunden.


Lageplan

Damit der Schalter zuverlässig funktioniert soll die Kraft F an der Führung so gewählt werden, dass die vorgegebene Kontaktkraft K eingestellt wird. Gesucht ist die analytische Lösung für ein Euler-Bernoulli-Modell der elastischen Struktur.

Wir lösen dazu das Gleichungssystem zur Bestimmung der Integrationskonstanten der Differentialgleichung des Euler-Bernoulli-Balkens und weiterer Unbekannten und geben die grafischen Lösungen für w, w', M und Q an.

Gegeben: K, ℓ, EI, m, g, k

Lösung mit Maxima

Koordinaten

Das System besteht aus den drei starren Stäben und dem elastischen Balken.

Die Verschiebungen und Verbiegungen der Stäbe verstehen wir am besten, wenn wir uns das System einmal im ausgelenkten Zustand ansehen. Gleichgewichtsbedingungen für die Führung erhalten wir aus einfachen Kräfte- und Momentenbilanzen, die Gleichgewichtsbedingungen für den Euler-Bernoulli-Balken ist die Feld-Differentialgleichung, deren Integrationskonstanten wir an die Rand- und Übergangsbedingungen anpassen.

Aus dem Bild sehen wir

,

und wir arbeiten im folgenden mit

.

Wir beginnen, indem wir Führung und Balken von einander freischneiden:

    

Damit legen wir die Schnittkräfte

  • Bz und
  • K

frei. K ist dabei gegeben, F und Bz sind unbekannt.

Wir müssen überlegen, wo wir Bedingungen für diese beiden Größen herbekommen.

Die Lösung der Teilaufgabe für den Euler-Bernoulli-Balken ist ein klassisches Randwertproblem mit

  1. zwei Gebieten, in denen ein Euler-Bernoulli-Balken in AB (Bereich I) und BC (Bereich II) durch eine Streckenlast q0 belastet ist sowie
  2. Rand- und Übergangsbedingungen in den Punkten A, B, C.

Die Biegesteifigkeit des Balkens ist konstant (also nicht von "x" abhängig), die Lösungen in beiden Bereichen ist also bis auf die Integrationskonstanten gleich.

Wir verwenden ein xi bzw. ξi als Ortskoordinaten je Bereich.

Zusätzlich zum "klassischen" Randwertproblem haben wir hier Zwangsbedingung durch die starren Stäbe und die Kontaktkraft.

Header

In dieser Lösung arbeiten wir mit dimensionslosen Koordinaten für die unabhängige Koordinate x und die abhängige Koordinate w(x).


/*******************************************************/
/* MAXIMA script                                       */
/* version: wxMaxima 18.10.1                           */
/* author: Andreas Baumgart                            */
/* last updated: 2019-12-07                            */
/* ref: TMC, Labor 1                                   */
/* description: analytical solution for beam with      */
/*              with parallel rod guidance             */
/*******************************************************/



Declarations

Die Streckenlast auf den Balken ist seine Gewichtskraft, also

.

Später werden wir noch die dimensionslose Koordinate

gebrauchen.


declare("ℓ", alphabetic);
 
/* system parameters                                  */
params: [k[B]     = gamma*EI/ℓ^3,
         q[0]     = m*g/ℓ,
         K        = kappa*Q[ref],
         d        = delta*W[ref],
         gamma    = 5,
         kappa    = 3,
         delta    = 1];

geometry : [beta  = 1-alpha,
            alpha = 1/3];

/* reference deflection selected:                         */
/* beam with doublesided bearing under constant line load */
dimless : [M[ref]=Q[ref]*ℓ, Q[ref]= 3*EI/ℓ^3*W[ref],
           F = Q[ref]*f, B[z] = Q[ref]*b,
           W[ref] = 5*q[0]*ℓ^4/(EI*384)]; 

dimless: subst(params, dimless);
params : subst(dimless, params);




Integration Of Differential Equation

In Bereich I und II gilt dieselbe Bewegungs-Differentialgleichung

,

die wir durch Integration lösen und dann bereichsweise an Rand- und Übergangsbedingungen anpassen. Die Bedeutung der gesuchten Auslenkung w und seiner Ableitungen sind

Einfacher wird es, wenn wir mit der dimensionslosen Koordinate

.

arbeiten. Dann ist

Praktisch ist es nun, Wref so zu wählen, dass wir die dimensionslose Koordinate für w einfach interpretieren können. Aus der Musterlösung für den beidseitig gelenkig gelagerten Euler-Bernoulli-Balken unter konstanter Streckenlast nehmen wir

und damit

.

Ein Aufintegrieren der Differentialgleichung liefert dann


/******************************************************/
/* Boundary Value Problem Formulation                 */
/* part I: field                                      */
dgl : EI*W[ref]*diff(w(xi),xi,4)/ℓ^4 = q[0];
dgl: solve(subst(dimless,subst(params,dgl)),'diff(w(xi),xi,4))[1];

params : append(params, [mu = rhs(dgl)]);
dgl: lhs(dgl) = mu;

/* generic solution */
displ : solve(integrate(integrate(integrate(integrate(dgl,xi),xi),xi),xi),w(xi));
sections: [[i=1, %c4=C[1,0], %c3=C[1,1], %c2=C[1,2], %c1=C[1,3]],
           [i=2, %c4=C[2,0], %c3=C[2,1], %c2=C[2,2], %c1=C[2,3]]];
  
/* section I */
define(  w[1](xi),  ratsimp(subst(sections[1],W[ref]*subst(displ,w(xi)))));
define(Phi[1](xi),     diff(w[1](xi),xi  )/ℓ  );
define(  M[1](xi), -EI*diff(w[1](xi),xi,2)/ℓ^2);
define(  Q[1](xi), -EI*diff(w[1](xi),xi,3)/ℓ^3);
  
/* section II */
define(  w[2](xi),  ratsimp(subst(sections[2],subst(displ,W[ref]*w(xi)))));
define(Phi[2](xi),     diff(w[2](xi),xi  )/ℓ  );
define(  M[2](xi), -EI*diff(w[2](xi),xi,2)/ℓ^2);
define(  Q[2](xi), -EI*diff(w[2](xi),xi,3)/ℓ^3);




Boundary Conditions

Für die 2*4 = 8 Integrationskonstanten

und die weiteren gesuchten Größen

suchen wir jetzt die passenden Gleichungen aus Rand- und Übergangsbedingungen des Systems.

Die Normalkräfte N brauchen wir dabei nicht auszuwerten.

Rand A.

Aus Rand A:

Geometrische Randbedingungen

Kraft- und Momenten-Randbedingungen

keine

Aus Übergang B:

Übergang B

Geometrische Randbedingungen

Kraft- und Momenten-Randbedingungen

  1. .

Aus Rand C:

Rand C

Geometrische Randbedingungen

Kraft- und Momenten-Randbedingungen

Gleichgewicht an der Führung

Führung

Für die starren Stäbe stellen wir die Summe aller Kräfte im senkrechten Stab auf:

.

Ein bisschen exotisch ist, dass wir nun zehn Unbekannte haben, nämlich

Aber mit den Randbedingungen oben steht uns ein vollständiges Gleichungssystem für diese zehn Unbekannten - die Integrationskonstanten und die Kräfte F, Bz - zur Verfügung.

Auch hier arbeiten wir mit dimensionslosen Größen, hier

,

wobei wir

setzen.


/******************************************************/
/* Boundary Value Problem Formulation                 */
/* part II: boundary conditions                       */
node[A]: [ w[1](0) = 0,
           M[1](0) = 0];
node[B]: [ w[1](alpha) = w[2](0),
           Phi[1](alpha) = Phi[2](0),       
          -Q[1](alpha) + Q[2](0) + B[z] = 0,
          -M[1](alpha) + M[2](0) = 0];
node[C]: [ w[2](beta) = d,
          -Q[2](beta) - K = 0,
          -M[2](beta) = 0];
node[D]: [F - k[B]*w[2](0)-B[z] = 0];

BCs : append(node[A],node[B],node[C],node[D]);
scale: [W[ref], -EI*W[ref]/(ℓ^2), 4*W[ref]/(24), 4*W[ref]/(24*ℓ),
        EI*W[ref]/(ℓ^3), EI*W[ref]/(ℓ^2), 4*W[ref]/24, EI*W[ref]/(ℓ^3),
        EI*W[ref]/(ℓ^2), 5*m*g/128];
BCs : makelist(expand(BCs[i]/scale[i]),i,1,length(BCs));
BCs : subst(dimless,subst(params,subst(dimless,subst(params,BCs))));




Prepare for Solver

Das Gleichungssystem wollen wir als

schreiben. Um die Ergebnisse kompakt darstellen zu können, wählen wir als Parameter

mit den dimensionslosen Größen

.

Damit erhalten wir - hier nach Einsetzen der gegeben System-Parameter:

.

/* integration constants = unknowns */
ICs : [C[1,0],C[1,1],C[1,2],C[1,3],C[2,0],C[2,1],C[2,2],C[2,3],f,b];
ACM: augcoefmatrix(BCs,ICs)$
/* system matrix and rhs */
AA :   submatrix(ACM,10)$
bb : - col(ACM,10)$
/* print OLE */
print(subst(params,AA),"*",transpose(ICs),"=",subst(params,bb))$




Solving

Das Lösen des Gleichungssystems liefert dann

und für

die dimensionslosen Werte

.
Ratio F/Q

Die Schaltkraft F ist negativ! Die Gewichtskraft des Balkens ist für die gewählte Konstellation zu groß, um die Kontaktkraft ohne F einzustellen. Wir schauen uns an, wie die erforderliche Kraft F über α verläuft:



/* ... solve .....*/
sol : ratsimp(solve(BCs, ICs))[1]$

/* check for required force */
part: expand(subst(geometry[1],subst(sol,f)));
plot2d(part, [alpha,0,1], [y,-50,0], [xlabel, "α ->"], [ylabel, "F/Q[ref] ->"])$




Post-Processing

Die Grafen der gesuchten Funktionen tragen wir nun direkt dimensionslos auf:

... für w(x):

w(x)


... für ϕ(x):

ϕ(x)


... für M(x):

M(x)


... für Q(x):

Q(x)


... für die Lager-Reaktionskräfte:


/******************************************************/
/* Post - Processing                                  */
/* part II: boundary conditions                       */

reactForces: [A[z] = Q[1](0)];
expand(subst(geometry,subst(params,subst(dimless,subst(params,subst(sol, subst(sol, reactForces)))))));
 
/* plot displacements */
 
fcts: [[ w [1](xi), w [2](xi)],
       [Phi[1](xi),Phi[2](xi)],
       [ M [1](xi), M [2](xi)],
       [ Q [1](xi), Q [2](xi)]];
facts: [1/W[ref], 1/(W[ref]/ℓ), 1/M[ref], 1/Q[ref]];
 
textlabels : ["w(x)/W[ref]→", "w'(x)/(W[ref]/ℓ)→", "M(x)/M[ref]→", "Q(x)/Q[ref]→"];
for i: 1 thru 4 do(
  f : expand(subst(geometry,subst(dimless,subst(params,facts[i]*[subst(sol, fcts[i][1]),
                                                                 subst(sol, fcts[i][2])])))),
  preamble: if i<=2 then "set yrange [] reverse" else "set yrange []",
  plot2d([[parametric,       t, subst(t,xi,f[1]), [t,0,subst(geometry,alpha)]],
          [parametric, subst(geometry,alpha)+t, subst(t,xi,f[2]), [t,0,subst(geometry,beta)]]],
                             [legend, "sec. I", "sec. II"],
                             [gnuplot_preamble, preamble],
                             [xlabel, "x/ℓ ->"],
                             [ylabel, textlabels[i]]))$






Links

  • ...

Literature

  • ...