Iterationsproblem

[Bjoerk]

Ich hann folschendes Iterationsproblem (Skizze Siehe Anlage). Es soll eine Fläche A gefunden werden, die gewisse Kriterien erfüllt.

Vorgabe ist ein Rechteck b * d sowie eine Belastung Nx, My und Mz.

Das Koordinatensystem x,y,z befindet sich im Schwerpunkt des Querschnitts b * d.

In Richtung der x-Achse drückt jemand (Kraft Nx), an den beiden anderen Achsen drehen welche (Momente My, Mz).

Der Querschnitt soll keine Zugspannungen aufnehmen können, d. h. der Querschnitt reißt auf und es stellt sich eine Fläche A ein, deren Flächeninhalt gleich Nx ist (In Wirklichkeit ein räumlicher Spannungspolyeder, hier aber wurscht).

Damit die Fläche A auch die Momente My und Mz im Gleichgewicht halten kann, muß sie sich wie folgt einstellen (S = Schwerpunkt der Fläche A):

Delphi-Quellcode:

repeat
  A := GetA .. // ??
  S:= GetSchwerpunkt(A);
  C:= GetC; // Konstante
until (A * C = NX) and (A * C * S.y = Mz) and (A * C * S.z = My);

Ich hab das mal vor 20 Jahren mit einer doppelten Repeat Schleife programmiert, kann diese procedure (1000 Zeilen) aber beim besten Willen nicht mehr nachvollziehen.

Anybody?

Crosspost

[Bjoerk]

Falls es jemand konkret durchspielen möchte, folgender Record berechnet die Fläche und den Schwerpunkt:

Delphi-Quellcode:

unit uNEck; // (c) 28.04.1985, 2.05.2012 Thomas Abel, Edgar Zocher

interface

const
  MaxFloatPoints = 100;

type
  TFloatPoint = record
    X, Y: double;
  end;
  TFloatPoints = array [0..MaxFloatPoints - 1] of TFloatPoint;
  TNEck = record
  private
    FCount: integer; // Anzahl der Eckpunkte
    FPoints: TFloatPoints;
    function Determinante(const I: integer): double;
    function GetPoint(Index: integer): TFloatPoint;
    procedure SetPoint(const Index: integer; const X, Y: double);
    function NextPointNumber(const Index: integer): integer;
  public
    procedure Clear;
    procedure AddPoint(const X, Y: double); overload;
    procedure AddPoint(const Index: integer); overload;
    procedure DelPoint(const Index: integer);
    procedure InsPoint(const Index: integer; const X, Y: double);
    function Schwerpunkt: TFloatPoint;
    function Flaeche: double;
    property Points[Index: integer]: TFloatPoint read GetPoint;
    property Count: integer read FCount;
  end;

implementation

function TNEck.GetPoint(Index: integer): TFloatPoint;
begin
  Result := FPoints[Index];
end;

procedure TNEck.SetPoint(const Index: integer; const X, Y: double);
begin
  FPoints[Index].X := X;
  FPoints[Index].Y := Y;
end;

procedure TNEck.AddPoint(const X, Y: double);
begin
  Inc(FCount);
  SetPoint(FCount - 1, X, Y);
end;

procedure TNEck.AddPoint(const Index: integer);
begin
  AddPoint(FPoints[Index].X, FPoints[Index].Y);
end;

procedure TNEck.DelPoint(const Index: integer);
var
  I: integer;
begin
  for I := Index to FCount - 2 do
    FPoints[I] := FPoints[I + 1];
  Dec(FCount);
end;

procedure TNEck.InsPoint(const Index: integer; const X, Y: double);
var
  I: integer;
begin
  Inc(FCount);
  for I := FCount - 1 downto Index + 1 do
    FPoints[I] := FPoints[I - 1];
  SetPoint(Index, X, Y)
end;

function TNEck.NextPointNumber(const Index: integer): integer;
begin
  if Index = FCount - 1 then
    Result := 0
  else
    Result := Index + 1;
end;

function TNEck.Determinante(const I: integer): double;
var
  J: integer;
begin
  J := NextPointNumber(I);
  Result := FPoints[I].X * FPoints[J].Y - FPoints[I].Y * FPoints[J].X;
end;

function TNEck.Flaeche: double;
var
  I: integer;
begin
  Result := 0;
  for I := 0 to FCount - 1 do
    Result := Result + Determinante(I) / 2;
end;

function TNEck.Schwerpunkt: TFloatPoint;
var
  I, J: integer;
  M1, M2: double;
