Основы компьютерного моделирования физических и технических систем. Основы моделирования физических и технических систем, страница 4

//Влияние перемещения DX на координату X нашего узла

CurEqX[DXOffset(T.Points[PointInTria].ID)] :=

CurEqX[DXOffset(T.Points[PointInTria].ID)]+

K[DXOffset(VertexIndex-1), DXOffset(PointInTria-1)];

//Влияние перемещения DY на координату X нашего узла

CurEqX[DYOffset(T.Points[PointInTria].ID)] :=

CurEqX[DYOffset(T.Points[PointInTria].ID)]+

K[DXOffset(VertexIndex-1), DYOffset(PointInTria-1)];

//Влияние перемещения DX на координату Y нашего узла

CurEqY[DXOffset(T.Points[PointInTria].ID)] :=

CurEqY[DXOffset(T.Points[PointInTria].ID)]+

K[DYOffset(VertexIndex-1), DXOffset(PointInTria-1)];

//Влияние перемещения DY на координату Y нашего узла

CurEqY[DYOffset(T.Points[PointInTria].ID)] :=

CurEqY[DYOffset(T.Points[PointInTria].ID)]+

K[DYOffset(VertexIndex-1), DYOffset(PointInTria-1)];

end;

end;

//Добавим сформированные уравнения в систему

SL.AddFullEquation(CurEqX, P.ExtForceX);

SL.AddFullEquation(CurEqY, P.ExtForceY);

end;

//Проверка - уравнений должно было получится столько же, сколько неизвестных

assert(SL.CurEquation = SL.Size+1);

//Решаем СЛАУ одним из методов, записывая результат в Х

SL.Solve(X, Method);

//Сохраним решение в свойствах DX DY узлов объекта 

for VertexID := 0 to Vertexes.Count - 1 do

begin

P := TRealPoint(Vertexes[VertexID]);

P.DX := X[DXOffset(VertexID)];

P.DY := X[DYOffset(VertexID)];

{    P.X := P.X+P.DX;

P.Y := P.Y+P.DY;}{Раскоментить для отображения перемещений}

end;

//Рассчитаем максимальные напряжения в конечных элементах

//По формуле e=B*r, где e - вектор напряжений, B - матрица градиентов,

//r - вектор перемещений (т.е. найденное решение СЛАУ)

Kmax := 0;

SetLength(DX, 6+1);

SetLength(EPS, 3+1);

for TriaID := 0 to Trias.Count - 1 do

begin

T := TTria(Trias[TriaID]);

T.BMatrix(E, v, d, B);

DX[1] := T.Points[1].DX;

DX[2] := T.Points[1].DY;

DX[3] := T.Points[2].DX;

DX[4] := T.Points[2].DY;

DX[5] := T.Points[3].DX;

DX[6] := T.Points[3].DY;

MultMatrixByVector(B, DX, EPS);

T.EpsX := EPS[1];

T.EpsY := EPS[2];

T.GammaX := EPS[3];

if T.EpsSum > Kmax then

Kmax :=T.EpsSum;

end;

end;

{ TRealLine }

function TRealLine.Same(L: TRealLine): boolean;

begin

Result :=

((FirstPoint = L.FirstPoint)and(EndPoint = L.EndPoint))

or

((FirstPoint = L.EndPoint)and(EndPoint = L.FirstPoint))

end;

{ TTria }

function TTria.Area: Extended;

begin

Result :=Abs(0.5*((x(1)-x(3))*(y(2)-y(3)) - (x(2)-x(3))*(y(1)-y(3))));

end;

procedure TTria.BMatrix(E, v, d: Extended; var B: TMatrix);

begin

SetSize(3, 6, B);

B[1, 1] := y(2)-y(3);

B[1, 3] := y(3)-y(1);

B[1, 5] := y(1)-y(2);

B[2, 2] := x(3)-x(2);

B[2, 4] := x(1)-x(3);

B[2, 6] := x(2)-x(1);

B[3, 1] := x(3)-x(2);

B[3, 2] := y(2)-y(3);

B[3, 3] := x(1)-x(3);

B[3, 4] := y(3)-y(1);

B[3, 5] := x(2)-x(1);

B[3, 6] := y(1)-y(2);

MultByReal(1/2/Area, B);

end;

function TTria.EpsSum: Extended;

begin

Result := Sqr(EpsX) + Sqr(EpsY) + Sqr(GammaX/3);

end;

function TTria.FindPoint(P: TRealPoint): Integer;

var

I: Integer;

begin

for I := 1 to 3 do

if Points[I] = P then

begin

Result := I;

exit

end;

Result := -1;

end;

//ГЛАВНАЯ ПРОЦЕДУРА ФОРМИРОВАНИЯ МАТРИЦЫ ЖЕСТКОСТИ

procedure TTria.KMatrix(E, v, d: Extended; var K: TMatrix);

var

Alpha: Extended;

B, C, BT, R: TMatrix;

//  I, J: Integer;

begin

Alpha := E*d/4/Area/(1-v*v);

SetSize(3, 6, B);

B[1, 1] := y(2)-y(3);

B[1, 3] := y(3)-y(1);

B[1, 5] := y(1)-y(2);

B[2, 2] := x(3)-x(2);

B[2, 4] := x(1)-x(3);

B[2, 6] := x(2)-x(1);

B[3, 1] := x(3)-x(2);

B[3, 2] := y(2)-y(3);

B[3, 3] := x(1)-x(3);

B[3, 4] := y(3)-y(1);

B[3, 5] := x(2)-x(1);

B[3, 6] := y(1)-y(2);

SetSize(3, 3, C);

C[1, 1] := 1;

C[2, 2] := 1;

C[1, 2] := v;

C[2, 1] := v;

C[3, 3] := (1-v)/2;

Transpose(B, BT);

MultMatrix(BT, C, R);

MultMatrix(R, B, K);

MultByReal(Alpha, K);

{  for I := 1 to 6 do

for J := 1 to 6 do

K[I, J] := 1;

}

end;

function TTria.Same(T: TTria): boolean;

var

I, J: Integer;

Match: set of 1..3;

begin

Match := [];

for I := 1 to 3 do

begin

for J := 1 to 3 do

if (not (J in Match)) and (T.Points[I] = Points[J]) then

begin

Include(Match, J);

break;

end;

end;

Result := Match = [1..3];

end;

function TTria.X(I: Integer): Extended;