Kielce University of Technology
1 e-mail: ksuckert@tu.kielce.pl 2 e-mail: zmbul@tu.kielce.pl
VALIDATION OF THE ADEQUACY OF FEM
REPRESENTATION OF RECTANGULAR MEMBRANE
ELEMENT WITH AN ADDITIONAL ROTATIONAL
DEGREE OF FREEDOM
A b s t r a c t
The paper provides an analysis of a membrane element with an additional rotational degree of freedom. Similar elements, the stiffness matrices obtained in an exact procedure, were investigated in other studies. The paper, however, analyses the membrane element stiffness matrix obtained in an approximate procedure with the use of reduced integration. Such an approach is an attempt at eliminating shear locking that occurs when a membrane element with additional rotational degrees of freedom is applied.
Keywords: finite element, membrane, rotational degree of freedom, stiffness matrix.
1. Introduction
With the widespread use of the finite element method in engineering calculations, fast advances in computer methods, and design of increasingly complex structures, it is necessary to develop new finite elements and to minimise computation errors. The aim can be achieved in two ways, i.e. by using a finer mesh and by introducing additional degrees of freedom in the element.
In some structures, none of the standard elements sufficiently reflects the character of the structure performance. Attempts are therefore made to create increasingly perfect elements, i.e. non-standard ones.
Fig. 1. Membrane element with an additional rotational
Such elements are employed when elements of different type, e.g. membrane and bar ones are jointed together. Then, additional degrees of freedom occur in a common node (nodes). The paper analyses a membrane element with an additional rotational degree of freedom, shown in Figure. 1. Also, the adequacy of the element representation is assessed. 2. Problem formulation
The stiffness matrix of the membrane element with an additional rotational degree of freedom is subjected to analysis to assess the adequacy of the finite element. For the element presented above, the stiffness matrix derived in an exact procedure can be found in studies [2-4].
Present considerations involve a matrix derived in a reduced procedure. The concept of reduced integration is applied to final element analysis because such an approach sometimes produces better results when compared with the exact solution.
The stiffness matrix analysed in the paper was obtained using four-point Gauss quadrature. The matrix derived in a reduced procedure was only slightly different from the exact one. Exemplary expressions k11 obtained with both methods are juxtaposed in the Table 1.
Table 1. Elements of the matrices obtained in an exact and reduced procedures
Element of the matrix obtained
in an exact procedure in an approximate (reduced) procedureElement of the matrix obtained
11 2 2
0.48571 0.2 0.2
k = a + −b ν 11 2 2
0.485625 0.2 0.2
k = a + b− ν
The matrix derived while using reduced integration was validated by verifying static equilibrium equations and employing the criteria stated by Gilewski [1], which should be satisfied by a correctly formulated finite element. The criteria include an ellipticity condition and a congruence condition (energy criterion). The first condition is obligatory, so elements that do not satisfy it should not be used at all, but an element does not need to be eliminated because of not fulfilling the congruence condition. 3. Validation of the adequacy of the formulation
of the stiffness matrix of a membrane element with an additional rotational degree of freedom 3.1. Verifying static equilibrium equations
The verification is only basic in character and it is the simplest one, yet it allows a quick validation of the correctness of the stiffness matrix. It should be noted that not fulfilling this condition indicates an incorrect formulation of the matrix, meeting the condition, however, does not ensure correctness, so further verification is required. Thus, such verification provides only rough assessment.
In the stiffness matrix, each column describes nodal forces in displacement unit state. The equilibrium equations are verified by applying those forces to appropriate nodes. In accordance with this criterion, the correctness of the matrix is indicated by zeroing of the sum of the projections of the forces on x and y axes and also of the moment with respect to an arbitrary point. As regards the matrix under consideration, the criterion mentioned above was fulfilled for each column.
3.2. Verifying the ellipticity condition
To verify the ellipticity condition, it is necessary to find eigenvalues and eigenvectors of the stiffness matrix. In accordance with [1], the ellipticity condition is equivalent to satisfying the following points:
1. Eigenvalues of the stiffness matrix must be non-negative
2. The spectrum of eigenvalues of the stiffness matrix must contain exactly the same number of zero eigenvalues as the number of rigid motions of elements.
3. Eigenvectors describing rigid motion must correspond to zero eigenvalues.
4. Eigenvectors that describe deformed shapes of the element must correspond to non-zero eigenvalues.
To verify the dependences presented above, it is necessary to solve an eigenproblem, which is formulated as follows:
(
K−λI q)
=0 (1) Using the Mathematica software and assuming the element dimensions: a = 0.4m, b = 0.2m, t = 0.01m, and materials characteristics E = 192.6GPa and υ = 0.3, eigenvalues and eigenvectors for the matrix derived in a reduced procedure were computed. For the matrix under consideration, all the eigenvalues are specified positive or are zero. In addition, three zero values are found, which corresponds to the number of rigid motions for 2D elements. Thus, points 1 and 2 given above are met, which is necessary to satisfy the ellipticity criterion. The verification of the two remaining points will be possible owing to the Figure 2, which shows eigenvectors for individual eigenvalues. The obtained eigenvalues are as follows:λ={26.9158, 11.7113, 9.69056, 6.15919, 4.89254, 0.157527, 0, 0, 0}
It can be observed that the vectors corresponding to non-zero eigenvalues (λ1 – λ6) describe deformed shapes of the element, whereas the vectors corresponding to zero eigenvalues (λ7 – λ9) describe rigid motions. Thus, points 3 and 4 are met, which is necessary to satisfy the ellipticity criterion.
In addition, points 3 and 4 can be verified analytically by computing the strains corresponding to eigenvectors – εi.
i i=Bf
ε (2)
where B is the strain matrix, and fi – an eigenvector corresponding to a given eigenvalue. Strain vectors corresponding to rigid motions are zero, whereas the vectors describing the deformed shapes are described by higher-order polynomials. For instance, vectors ε1 andε8have the form:
2 3 2 1 1 3 2 2 2 2 0.503614 0.100623 0.28696 0.095633 5.68815 0.117333 2.63317 ( ) ( 0.877725 , 1.08917 0.57392 2.83442 3.78101 0.57392 1.31659 1.31659 ) η η η ξ ξ ε λ ξ η η ξ ηξ η ξ ξ ηξ − − + − − + + = − − + − + + − + − (3) = 0 0 0 ) ( 8 8 λ ε
Like ε1 given above, any other vector corresponding to the element deformed shape, has a term, which is a constant, and which provides permanency of deformations.
The analytical calculations that were performed make it possible to draw conclusions identical to those for the Figure 2. The analysis of eigenvalues and eigenvectors of the stiffness and strain matrices provides an unambiguous proof that all four points corresponding to the ellipticity criterion were satisfied.
3.3. Verifying the congruence condition
This verification, also termed the energy criterion, involves the comparison of the energy accumulated
in an element being considered computed from the FEM representation with the exact solution:
{
− + + + − = 2 2 2 2 0 2 1 2 1 2 1 xy y x y x Et U ε ε νε ε ν γ ν (4)Where: (4), εx, εy – denote linear strains, γxy – shape
deformation.
To verify the correctness of a given element, the elastic strain energy was computed with respect to its gravity centre, in accordance with the energy criterion.
The computations performed with the Mathematica software for the matrix obtained by reduced integration produced an approximation of the energy density. After taking into account Cauchy’s equations and making mathematical transformations, it takes on the following form:
(
)
(
)
2 2 2 0 2 2 1 1 2 2 1 2 1 0,03 0,0128 8 x xy y x y xy xy x y Et U ε νγ ε νε ε ν ν ν γ γ ε ε − = + + + + − + + − + + (5)The expression above contains two additional terms, which are double underlined. It means that the values
of the energy density in the element with an additional rotational degree of freedom will be overestimated in relation to the ‘exact’ value. The phenomenon is known. It always occurs when degrees of freedom are added. In this case, a rigid node was added, which results in increased stiffness of the structure.
4. Conclusions
Energy criterion is not satisfied. That does
not, however, eliminate the element, because
the congruence condition is not obligatory, and
the overriding ellipticity condition is satisfied.
Applying reduced integration aimed to improve
the results by approximating the obtained energy
value to the exact solution. In the case under
consideration, however, such an approach does
not affect the results, which are the same as for
the exact solution.
References
[1] Gilewski G.: On Finite Element Assessment Criteria.
From Timoshenko Beam to Medium Thickness Plate (in Polish). Oficyna Wydawnicza Politechniki
Warszawskiej. 2005
[2] Laskowska U.: Transition elements with rotational
degrees of freedom (in Polish), Scientific Exercise
Books Silesian University of Technology, Building z.93, Gliwice 2001, pp. 269–276.
[3] Pawlak-Laskowska U.: Doctoral thesis: Non-standard
Finite Elements in Structure Modelling (in Polish).
Kielce 2004.
[4] Rakowski G.: New Frame – Membrane Transition
Element (in Polish). XLVI Konferencja Naukowa
Inżynierii Lądowej i Wodnej PAN i Komitetu Nauki PZITB, Wrocław-Krynica 2000.
[5] Rakowski G., Kacprzyk Z.: Finite Element Method in
Structure Mechanics (in Polish). Oficyna Politechniki
Warszawskiej. Warszawa 2005.
Katarzyna Suckert Urszula Pawlak
Weryfikacja poprawności sformułowania MES
prostokątnego elementu tarczowego z dodatkowym
rotacyjnym stopniem swobody
1. Wprowadzenie
Powszechność stosowania metody elementów skończonych w obliczeniach inżynierskich, szybki rozwój metod komputerowych oraz projektowanie coraz bardziej złożonych konstrukcji stwarza ko-nieczność opracowywania nowych elementów skoń-czonych oraz minimalizacji błędów obliczeń. Naj-częściej jest to realizowane na dwa sposoby: poprzez zagęszczenie siatki podziału lub przez wprowadzenie dodatkowych stopni swobody w elemencie.
W niektórych konstrukcjach żaden ze standardo-wych elementów nie odzwierciedla w sposób wy-starczający charakteru pracy konstrukcji. W związku z tym stale prowadzone są prace nad stworzeniem co-raz doskonalszych elementów, to znaczy elementów niestandardowych.
Elementy takie stosuje się kiedy mamy do czynie-nia z elementami różnego typu łączącymi się ze sobą,
np. tarczowy z prętowym. Wtedy we wspólnym węź-le (węzłach) pojawiają się dodatkowe stopnie swobo-dy. Przedmiotem dalszych rozważań będzie element tarczowy z dodatkowym rotacyjnym stopniem swo-body, przedstawiony na rysunku 1. W pracy dokona-no oceny poprawdokona-ności sformułowania wyżej wymie-nionego elementu.
2. Istota zagadnienia
W celu oceny poprawności elementu skończone-go analizie poddano macierz sztywności elementu tarczowego z dodatkowym rotacyjnym stopniem swobody. Dla przedstawionego powyżej elementu macierz sztywności wyprowadzona w sposób ścisły znajduje się w pracach [2-4].
Przedmiotem rozważań będzie macierz wypro-wadzona w sposób zredukowany. Koncepcję całko-wania zredukowanego, wykorzystuje się w analizie elementów skończonych, ponieważ takie podejście
niekiedy daje lepsze wyniki w stosunku do rozwią-zania ścisłego.
Analizowaną w pracy macierz sztywności otrzyma-no stosując czteropunktową kwadraturę Gaussa. Ma-cierz uzyskana w sposób zredukowany nieznacznie różniła się od macierzy ścisłej. Dla przykładu w tabe-li 1 zestawiono wyrażenia k11 uzyskane obydwoma
sposobami.
Macierz otrzymaną przy zastosowaniu całkowania zredukowanego sprawdzono wykonując weryfikację równań równowagi statyki oraz wykorzystując kryte-ria podane przez W. Gilewskiego [1], jakie poprawny element skończony powinien spełniać, tj.: warunek eliptyczności oraz warunek zgodności (kryterium energetyczne). Pierwszy warunek jest obligatoryjny i elementy nie spełniające go nie powinny być w ogóle stosowane, natomiast nie spełnienie warunku zgod-ności nie dyskwalifikuje elementu.
3. Sprawdzenie poprawności sformułowania macierzy sztywności elementu tarczowego z dodatkowym rotacyjnym stopniem swobody
3.1. Sprawdzenie równań równowagi statyki
Sprawdzenie to ma najbardziej podstawowy cha-rakter i jest najłatwiejsze, jednak umożliwia szybką weryfikację poprawności macierzy sztywności. Na-leży zwrócić uwagę, że niespełnienie tego warunku świadczy o błędnym sformułowaniu macierzy, nato-miast spełnienie nie gwarantuje prawidłowości i wy-maga dalszej weryfikacji. Zatem to sprawdzenie ma jedynie charakter orientacyjny.
W macierzy sztywności każda kolumna opisuje siły węzłowe w jednostkowym stanie przemieszczenia. Przykładając te siły do odpowiednich węzłów spraw-dza się równania równowagi. O poprawności macierzy w myśl tego kryterium świadczy zerowanie się sumy rzutów sił na oś x i y oraz momentu względem dowol-nego punktu. Dla analizowanej macierzy powyższy warunek spełniony został dla każdej z kolumn.
3.2. Sprawdzenie warunku eliptyczności
W celu sprawdzenia warunku eliptyczności należy znaleźć wartości i wektory własne macierzy sztyw-ności. Według [1] warunek eliptyczności jest jedno-znaczny ze spełnieniem następujących punktów:
1. Wartości własne macierzy sztywności muszą być nieujemne.
2. Spektrum własności własnych macierzy sztyw-ności musi zawierać dokładnie tyle samo zero-wych wartości własnych ile jest ruchów sztyw-nych elementów (LRSE).
3. Zerowym wartościom własnym muszą odpo-wiadać wektory własne, opisujące ruch sztywny. 4. Niezerowym wartościom własnym muszą odpo-wiadać wektory własne, które opisują odkształ-cone postacie elementu.
Aby sprawdzić powyższe zależności należy roz-wiązać zagadnienie własne (1).
Przy wykorzystaniu programu Mathematica i przyjęciu wymiarów elementu: a = 0.4 m, b = 0.2 m, t = 0.01 m, modułu sprężystości E = 192.6GPaoraz współczynnika Poissona υ = 0.3 obliczono wartości i wektory własne dla macierzy uzyskanej w sposób zredukowany. Dla analizowanej macierzy wszystkie wartości własne są określone dodatnio lub są zero-we, ponadto wartości zerowych jest 3 co odpowiada liczbie ruchów sztywnych dla elementów 2D. Zatem spełniony jest punkt 1 i 2, konieczny do spełnienia kryterium eliptyczności. Weryfikację dwóch pozo-stałych punktów umożliwi rysunek 2, ilustrujący wektory własne dla poszczególnych wartości wła-snych. Otrzymane wartości własne są następujące: λ={26,9158, 11,7113, 9,69056, 6,15919, 4,89254, 0,157527, 0, 0, 0}.
Zauważono, że wektory odpowiadające niezero-wym wartościom własnym (λ1– λ6) opisują odkształ-cone postacie elementu, natomiast wektory odpowia-dające zerowym wartościom własnym (λ7– λ9) opisują ruchy sztywne. Zatem spełnione są punkty 3 i 4, ko-nieczne do spełnienia kryterium eliptyczności.
Dodatkowo punkty 3 i 4 można sprawdzić anali-tycznie, obliczając odkształcenia odpowiadające wektorom własnym – εi (2).
Odpowiadające ruchom sztywnym wektory od-kształceń są zerowe, natomiast wektory opisujące postacie odkształcone są opisane wielomianami wyższego rzędu. Dla przykładu wektory ε1 ε8 mają
postać (3).
Podobnie jak przytoczony wektor ε1, każdy inny odpowiadający postaci odkształconej elementu po-siada człon, będący stałą, który zapewnia stałość od-kształceń.
Przeprowadzone obliczenia analityczne pozwalają wyciągnąć identyczne wnioski, jak dla przedstawionego rozwiązania graficznego na rysunku 2. Analiza wartości i wektorów własnych macierzy sztywności oraz odkształ-ceń jednoznacznie udowadnia, że spełnione są wszystkie 4 punkty odpowiadające kryterium eliptyczności. 3.3. Sprawdzenie warunku zgodności
Sprawdzenie to, zwane również kryterium energe-tycznym polega na przyrównaniu energii
nagroma-dzonej w rozpatrywanym elemencie wyliczonej ze sformułowania MES do rozwiązania dokładnego wg wzoru (4).
W celu sprawdzenia poprawności rozpatrywanego elementu w myśl kryterium energetycznego obliczo-no energię sprężystą względem jego środka ciężkości. W wyniku obliczeń wykonanych w programie Ma-thematica dla macierzy uzyskanej poprzez całkowanie zredukowane uzyskano przybliżoną postać gęstości energii, która po uwzględnieniu równań Cauchy’ego oraz dokonaniu przekształceń matematycznych przyjmuje postać (5).
Powyższe wyrażenie zawiera dwa dodatkowe czło-ny, które podkreślono podwójną linią. Oznacza to, że wartość gęstości energii w elemencie z dodatkowym rotacyjnym stopniem swobody będzie zawyżona w stosunku do wartości „dokładnej”. Zjawisko to jest znane, występuje zawsze przy dodawaniu stopni swobody [3]. W rozważanym przypadku ma to zwią-zek z tym, że dodano sztywny węzeł, co skutkuje zwiększeniem sztywności konstrukcji.
4. Wnioski
Kryterium energetyczne nie jest spełnione. Nie dyskwalifikuje to jednak elementu, gdyż warunek zgodności nie jest obligatoryjny, a nadrzędny waru-nek – eliptyczności został spełniony. Zastosowanie całkowania zredukowanego miało na celu polep-szenie otrzymanych wyników poprzez przybliżenie otrzymanej wartości energii do rozwiązania ścisłego. Jednak w rozważanym przypadku takie podejście nie ma wpływu na wyniki, są one identyczne jak dla roz-wiązania ścisłego.
