» Blog » Mechanika » Odczyt danych z pliku RST i transformacja tensora naprężeń

Odczyt danych z pliku RST i transformacja tensora naprężeń

Pisanie skryptów w środowisku Mechanical nie musi ograniczać się do prostej automatyzacji czynności niezbędnych do wykonania analizy MES. Dzięki narzędziom takim jak np. ResultReader możemy uzyskać dostęp do surowych wyników w pliku RST i przetworzyć wyniki w sposób specyficzny dla naszych potrzeb, co zwykle wymaga implementacji szeregu funkcji matematycznych. Niniejszy artykuł przybliży dwa zaimplementowane w Mechanical API obiekty, których wykorzystanie pozwoli zaoszczędzić nieco czasu: Vector3D i Matrix4D.

Z uwagi na zakres zagadnienia, artykuł został podzielony na 3 części, z których każda porusza kolejny etap pracy:

Część I – Obiekty Vector3D i Matrix4D: zawiera krótkie wprowadzenie do programowania obiektowego – jak je tworzyć i jak z nich korzystać.

Część II – Transformacja koordynat do lokalnego układu współrzędnych: pełen opis przykładu kodu transformującego koordynaty punktu 3D do lokalnego układu współrzędnych z wykorzystaniem opisanych wcześniej obiektów Vector3D i Matrix4D.

Część III – Odczyt danych z pliku RST i transformacja tensora naprężeń

Odczyt danych z pliku RST i transformacja tensora naprężeń

Opis procesu odczytu danych z pliku RST (na bazie elementów powłokowych) oraz ich transformacji do lokalnego układu współrzędnych z wykorzystaniem obiektów Vector3D i Matrix4D.

  1. Odczyt danych z pliku RST
  2. Transformacja tensora naprężeń
  3. Finalny kod

Odczyt danych z pliku RST

Najprostszym sposobem na uzyskanie danych numerycznych, jest utworzenie obiektu typu Result, a następnie dostęp do danych z wykorzystaniem własności PlotData. Metoda ta wymaga uprzedniego utworzenia i konfiguracji odpowiednich obiektów. Proces ten jest intuicyjny, ale w niektórych sytuacjach może wymagać utworzenia dużej ilości obiektów.

W sytuacjach, kiedy zależy nam na dobrej optymalizacji możemy pominąć tworzenie obiektów typu Result i odczytać dane bezpośrednio z pliku RST. Metoda ta jest szybsza, ale bardziej wymagająca od strony przetwarzania danych. W tym akapicie postaram się przybliżyć kwestię transformacji tensora naprężeń do lokalnego układu współrzędnych. W artykule będę posiłkował się prostą analizą wykonaną w ANSYS Mechanical 2023R2.

Transformacja tensora naprężęń: model testowy
Rys. 1 Model testowy

Zacznijmy od odczytu danych z pliku RST – API Mechanical udostępnia klasę ResultReader dostępną z poziomu dowolnej analizy. Aby uniknąć problemów z prawami dostępu, klasę najlepiej inicjalizować w bloku with:

meshData = DataModel.MeshDataByName('Global')
eList = [57]
eInfo = [meshData.ElementById(e) for e in eList]
with ExtAPI.DataModel.AnalysisList[0].GetResultsData() as r:

stress = r.GetResult('S')

sComponents = ['X','Y','Z','XY','YZ','XZ']

noComponents = len(sComponents)

stress.SelectComponents(sComponents)

Powyższy kod tworzy obiekt umożliwiający dostęp do naprężeń (S). Metoda SelectComponents umożliwia wybranie interesujących nas komponentów (w tym przypadku wszystkich). Istotna jest kolejność symboli komponentów w liście. Na tym etapie natrafiamy na pierwszy element, który może sprawić trudności mniej zaawansowanym użytkownikom.

vs = stress.GetElementValues(eList, False)

Powyższa metoda pobiera komponenty naprężeń wg. listy sComponents dla elementów, których ID znajdują się w zmiennej eList. Wynikiem jest lista, która zawiera wartości z podziałem na elementy, warstwy, węzły oraz komponenty (w tej kolejności). Aby lepiej zobrazować kolejność wartości posłużę się przykładem. Lista wynikowa dla jednego elementu powłokowego definiowanego przez 4 węzły zawiera 72 wartości: 1 element x 3 warstwy x 4 węzły x 6 komponentów = 72.

Transformacja tensora naprężeń: Lokalizacja elementu testowego
Rys. 2 Lokalizacja elementu testowego

Dla przejrzystości, zamieszczamy pełną strukturę listy:

i
ID
Layer
Node
S
i
ID
Layer
Node
S
i
ID
Layer
Node
S
0
1
Bottom
N1
X
24
1
Top
N1
X
48
1
Middle
N1
X
1
1
Bottom
N1
Y
25
1
Top
N1
Y
49
1
Middle
N1
Y
2
1
Bottom
N1
Z
26
1
Top
N1
Z
50
1
Middle
N1
Z
3
1
Bottom
N1
XY
27
1
Top
N1
XY
51
1
Middle
N1
XY
4
1
Bottom
N1
YZ
28
1
Top
N1
YZ
52
1
Middle
N1
YZ
5
1
Bottom
N1
XZ
29
1
Top
N1
XZ
53
1
Middle
N1
XZ
6
1
Bottom
N2
X
30
1
Top
N2
X
54
1
Middle
N2
X
7
1
Bottom
N2
Y
31
1
Top
N2
Y
55
1
Middle
N2
Y
8
1
Bottom
N2
Z
32
1
Top
N2
Z
56
1
Middle
N2
Z
9
1
Bottom
N2
XY
33
1
Top
N2
XY
57
1
Middle
N2
XY
10
1
Bottom
N2
YZ
34
1
Top
N2
YZ
58
1
Middle
N2
YZ
11
1
Bottom
N2
XZ
35
1
Top
N2
XZ
59
1
Middle
N2
XZ
12
1
Bottom
N3
X
36
1
Top
N3
X
60
1
Middle
N3
X
13
1
Bottom
N3
Y
37
1
Top
N3
Y
61
1
Middle
N3
Y
14
1
Bottom
N3
Z
38
1
Top
N3
Z
62
1
Middle
N3
Z
15
1
Bottom
N3
XY
39
1
Top
N3
XY
63
1
Middle
N3
XY
16
1
Bottom
N3
YZ
40
1
Top
N3
YZ
64
1
Middle
N3
YZ
17
1
Bottom
N3
XZ
41
1
Top
N3
XZ
65
1
Middle
N3
XZ
18
1
Bottom
N4
X
42
1
Top
N4
X
66
1
Middle
N4
X
19
1
Bottom
N4
Y
43
1
Top
N4
Y
67
1
Middle
N4
Y
20
1
Bottom
N4
Z
44
1
Top
N4
Z
68
1
Middle
N4
Z
21
1
Bottom
N4
XY
45
1
Top
N4
XY
69
1
Middle
N4
XY
22
1
Bottom
N4
YZ
46
1
Top
N4
YZ
70
1
Middle
N4
YZ
23
1
Bottom
N4
XZ
47
1
Top
N4
XZ
71
1
Middle
N4
XZ
Tab. 1 Struktura danych uzyskana funkcją GetElementValues(..) dla elementu powłokowego

Struktura powtarza się dla każdego kolejnego elementu. Należy zaznaczyć, że struktura zmienia się w zależności od ilości węzłów oraz typu elementu. Dodatkowo, dla elementów belkowych znaczenie ma ustawienie opcji Beam Section Results dla obiektu Solution. Dla dociekliwych, pełną dokumentację można znaleźć w pomocy ANSYS.

*Dla osób bez konta ANSYS (Customization Suite 2023 R2 -> ACT Customization Guide for Mechanical -> MECHANICAL FEATURE CREATION -> POSTPROCESSING CAPABILITIES IN MECHANICAL -> Retrieving Mechanical Results More Efficiently).

Wracając do naszego skryptu, dla uproszczenia będziemy operować na wartościach uśrednionych dla elementu (Elemental Mean) w środkowej warstwie pojedynczego elementu. Ułatwi to weryfikację oraz wyeliminuje z kodu niepotrzebne pętle.

Na tym etapie mamy zmienną v, która zawiera listę 72 wartości dla elementu o ID=57.

nCnt = len(e.CornerNodeIds)
    iC = 0
    for component in sComponents:
        result[component] = sum(v[iC::noComponents][nCnt*2:nCnt*3])/nCnt
        iC += 1

Powyższy skrypt iteruje po komponentach po czym dwukrotnie zawęża listę do wartości dla danych komponentów na wymaganej warstwie. Na początku, zawężamy listę do wartości dla danego komponentu ze wszystkich warstw. Rzut oka na tabelkę i widzimy, że co 6 element w liście reprezentuje wartość danego komponentu. Wyrażenie:

v[0::6]

pozwala na utworzenie nowej listy każdego 6 elementu zaczynając od indeksu 0. Poniżej prosty przykład obrazujący wybieranie n-tego elementu z listy:

l = range(0,11)
# [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
l[0::3]
# [0, 3, 6, 9]
l[1::3]
# [1, 4, 7, 10]

Odpowiednio – v[1::6] zwróci co 6 wartość poczynając od drugiego komponentu (Y).

