12PosobKR2015 (В.А. Столярчук. Решение задач в Sigma (учебное пособие)), страница 10

Но прежде всего надо «отловить» номера узлов, лежащих на нагруженной стороне. Отыскиваем нужные номера узлов по координатам.

В массиве CORD узлы располагаются в порядке возрастания номеров.

Последовательно пробегая массив CORD:

DO 10 I=1,NP

(где NP – число узлов) мы, понятно, сначала найдем узел (принадлежащий линии нагружения) с меньшим номером и поместим номер этого узла в массив отобранных номеров IFORCE1(I1)=I, потом с большим номером и его тоже поместим в массив номеров отобранных узлов и т.д. Следовательно, в массиве IFORCE1(I1)=I номера отобранных узлов будут возрастать:

2, 5, 7, 8, 9, 10, 16

Но относительно функции нагрузки они должны стоять в обратном порядке, а именно: 16, 10, 9, 8, 7, 5, 2

Понятно, что для такой операции сложная сортировка не нужна: задача решается очень просто, так как последний ряд чисел является записанным первым рядом, но в обратном порядке.

Произошло это потому, что в данном случае получился простейший вариант нумерации узлов, так как в исходной геометрической модели зоны, образующие нагруженную линию, выстроились слева направо по порядку 7,6,5, из-за того, что автор зону 5 образовал со стороны зоны 3, затем со стороны зоны 5 создал зону 6, после чего со стороны зоны 6 образовал зону 7.

Очевидно, что порядок образования зон мог быть иным. Пример другого образования геометрической модели и результата расчета на ней приведён ниже (рис. 44).

Как видим, результаты расчета не изменились, а вот нумерация узлов в интересующей нас области наоборот: узлы на нагруженной верхней стороне выстроились в ином порядке, а именно: 2, 5, 9, 8, 7, 10, 16. И здесь уже ни о каком возрастании в порядке увеличения номеров не идёт речь.

Очевидно, что после отбора узлов, принадлежащих верхней горизонтальной стороне, в массиве IFORCE1 они выстроятся в порядке возрастания номеров: 2, 5, 7, 8, 9, 10, 16, хотя для решения задачи их необходимо выстроить в порядке 16, 10, 7, 8, 9, 5, 2.

И здесь простейшей перезаписью ряда чисел не обойдёшься.

Именно поэтому здесь необходима сортировка массива IFORCE1 в порядке убывания координаты Х.

Итак, алгоритм подсчета и распределения сил по узлам в случае распределенной нагрузки распадается на три части. В первой отыскиваются по координатам номера узлов к которым приложены силы согласно действующей на этом участке функции нагрузки и заносятся в отдельный массив. Во второй происходит сортировка массива с номерами нагруженных узлов. В третьей находятся, собственно, силы, действующие в узлах.

Есть ещё один аспект, который необходимо цчитывать при составлении программы: проверка правильности вычислений. Для этого приходится сначала вычислить суммарную силу от каждой нагрузки по участкам с помощью какой-либо посторонней программы и сравнить с теми значениями, которыне дает подпрограмма Force при разных NRC. Дело в том, что программа должна подсчитывать нагрузку на каждом из участков с одной и той же степенью точности при разных NRC. Иначе получится, что при одних NRC рассчитывается одна задача, при других – другая. Следовательно в программе должно быть предусмотрен вывод информации отдельно для каждого из участков. Для упрощения программы лучше всего проводить расчет для каждого из участков отдельно при разных NRC, отключая нагрузку по другим участкам. Для отключения/включения нагрузки по участкам в Sigma используется аппарат свободных параметров, с помощью которых осуществляются эти функции. Значения парамеров PRM1, PRM2, PRM3, PRM4, PRM5, PRM6 задаются в окне основных параметров программы (при этом можно использовать только их часть и можно, вообще, обойтись одним параметром, задавая ему разные значения) (рис. 45).

В рассматриваемом примере автор решил задавать управляющие свободные параметры в следующем виде:

Свободные параметры.

Понятно, что значения этих параметров должны передаваться в подпрограмму Force.