begin
  M1 := 0;
  M2 := 0;
  for I := 0 to FCount - 1 do
  begin
    J := NextPointNumber(I);
    M1 := M1 + (FPoints[I].X + FPoints[J].X) * Determinante(I) / 6;
    M2 := M2 + (FPoints[I].Y + FPoints[J].Y) * Determinante(I) / 6;
  end;
  Result.X := M1 / Flaeche;
  Result.Y := M2 / Flaeche;
end;

procedure TNEck.Clear;
begin
  FCount := 0;
end;

end.

[Jens01]

Booh, Statik.
Du willst das A jetzt bei jeder Iteration verändern, so dass irgendwann alle 3 Bedingungen gültig sind? Dafür mußt Du doch irgendwelche Kriterien haben, damit das Ganze konvergieren kann.
Aus meiner Sicht fehlt noch irgendwas für A.

Gruss Jens

Mich würde ja Dein ursprünglicher Code interessieren

Edit: Habe Deinen 2. Post eben erst gesehen...
Also so soll die Fläche ermittelt werden. Willst Du jetzt irgendwelche Flächen ausprobieren, bis die 3 Bedingungen stimmen? Da fehlt noch was. So ist es ein einfaches Durchprobieren, das fast unendlich wird.

[Jens01]

Habe gerade geguckt.
Wenn es sich um Kernweiten handelt -und ich nicht ganz falsch liege- dann hängt die Form für diesen spannungsfreien Bereich von der Form des Bauteils ab. Bei so einem Rechteck ist da ja so ein Karo. Das Ding könnte man wahrscheinlich verändern und dann nach den Kriterien abprüfen.
Gruss Jens

[Bjoerk]

Hi Jens,

Ja, es gibt dann eine Lösung, sofern S innerhalb der zweiten Kernfläche liegt. Das ist ja der Kippnachweis, der zum Beispiel bei einem Einzelfundament zu führen ist. Bei einer Stahlbetonbemessung macht man auch nichts anderes, nur daß der Stahl hier eventuelle Zugkräfte aufnimmt, es also keine Beschränkung auf irgendwelche Kernweiten gibt.

Den Code kann ich leider nicht veröffentlichen, erstens zu peinlich und zweites Betriebsinterna.

Ich will aber eigentlich nicht fachsimpeln. Mir geht es darum, ob jemand dem Dingen ansieht, "ah, das könnte eine doppelte Newton Iteration werden" oder "mit XY funktionieren".

Gruß
Thomas

[Jens01]

Mir gings nicht ums Fachliche. Bin nur der Meinung, dass noch Kriterien für eine Iteration fehlen, die ich mit meinem vorhergehenden Post erörtern wollte. Ich mag aber auch völlig falsch liegen...

Zitat:

Den Code kann ich leider nicht veröffentlichen, erstens zu peinlich und zweites Betriebsinterna.

Schade, man kann ja mal fragen

[Bjoerk]

Ich weiß jetzt wieder, wie ich’s damals gemacht habe. Oh Mann, eh.

Wie wir schon sagten, es gibt eine Lösung, wenn Nx < 0 und (My/Nx/d)^2 + (Mz/Nx/b)^2 <= 1/9

Am Anfang hat man A = b*d. Dann berechnet man in den Eckpunkten die Spannungen. Wechselt in 2 betrachteten Punkten das Vorzeichen, dann berechnet man den Nullpunkt der Line zwischen diesen beiden Punkten und das ist dann ein korrigierter Eckpunkt der Fläche. Abbruchkriterium ist, wenn sich die Koordinaten der Punkte der Fläche nicht mehr ändern.

Kleiner Ausschnitt aus’m Katastrophencode:

Punkte sind array [1..2, 1..Count]

Delphi-Quellcode:

  J:= 1;
  for K:= 1 to CountAlt do
  begin
    if K = CountAlt then
      I:= 1
    else
      I:= K + 1;
    if Sigma[K] <= 0 then
    begin
      PNeu[1, J]:= P[1, K];
      PNeu[2, J]:= P[2, K];
      Inc(J);
    end;
    if (Sigma[K] * Sigma[I]) < 0 then
    begin
      PNeu[1, J]:= P[1, K] + Sigma[K] * (P[1, I] - P[1, K]) / (Sigma[K] - Sigma[I]);
      PNeu[2, J]:= P[2, K] + Sigma[K] * (P[2, I] - P[2, K]) / (Sigma[K] - Sigma[I]);
      Inc(J);
    end;
  end;
  CountNeu:= J;

Da Spannungen immer bezüglich Xsi/Eta berechnet werden müssen, ist das eine blöde Transponiererei von X/Y nach jeweils Xsi/Eta, das ist alles.

Alles?

Wirklich alles?

Ich fürchte nicht, ich hab' jetzt aber erst mal keinen Bock mehr..

Gruß und schönen Sonntag

Thomas