Ostatecznie, dla komponentu X otrzymamy listę 12 wartości (dla 4 węzłów oraz 3 warstw). Interesująca nas warstwa jest ostatnia. Musimy również pamiętać, że długość listy zależy od liczby węzłów danego elementu. Jeśli zmienna nCnt zawiera informację o liczbie węzłów to:

#l = ['N1_B', 'N2_B', 'N3_B', 'N4_B', 'N1_T', 'N2_T', 'N3_T', 'N4_T', 'N1_M', 'N2_M', 'N3_M', 'N4_M']
l[nCnt*2:nCnt*3]
#['N1_M', 'N2_M', 'N3_M', 'N4_M']

Na koniec sumujemy wartości w liście metodą sum oraz dzielimy przez liczbę węzłów. Tak uśrednioną wartość wstawiamy do słownika używając symbolu komponentu jako klucza.

Poniższa tabelka prezentuje porównanie wartości pobranych bezpośrednio z pliku RST z wartościami uzyskanymi z wykorzystaniem obiektów Result (dla elementu o ID=57 załączonego przykładu).

Komponent
ResultReader
Result
X
-47.4161763191223
-47.416
Y
7.89010169316312E-20
7.8901e-020
Z
-1.65947604179382
-1.6595
XY
-0.000511929980348214
-5.1193e-004
YZ
-0.000644275627564639
-6.4428e-004
XZ
-31.5254945755005
-31.525
Tab. 2 Porównanie wartość uzyskanych obiektem Result oraz skryptem

Transformacja tensora naprężeń

Po odczytaniu wartości naprężeń dla elementu możemy utworzyć obiekt Matrix4D reprezentujący tensor naprężeń w postaci:

 

Poniższy kod obrazuje wykorzystanie utworzonego wcześniej słownika oraz konstruktora kolumnowego:

sigma = Matrix4D.CreateSystem(

Vector3D(result['X'], result['XY'], result['XZ']),

Vector3D(result['XY'], result['Y'], result['YZ']),

Vector3D(result['XZ'], result['YZ'], result['Z']),

)

Aby uzyskać wartości naprężeń w lokalnym układzie współrzędnych, wykorzystamy macierz obrotu Q. Operacje mnożenia macierzy, które należy wykonać możemy zapisać w następujący sposób:

 

Biorąc pod uwagę, że tworząc macierz obrotu za pomocą metody CreateSystem tworzymy macierz już transponowaną, musimy ją skopiować i transponować, aby otrzymać oryginalną macierz obrotu.

local_cs = DataModel.GetObjectsByName('newCS_1')[0]
Q_T = Matrix4D.CreateSystem(

Vector3D(*local_cs.XAxis),

Vector3D(*local_cs.YAxis),

Vector3D(*local_cs.ZAxis))

Q = Q_T*Matrix4D.Identity
Q.Transpose()
sigma2 = Q*sigma*Q_T
print(sigma2)
#|-0.402686442266813, 0.0128198667790177, 0.515360592656693, 0 |
#|0.0128198667790177, 0.179480241292891, 0.308007133267354, 0 |
#|0.515360592656693, 0.308007133267354, -0.353742978198594, 0 |
#| 0, 0, 0, 1 |

Finalny kod

Ostateczny kod możemy podzielić na 2 części, część odczytu danych z pliku RST:

meshData = DataModel.MeshDataByName('Global')
eList = [57]
eInfo = [meshData.ElementById(e) for e in eList]
 
noLayers = 3
sComponents = ['X','Y','Z','XY','YZ','XZ']
noComponents = len(sComponents)
noElements = len(eList)
 
result = {}
 
with ExtAPI.DataModel.AnalysisList[0].GetResultsData() as r:

stress = r.GetResult('S')

stress.SelectComponents(sComponents)

vs = stress.GetElementValues(eList, False)

e = eInfo[0]

nCnt = len(e.CornerNodeIds)

iC = 0

for component in sComponents:

result[component] = sum(vs[iC::noComponents][nCnt*2:nCnt*3])/nCnt

iC += 1

Oraz część transformująca naprężenia w zmiennej results:

sigma = Matrix4D.CreateSystem(

Vector3D(result['X'], result['XY'], result['XZ']),

Vector3D(result['XY'], result['Y'], result['YZ']),

Vector3D(result['XZ'], result['YZ'], result['Z']),

)
local_cs = DataModel.GetObjectsByName('newCS_1')[0]
Q_T = Matrix4D.CreateSystem(

Vector3D(*local_cs.XAxis),

Vector3D(*local_cs.YAxis),

Vector3D(*local_cs.ZAxis))

Q = Q_T*Matrix4D.Identity
Q.Transpose()
sigma2 = Q*sigma*Q_T

Autor: Bartosz Płochocki, MESco sp. z o.o.

Zainteresował Cię temat? Zobacz nasze szkolenia poświęcone skryptowaniu na poziomie podstawowym i zaawansowanym.