Обратите внимание, что в окне рис. 45 задаются и характеристики материала из которого изготовлена пластина. Харатеристики материала вместе с толщиной называют свойством конечного элемента. Если характеристики материала или толщина меняются в разных подобластях рассчитываемого объекта, то задается другое свойство КЭ. Программа расчитана на использование трёх свойств.

Текст продпрограммы Force для рассматриваемого примера приведен ниже.

* Содержит подпрограмму FORCE

* ПOДПPOГPAMMA PAЗHECEHИЯ HAГPУЗOK ПО УЗЛАМ

* BЫЗЫBAETCЯ ИЗ MAIN, вызываемых модулей нет.

====================================================

*PLOSH-подпрограмма вычисления интеграла методом трапеций

*HEIGH-функция вычисления интеграла методом трапеций

*IPR массив управления печатью

*IFORCE_ - массив состоящий из номеров точек и координат точек к которым приложена распределенная сила

*CORD - одномерный массив глобальных координат узлов (первые два элемента этого массива являются координатами узла 1 по оси Х и Y соответственно, вторая пара элементов - узла 2 и т.д.).

*R - массив значений усилий в узлах (аналогичен по структуре массиву CORD, только вместо координат в нем указываются значения нагрузок по оси Х и Y в узлах);

===================================================

* начало кода FORCE

===================================================

SUBROUTINE FORCE (NB, CORD, NRC, DB, IPR, NP, RSUM, R,

> NDF, DH, NR, PRM1, PRM2, PRM3, PRM4, PRM5, PRM6)

LOGICAL*1 IPR(50)

DIMENSION CORD(1), R(1), IFORCE2(100), IFORCE1(100), IFORCEO(100)

* функция вычисления интеграла методом трапеций для прямой Y=60.

* входные значения: (шаг(SH), начальная точка отрезка(XP1), конечная точка отрезка(XH))

RFORSES(SH,XP1,XH,PAR)=(PAR*SH*((XP1**2.+600.)+(XH**2.+600.))/2.)

* NR - количество нагруженных узлов. * NP - количество узлов

I1=0

I2=0

DO 10 I=1,NP

XT=CORD(NDF*(I-1)+1)

YT=CORD(NDF*(I-1)+2)

* Обнуление массива сил

R(NDF*(I-1)+1)=0

R(NDF*(I-1)+2)=0

* Узлы на прямой Y=60

IF ((ABS(YT-60).LT.0.1).AND.(ABS(PRM1-1). LT.0.1)) THEN

I1=I1+1

* Массив IFORCE1 содержит номера нагруженных узлов на прямой Y=60

IFORCE1(I1)=I

NR=NR+1

ENDIF

* Узлы на прямой Y+3X-180=0, psi[0;31.42] x[40; 50.06]

IF (ABS(PRM2-1).LT.0.1) THEN

IF ((ABS(YT+3*XT-180).LT.0.1).AND.(XT-50.06-0.1).LT.0) THEN

I2=I2+1

*Массив IFORCE2 содержит нагруженные узлы на прямой 3X+Y-180=0

IFORCE2(I2)=I

NR=NR+1

ENDIF

ENDIF

* Узлы на прямой Y+3X-180=0, psi[31.42;63.25] x[50.06; 60]

IF (ABS(PRM2-2).LT.0.1) THEN

IF ((ABS(YT+3*XT-180).LT.0.1).AND.((XT-50.06+0.1).GT.0)) THEN

I2=I2+1

IFORCE2(I2)=I

NR=NR+1

ENDIF

ENDIF

* Узлы на прямой Y+3X-180=0, psi[0;63.25] x[40;60]

IF (ABS (PRM2-3). LT. 0.1) THEN

IF (ABS (YT+3*XT-180). LT.0.1) THEN

I2=I2+1

IFORCE2(I2)=I

NR=NR+1

ENDIF

ENDIF

10 CONTINUE

* Сортировка массивов. Во всех трех случаях используется сортировка пузырьком

* Массив IFORCE1 содержит номера нагруженных узлов на прямой Y=60

DO 201 I=I1,2,-1

DO 201 J=1,I-1

IF (CORD(NDF*(IFORCE1(J+1)-1)+1).LT.

CORD(NDF*(IFORCE1(J)-1)+1)) THEN

TMPI=IFORCE1(J+1)

IFORCE1(J+1)=IFORCE1(J)

