KerMor  0.9
Model order reduction for nonlinear dynamical systems and nonlinear approximation
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
CurvedBeam.m
Go to the documentation of this file.
1 namespace models{
2 namespace beam{
3 
4 
5 /* (Autoinserted by mtoc++)
6  * This source code has been filtered by the mtoc++ executable,
7  * which generates code that can be processed by the doxygen documentation tool.
8  *
9  * On the other hand, it can neither be interpreted by MATLAB, nor can it be compiled with a C++ compiler.
10  * Except for the comments, the function bodies of your M-file functions are untouched.
11  * Consequently, the FILTER_SOURCE_FILES doxygen switch (default in our Doxyfile.template) will produce
12  * attached source files that are highly readable by humans.
13  *
14  * Additionally, links in the doxygen generated documentation to the source code of functions and class members refer to
15  * the correct locations in the source code browser.
16  * However, the line numbers most likely do not correspond to the line numbers in the original MATLAB source files.
17  */
18 
20  :public models.beam.Beam {
39  public:
40 
41  pc;
42 
43  R;
44 
46 
48 
49  B;
50 
51  B3;
52 
53  B4;
54 
55 
56  public:
57 
59 
61 
62 
63  private:
64 
65  TM;
66 
67 
68  public:
69 
70 
71  CurvedBeam(models.BaseFullModel model,material,pointsidx) {
72  this = this@models.beam.Beam(model, material, pointsidx(1:2));
73  this.pc= pointsidx(3);
74  this.initialize;
75  }
76 
77 
78  function M = getLocalMassMatrix() {
79 
80  c = this.c;
81 
82  /* Ansatzfunktionen N und Verbindungsmatrix für einen Viertelkreis holen */
83 
84  /* Matrizen aufsetzen */
85  rhoJ = blkdiag(c(7), c(6), c(6));
86 
87  M = zeros(12);
88 
89  L = this.Length;
90  /* Stützstellen für Gaußquadratur */
91  s = 0.5*L * [ (1-sqrt(3/5)), 1, (1+sqrt(3/5)) ];
92  /* Gewichte für Gaußquadratur */
93  w = 0.5*L * 1/9 * [5 8 5];
94 
95  for i=1:3
96  N = this.circle_shape_functions(s(i), this.B);
97  N1 = N(1:3,:);
98  N2 = N(4:6,:);
99  M = M + w(i) * ( c(5)*(N1" *N1) + N2 "*rhoJ*N2 );
100  end
101 
102  M = this.TG * M * this.TG^t;
103  }
120  function K = getLocalStiffnessMatrix() {
121 
122  c = this.c;
123  /* Ansatzfunktionen N und Verbindungsmatrix für einen Viertelkreis holen */
124 
125  /* Matrizen aufsetzen */
126  D = blkdiag(c(3), c(4), c(4));
127  E = blkdiag(c(2), c(1), c(1));
128  rhoJ = blkdiag(c(7), c(6), c(6));
129 
130  M = zeros(12);
131 
132  L = this.Length;
133  /* Stützstellen für Gaußquadratur */
134  s = 0.5*L * [ (1-sqrt(3/5)), 1, (1+sqrt(3/5)) ];
135  /* Gewichte für Gaußquadratur */
136  w = 0.5*L * 1/9 * [5 8 5];
137 
138  for i=1:3
139  N = this.circle_shape_functions(s(i), this.B);
140  N1 = N(1:3,:);
141  N2 = N(4:6,:);
142  M = M + w(i) * ( c(5)*N1" *N1 + N2 "*rhoJ*N2 );
143  end
144 
145  B3 = this.B(7:9,:);
146  B4 = this.B(10:12,:);
147 
148  K = L * (B4" * E * B4 + B3 " * D * B3);
149 
150  K = this.TG * K * this.TG^t;
151  }
168  function f = getLocalForceMatrix() {
169 
170  f = zeros(12,3);
171 
172  L = this.Length;
173  q_lok = (this.c(5) + this.Material.q_plus) * this.T^t;
174  /* Stützstellen für Gaußquadratur */
175  s = 0.5*L * [(1-sqrt(3/5)), 1, (1+sqrt(3/5))];
176  /* Gewichte für Gaußquadratur */
177  w = 0.5*L * 1/9 * [5 8 5];
178  for i=1:3
179  N = this.circle_shape_functions(s(i), this.B);
180  f = f + w(i) * N(1:3,:)" *(this.Fren(s(i)/this.R) "*q_lok);
181  end
182 
183  f = this.TG * f;
184  }
201  function [K , R , U_pot ] = getLocalTangentials(u) {
202 
203  /* Wertet tangentielle Steifigkeitsmatrix, "Residuumsvektor" (= Vektor der virtuellen internen Energie) und potenzielle Energie eines
204  * Timoshenko-Balkens (Kreisbogen mit Öffnungswinkel [angle]) (numerisch mit speziell gewonnenen Ansatzfunktionen) aus */
205 
206  /* c1 = E*I
207  * c2 = G*It
208  * c3 = E*A
209  * c4 = G*As
210  * c5 = rho*A
211  * c6 = rho*I
212  * c7 = rho*It */
213 
214  c = this.c;
215  B = this.B;
216 
217  /* Potenzielle Energie */
218  U_pot = 0;
219  /* Residuumsvektor */
220  R = zeros(12,1);
221  /* Tangentielle lokale Steifigkeitsmatrix */
222  K = zeros(12);
223 
224  /* lokale Verschiebungen des Elements */
225  u_lok = this.TG^t * u;
226 
227  /* Stoff-Matrizen aufsetzen */
228  D = blkdiag(c(3), c(4), c(4));
229  E = blkdiag(c(2), c(1), c(1));
230 
231  L = this.Length;
232  /* % Stützstellen und Gewichte für Gaußquadratur */
233  num_gauss_points = 3;
234  switch (num_gauss_points)
235  case 1
236  /* Exakt bis Grad 1 */
237  s = 0.5*L;
238  w = L;
239  case 2
240  /* Exakt bis Grad 3 */
241  s = 0.5*L * [ (1-sqrt(1/3)), (1+sqrt(1/3)) ];
242  w = 0.5*L * [1 1];
243  case 3
244  /* Exakt bis Grad 5 */
245  s = 0.5*L * [ (1-sqrt(3/5)), 1, (1+sqrt(3/5)) ];
246  w = 0.5*L * 1/9 * [5 8 5];
247  case 4
248  /* Exakt bis Grad 7 */
249  s = 0.5*L * [ ( 1 - sqrt( 3/7 + 2/7*sqrt(6/5) ) ),
250  ( 1 - sqrt( 3/7 - 2/7*sqrt(6/5) ) ),
251  ( 1 + sqrt( 3/7 - 2/7*sqrt(6/5) ) ),
252  ( 1 + sqrt( 3/7 + 2/7*sqrt(6/5) ) )];
253  w = 0.5*L * 1/36 * [18-sqrt(30), 18+sqrt(30), 18+sqrt(30), 18-sqrt(30)];
254  case 5
255  /* Exakt bis Grad 9 */
256  s = 0.5*L * [ 1 - 1/3 * sqrt( 5 + 2*sqrt(10/7) ),
257  1 - 1/3 * sqrt( 5 - 2*sqrt(10/7) ),
258  1,
259  1 + 1/3 * sqrt( 5 - 2*sqrt(10/7) ),
260  1 + 1/3 * sqrt( 5 + 2*sqrt(10/7) )];
261  w = 0.5*L * 1/900 * [322 - 13*sqrt(70), 322 + 13*sqrt(70), 512, 322 + 13*sqrt(70), 322 - 13*sqrt(70)];
262  otherwise
263  /* Exakt bis Grad 5 */
264  s = 0.5*L * [ (1-sqrt(3/5)), 1, (1+sqrt(3/5)) ];
265  w = 0.5*L * 1/9 * [5 8 5];
266  end
267 
268  /* %
269  * Kreuzproduktsmatrix e1 x v = S_e1 * v */
270  S_e1 = [0 0 0; 0 0 -1; 0 1 0];
271  /* Matrizen, die die lokale Verschiebung auf ihre (über den Kreisbogen
272  * konstante!) *LINEARE* Verzerrung und Krümmung abbilden */
273  B3 = B(7:9,:);
274  B4 = B(10:12,:);
275 
276  for i=1:length(s)
277  /* Auswerten der Ansatzfunktionen */
278  N = this.circle_shape_functions(s(i), B);
279  /* N1 = N(1:3,:); % Verschiebungsansätze */
280  N2 = N(4:6,:); /* Verdrehungsansätze */
281 
282 
283  /* Vorberechnungen zur Auswertung der Verzerrungen */
284  E_lin = B3; /* = (für gerade Balken) N1_prime + S_e1 * N2; */
285 
286  E_skw = E_lin - 2 * S_e1 * N2; /* = N1_prime - S_e1 * N2; */
287 
288 
289  E_nl_1 = 0.125 * ( E_skw(2,:)" * E_skw(2,:) + E_skw(3,:) " * E_skw(3,:) );
290  E_nl_2 = 0.5 * N2(1,:)^t * E_skw(3,:);
291  E_nl_3 = 0.5 * N2(1,:)^t * E_skw(2,:);
292 
293  E_1 = 2 * E_nl_1; /* = E_nl_1 + E_nl_1'; */
294 
295  E_2 = E_nl_2 + E_nl_2^t;
296  E_3 = E_nl_3 + E_nl_3^t;
297 
298  E_nonlin = [u_lok" * E_nl_1; u_lok " * E_nl_2; u_lok^t * E_nl_3];
299  dE_nonlin = [u_lok" * E_1; u_lok " * E_2; u_lok^t * E_3];
300  dE = (E_lin + dE_nonlin);
301 
302  /* Verzerrung */
303  eps = (E_lin + E_nonlin) * u_lok;
304  /* Krümmung */
305  kap = B4 * u_lok; /* = N2_prime * u_lok; */
306 
307 
308  /* Kraft */
309  F = D * eps;
310  /* Moment */
311  M = E * kap;
312 
313  /* Integrations-Update potenzielle Energie */
314  U_pot = U_pot + w(i) * ( eps" * F + kap " * M);
315  /* Integration */
316  R = R + w(i) * ( dE" * F + B4 " * M);
317  K = K + w(i) * ( dE" * D * dE + B4 " * E * B4 + F(1)*E_1 + F(2)*E_2 + F(3)*E_3 );
318 
319  end
320 
321  /* Transformation in globales Koordinatensystem (nicht nötig bei Skalaren) */
322  K = this.TG * K * this.TG^t;
323  R = this.TG * R;
324  }
325 
326 
327  function B = circle_connect_matrix() {
328  C0 = this.circle_shape_functions(0, eye(12));
329  /* Ehemals: L = this.R * this.angle */
330  CL = this.circle_shape_functions(this.Length, eye(12));
331  B = inv([C0; CL]);
332  }
333 
334 
335  function N = circle_shape_functions(s,B) {
336 
337  /* => x(s) = C(s)*B * v =: N(s) * [x(0); x(L)] ! */
338 
339  /* Da M für jedes Element konstant ist (hängt nur von L ab!), könnte M im
340  * Voraus berechnet werden. */
341 
342  /* Anmerkung: Die Basisfunktionen sind so gebaut, dass die Ableitung
343  * (d/ds u = u' + Gu, etc) eben jene Konstanten sind, die in der zweiten
344  * Hälfte von c stehen (c2 s.o.): d/ds N(s)*v = (B*v)(7:12)
345  * (d/ds x(s) = d/ds (C(s))*c = c2 ?) */
346 
347 
348  /* C = @(s) [
349  * cos(s/R) sin(s/R) 0 0 0 R-R*cos(s/R) R*sin(s/R) R-R*cos(s/R) 0 0 0 R*s-R^2*sin(s/R);
350  * -sin(s/R) cos(s/R) 0 0 0 R*sin(s/R) R*cos(s/R)-R R*sin(s/R) 0 0 0 R^2-R^2*cos(s/R);
351  * 0 0 1 R-R*cos(s/R) -R*sin(s/R) 0 0 0 s R*s-R^2*sin(s/R) R^2*cos(s/R)-R^2 0;
352  * 0 0 0 cos(s/R) sin(s/R) 0 0 0 0 R*sin(s/R) R-R*cos(s/R) 0;
353  * 0 0 0 -sin(s/R) cos(s/R) 0 0 0 0 R*cos(s/R)-R R*sin(s/R) 0;
354  * 0 0 0 0 0 1 0 0 0 0 0 s]; */
355  co = cos(s/this.R);
356  si = sin(s/this.R);
357  RmRco = this.R-this.R*co;
358  R2mR2co = this.R*RmRco;
359  Rsi = this.R*si;
360  RsmR2si = this.R*s - this.R^2*si;
361  Cs = [
362  co si 0 0 0 RmRco Rsi RmRco 0 0 0 RsmR2si;
363  -si co 0 0 0 Rsi -RmRco Rsi 0 0 0 R2mR2co;
364  0 0 1 RmRco -Rsi 0 0 0 s RsmR2si -R2mR2co 0;
365  0 0 0 co si 0 0 0 0 Rsi RmRco 0;
366  0 0 0 -si co 0 0 0 0 -RmRco Rsi 0;
367  0 0 0 0 0 1 0 0 0 0 0 s];
368  /* Cs = [
369  * co si 0 0 0 R-R*co R*si R-R*co 0 0 0 R*s-R^2*si;
370  * -si co 0 0 0 R*si R*co-R R*si 0 0 0 R^2-R^2*co;
371  * 0 0 1 R-R*co -R*si 0 0 0 s R*s-R^2*si R^2*co-R^2 0;
372  * 0 0 0 co si 0 0 0 0 R*si R-R*co 0;
373  * 0 0 0 -si co 0 0 0 0 R*co-R R*si 0;
374  * 0 0 0 0 0 1 0 0 0 0 0 s]; */
375  N = Cs*B;
376  }
394  function plot(p,u1,u2,col1,col2,plot_options) {
395 
396 
397  L = this.Length;
398  R = this.R;
399  T = this.T;
400  N = this.split;
401  T1 = this.T_block1;
402  T2 = this.T_block2;
403  T_Fren = this.Fren;
404  B = this.B;
405  pc = p(this.pc,:); /* #ok<*PROP> */
406 
407 
408  /* Auswertungsstellen */
409  x = (0:L/(N+1):L)^t;
410 
411  /* Ruhelage plotten */
412  xx = R * cos(x/R);
413  yy = R * sin(x/R);
414  zz = 0*x;
415  COR = T * [xx" ; yy "; zz^t];
416  /* plot3( pc(1) + COR(1,:), pc(2) + COR(2,:), pc(3) + COR(3,:), 'k:' ); */
417 
418  u_nodes = [u1;u2];
419  u_glob = zeros(6,N+2);
420  u_lok = zeros(6,N+2);
421 
422  u_glob(:,1) = [T1 * u1(1:3); T1 * u1(4:6)];
423  u_glob(:,N+2) = [T2 * u2(1:3); T2 * u2(4:6)];
424  u_lok(:,1) = u1;
425  u_lok(:,N+2) = u2;
426 
427  for i = 2:N+1
428  N_U = this.circle_shape_functions(x(i), this.B);
429  u_lok(:,i) = N_U * u_nodes;
430  T_i = T_Fren(x(i)/R);
431  u_glob(:,i) = [T * T_i * u_lok(1:3,i); T * T_i * u_lok(4:6,i)];
432  end
433 
434  cmap = colormap;
435  col = round(x/L*col2 + (L-x)/L*col1);
436 
437  /* % Mittellinie plotten */
438  if (plot_options.centerline)
439  /* Einfarbig, dünn
440  * plot3(pc(1) + COR(1,:) + u_glob(1,:), pc(2) + COR(2,:) + u_glob(2,:), pc(3) + COR(3,:) + u_glob(3,:), 'Color', 0.5 * (cmap(col(1),:)+cmap(col(N+2),:)))
441  * elseif (centerline == 2)
442  * Mehrfarbig, dick */
443  for i=1:N+1
444  plot3(pc(1) + COR(1,[i i+1]) + u_glob(1,[i i+1]), pc(2) + COR(2,[i i+1]) + u_glob(2,[i i+1]), pc(3) + COR(3,[i i+1]) + u_glob(3,[i i+1]), " LineWidth ", plot_options.centerline, " Color ", 0.5 * (cmap(col(i),:)+cmap(col(i+1),:)) )
445  end
446  end
447 
448  if (plot_options.endmarker)
449  /* Elementgrenzenmarkierungen plotten */
450  plot3( pc(1) + COR(1,1) + u_glob(1,1), pc(2) + COR(2,1) + u_glob(2,1), pc(3) + COR(3,1) + u_glob(3,1), " + ", " LineWidth ", plot_options.endmarker, " Color ", cmap(col(1),:) )
451  plot3( pc(1) + COR(1,N+2) + u_glob(1,N+2), pc(2) + COR(2,N+2) + u_glob(2,N+2), pc(3) + COR(3,N+2) + u_glob(3,N+2), " + ", " LineWidth ", plot_options.endmarker, " Color ", cmap(col(N+2),:) )
452  end
453 
454  /* % Querschnitte plotten */
455  if (plot_options.crosssection ~= 0)
456  N_angle = 18;
457  R_cross = 0.6285;
458  split_angle = 0 : (2*pi/N_angle) : 2*pi;
459 
460  xx = zeros(1, N_angle+1);
461  yy = R_cross * cos(split_angle);
462  zz = R_cross * sin(split_angle);
463 
464  COR_lok = [xx; yy; zz];
465 
466  for i = 1:N+2
467  phi = u_lok(4,i);
468  psi = u_lok(5,i);
469  theta = u_lok(6,i);
470  rot_angle = norm([phi psi theta]);
471  S_rot = [0 -theta psi; theta 0 -phi; -psi phi 0];
472 
473  if (plot_options.crosssection == 1)
474  /* Rotationen 1. Ordnung */
475  ROT = eye(3) + S_rot;
476  elseif (plot_options.crosssection == 2)
477  /* Rotationen 2. Ordnung */
478  ROT = eye(3) + S_rot + 0.5*S_rot*S_rot;
479  else
480  /* Exakte Rotationen (R = e^(S_rot)) */
481  if (rot_angle > 1e-9)
482  ROT = eye(3) + sin(rot_angle)/rot_angle * S_rot + (1-cos(rot_angle))/rot_angle^2 * S_rot*S_rot;
483  else
484  ROT = eye(3);
485  end
486  end
487 
488  COR_u = T * T_Fren(x(i)/R) * ROT*COR_lok;
489 
490  col_act = cmap(col(i),:);
491  plot3(pc(1) + COR(1,i) + u_glob(1,i) + COR_u(1,:), pc(2) + COR(2,i) + u_glob(2,i)+ COR_u(2,:), pc(3) + COR(3,i) + u_glob(3,i) + COR_u(3,:), " Color ", col_act);
492  end
493  end
494  }
510  protected:
511 
512  function initialize() {
513  initialize@models.beam.Beam(this);
514 
515  m = this.Model;
516  /* Lokales Koordinatensystem in globalem entwickeln */
517  e_x = ( m.Points(this.PointsIdx(1),:) - m.Points(this.pc,:) )^t;
518  this.R= norm(e_x);
519  e_x = e_x / this.R;
520  e_y = (m.Points(this.PointsIdx(2),:) - m.Points(this.pc,:))^t / this.R;
521 
522  /* Öffnungswinkel des Kreissegments */
523  this.angle= acos(e_x^t*e_y);
524 
525  /* Länge */
526  this.Length= this.R*this.angle;
527 
528  e_y = e_y - (e_x^t*e_y) * e_x;
529  e_y = e_y / norm(e_y);
530 
531  e_z = [e_x(2)*e_y(3) - e_x(3)*e_y(2);
532  e_x(3)*e_y(1) - e_x(1)*e_y(3);
533  e_x(1)*e_y(2) - e_x(2)*e_y(1)];
534  e_z = e_z / norm(e_z);
535  this.T= [e_x e_y e_z];
536 
537  /* Frenet-Basis (in _lokalen_ Koords) als anonyme Funktion */
538  this.Fren= @(s) [-sin(s) -cos(s) 0; cos(s) -sin(s) 0; 0 0 1];
539 
540  /* Globale transformationsmatrix berechnen
541  * Transformationsmatrix: natürliche Koords -> glob. Koords
542  * Sortierung der Variablen:
543  * u0, v0, w0, phi0, psi0, theta0, u1, v1, w1, phi1, psi1, theta1
544  * 1 2 3 4 5 6 7 8 9 10 11 12 */
545  this.T_block1= this.T * this.Fren(0);
546  this.T_block2= this.T * this.Fren(this.angle);
547  TG = zeros(12);
548  TG([1 2 3], [1 2 3]) = this.T_block1; /* Verschiebung Anfangsknoten */
549 
550  TG([7 8 9], [7 8 9]) = this.T_block2; /* Verschiebung Endknoten */
551 
552  TG([4 5 6], [4 5 6]) = this.T_block1; /* Winkel Anfangsknoten */
553 
554  TG([10 11 12], [10 11 12]) = this.T_block2; /* Winkel Endknoten */
555 
556  this.TG= TG;
557 
558  /* Für die Auswertung der Ansatzfunktionen und Umrechnung der Knotengrößen in Verzerrungen */
559  B = this.circle_connect_matrix;
560  this.B= B;
561  this.B3= B(7:9,:);
562  this.B4= B(10:12,:);
563 
564  /* Effektive Konstanten */
565  m = this.Material;
566  /* Berechnung des Flexibilitätsfaktors (verringerte Biegesteifigkeit auf Grund von Ovalisierung) */
567  dm_tmp = m.d_a - m.s;
568  h_tmp = 4 * this.R * m.s / dm_tmp^2;
569  flex_factor = 1.65 / ( h_tmp*(1 + 6*m.p/m.E*(0.5*dm_tmp/m.s)^2*(this.R/m.s)^(1/3)) );
570  if (flex_factor < 1)
571  flex_factor = 1;
572  end
573 
574  this.c(1) = m.E * m.Iy / flex_factor; /* c1 = E*I */
575 
576  this.c(2) = m.G * m.It; /* c2 = G*It */
577 
578  this.c(3) = m.E * m.A; /* c3 = E*A */
579 
580  this.c(4) = m.k * m.G * m.A; /* c4 = G*As */
581 
582  this.c(5) = m.rho * m.A; /* c5 = rho*A */
583 
584  this.c(6) = m.rho * m.Iy; /* c6 = rho*I */
585 
586  this.c(7) = m.rho * m.It; /* c7 = rho*It */
587 
588  this.c(8) = m.rho; /* c8 = rho */
589 
590  }
591 
592 
593 };
594 }
595 }
596 
Model
The model that contains the structure element.
function plot(p, u1, u2, col1, col2, plot_options)
Zeichnet Timoshenko Kreisbogen Die Ansatzfunktionen werden dabei an N Zwischenstellen ausgewertet...
Definition: CurvedBeam.m:394
TG
c_theta = []; kappa = []; alphaA = []; The global transformation matrix (?)
function M = getLocalMassMatrix()
Berechnet lokale Steifigkeits- und Massenmatrix eines gekrümmten Timoshenko-Balkens (Viertelkreis) (nu...
Definition: CurvedBeam.m:78
CurvedBeam(models.BaseFullModel model, material, pointsidx)
Definition: CurvedBeam.m:71
The base class for any KerMor detailed model.
Definition: BaseFullModel.m:18
PointsIdx
Punkt-Index-Array.
Beam:
Definition: Beam.m:19
function K = getLocalStiffnessMatrix()
Berechnet lokale Steifigkeits- und Massenmatrix eines gekrümmten Timoshenko-Balkens (Viertelkreis) (nu...
Definition: CurvedBeam.m:120
Length
Indizes (lokal pro Knoten) in die die lokalen Matrizen assembliert werden Konvention: Freiheitsgrade ...
#define F(x, y, z)
Definition: CalcMD5.c:168
c
Effektive Stoffkonstanten.
function N = circle_shape_functions(s, B)
Wertet die Basisfunktionen für ein Kreiselement aus: Auf dem Element werden konstante Strains angenomm...
Definition: CurvedBeam.m:335
Material:
Definition: Material.m:19
function B = circle_connect_matrix()
Definition: CurvedBeam.m:327
function [ K , R , U_pot ] = getLocalTangentials(u)
Definition: CurvedBeam.m:201
split
Speichert den anteil des jeweiligen Elements an der Gesamtsystemlänge.
Definition: Beam.m:42
function initialize()
Definition: CurvedBeam.m:512
function f = getLocalForceMatrix()
Berechnet lokale Steifigkeits- und Massenmatrix eines gekrümmten Timoshenko-Balkens (Viertelkreis) (nu...
Definition: CurvedBeam.m:168
T
Transformationsmatrix für das lokale Koordinatensystem.