» 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:

iIDLayerNodeSiIDLayerNodeSiIDLayerNodeS
01BottomN1X241TopN1X481MiddleN1X
11BottomN1Y251TopN1Y491MiddleN1Y
21BottomN1Z261TopN1Z501MiddleN1Z
31BottomN1XY271TopN1XY511MiddleN1XY
41BottomN1YZ281TopN1YZ521MiddleN1YZ
51BottomN1XZ291TopN1XZ531MiddleN1XZ
61BottomN2X301TopN2X541MiddleN2X
71BottomN2Y311TopN2Y551MiddleN2Y
81BottomN2Z321TopN2Z561MiddleN2Z
91BottomN2XY331TopN2XY571MiddleN2XY
101BottomN2YZ341TopN2YZ581MiddleN2YZ
111BottomN2XZ351TopN2XZ591MiddleN2XZ
121BottomN3X361TopN3X601MiddleN3X
131BottomN3Y371TopN3Y611MiddleN3Y
141BottomN3Z381TopN3Z621MiddleN3Z
151BottomN3XY391TopN3XY631MiddleN3XY
161BottomN3YZ401TopN3YZ641MiddleN3YZ
171BottomN3XZ411TopN3XZ651MiddleN3XZ
181BottomN4X421TopN4X661MiddleN4X
191BottomN4Y431TopN4Y671MiddleN4Y
201BottomN4Z441TopN4Z681MiddleN4Z
211BottomN4XY451TopN4XY691MiddleN4XY
221BottomN4YZ461TopN4YZ701MiddleN4YZ
231BottomN4XZ471TopN4XZ711MiddleN4XZ
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).

KomponentResultReaderResult
X-47.4161763191223-47.416
Y7.89010169316312E-207.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.