IFORCE1(J)=TMPI

ENDIF

201 CONTINUE

*Массив IFORCE2 содержит нагруженные узлы на прямой 3X+Y-180=0

DO 202 I=I2,2,-1

DO 202 J=1,I-1

IF (CORD(NDF*(IFORCE2(J+1)-1)+1).GT.

>CORD(NDF*(IFORCE2(J)-1)+1)) THEN

TMPI=IFORCE2(J+1)

IFORCE2(J+1)=IFORCE2(J)

IFORCE2(J)=TMPI

ENDIF

202 CONTINUE

* Задание распределенной силы для прямой Y=60

DO 300 I=1,I1

J=IFORCE1(I)

* Координата текущей точки

X1=CORD(2*(J-1)+1)

Y1=CORD(2*(J-1)+2)

* Координата предыдущей точки

IF (I.NE.1) THEN

X3=CORD(2*(IFORCE1(I-1)-1)+1)

Y3=CORD(2*(IFORCE1(I-1)-1)+2)

ELSE

X3=X1

Y3=Y1

ENDIF

* Координата следующей точки

IF (I.NE.I1) THEN

X2=CORD(2*(IFORCE1(I+1)-1)+1)

Y2=CORD(2*(IFORCE1(I+1)-1)+2)

ELSE

X2=X1

Y2=Y1

ENDIF

* Сила, приложенная к текущему узлу (вычисляется через площадь трапеции)

XP1=ABS((X1+X2)/2.-40.)

XP2=ABS((X1+X3)/2.-40.)

SH=ABS(XP1-XP2)/10.

XH = XP1+SH

R1=0.

DO 500 IS=1,10,1

R1=R1+ABS(RFORSES(SH,XP1,XH,PRM4))

XP1=XP1+SH

! PRINT*, R1, XP1, XP2

XH=XH+SH

500 CONTINUE

! поскольку сила действует лишь по оси Y, то нагрузка полностью идет на ось Y

R(NDF*(IFORCE1(I)-1)+2)=R(NDF*(IFORCE1(I)-1)+2)+R1

! блокируем вывод точки, на которую действует нагрузка с прямой Y=180-3X

IF(INT(PRM2).EQ.2) THEN

IF(IPR(23)) WRITE(6,21) J,CORD(NDF*(J-1)+1),

>CORD(NDF*(J-1)+2),R((J-1)*NDF+1),R((J-1)*NDF+2)

ELSE IF(X1.LT.39.9) THEN

IF(IPR(23)) WRITE(6,21) J,CORD(NDF*(J-1)+1),

>CORD(NDF*(J-1)+2),R((J-1)*NDF+1),R((J-1)*NDF+2)

ELSE

NR=NR-1

ENDIF

300 CONTINUE

* Задание распределенной силы для прямой Y+3X-180=0

DO 400 I=1,I2

J=IFORCE2(I)

J_PRED=IFORCE2(I+1)-1

* Координата текущей точки

X1=CORD(2*(IFORCE2(I)-1)+1)

Y1=CORD(2*(IFORCE2(I)-1)+2)

* Координата предыдущей точки

IF (I.NE.1) THEN

X3=CORD(2*(IFORCE2(I-1)-1)+1)

Y3=CORD(2*(IFORCE2(I-1)-1)+2)

ELSE

X3=X1

Y3=Y1

ENDIF

* Координата следующей точки

IF (I.NE.I2) THEN

X2=CORD(2*(IFORCE2(I+1)-1)+1)

Y2=CORD(2*(IFORCE2(I+1)-1)+2)

ELSE

X2=X1

Y2=Y1

ENDIF

* Сила приложенная к текущему узлу (вычисляется через площадь трапеции)

* R1-распределенная сила приложенная к текущему узлу

R1=RPLOSH(X3,X1,X2, PRM3) !(предыдущая, текущая, следующая)

! PRINT*, R1, R2

IF (ABS(CORD(2*(IFORCE2(I)-1)+1)-60.).GT.0.1) THEN

R(NDF*(IFORCE2(I)-1)+1)=R(NDF*(IFORCE2(I)-1)+1)

>+R1*SQRT(1-1/(1+(3.)**2)) ! коэф. 3 взят из тангенса(коэф. - k) уравнения прямой