Reduced basis approaches for parametrized bifurcation problems held by non-linear Von Kármán equations
RReduced basis approaches for parametrizedbifurcation problems held by non-linear VonKármán equations
Federico PichimathLab, Mathematics Area, SISSA,International School for Advanced Studiesvia Bonomea 265, 34136 Trieste, ItalyGianluigi RozzamathLab, Mathematics Area, SISSA,International School for Advanced Studiesvia Bonomea 265, 34136 Trieste, ItalyJune 25, 2019
Abstract
This work focuses on the computationally efficient detection of the buck-ling phenomena and bifurcation analysis of the parametric Von Kármánplate equations based on reduced order methods and spectral analysis. Thecomputational complexity - due to the fourth order derivative terms, the non-linearity and the parameter dependence - provides an interesting benchmarkto test the importance of the reduction strategies, during the construction ofthe bifurcation diagram by varying the parameter(s). To this end, togetherthe state equations, we carry out also an analysis of the linearized eigenvalueproblem, that allows us to better understand the physical behaviour near thebifurcation points, where we lose the uniqueness of solution. We test thisautomatic methodology also in the two parameter case, understanding theevolution of the first buckling mode. a r X i v : . [ m a t h . NA ] J un Introduction and motivation
In this work we are interested in the numerical approximation of parameter depen-dent non-linear structural problems governed by the Von Kármán plate equations[36]. The main issue, and so the most interesting feature, of this kind of model,is the bifurcation phenomenon that corresponds to the buckling behaviour of anelastic thin plate.From the mathematical point of view, we dealt with a problem in whichmultiple solutions for the same value of the parameter can arise. This led us tothe non-uniqueness of the solution that we investigated. Numerically speaking themodel presents three main problems that we had to face with: (i) non-linearity (ii)parameter dependency (iii) high order derivative terms.Thus we relied on the Reduced Basis (RB) method aiming at reducing the com-putational time in order to have a better understanding of the physical phenomenathat we were modelling. In practice we wanted to find an efficient and rapid wayto draw the bifurcation diagram and detect the buckling points, i.e. critical valuesof the parameter λ , which in our case controls the compression along the edges ofthe plate.The RB method is one of the main Reduced Order Modelling (ROM) techniquesand we applied it in conjunction with the Galerkin-Finite Element (FE) method[27, 20, 31, 21]. As the latter, the RB approach is a Galerkin projection over afinite dimensional subspace of the weak formulation of the model. Thus, we solvethe FE problem, called high order formulation because of the huge number ofdegrees of freedom involved, then we construct a subspace as the span of somebasis functions computed before, and finally we project again on this new smallersubspace. A recent review chapter focused on parametric elasticity problemssolved by RB method is [22].This technique allows us to efficiently study the entire behaviour of the plateunder compression, but still involves a huge amount of computations in the offlinehigh order phase. Thus, being inspired by the recent works in branching detection[19, 29], also in a more industrial framework [17], we supplemented the modelwith the eigenvalue analysis of the linearized equations.Indeed, in the past, many authors established the connection between thebuckling points and the eigenvalues behaviour, but the computational complexityof the problem was not affordable for that time [6, 8]. Hence we propose a newreduced order approach, in order to avoid the huge computational cost. Previousworks on model order detection for non-linear elasticity could be found in [34, 37],as well as in preliminary works by Noor and Peters [25, 24, 26].2he structure of the work is the following. In Section 2 we provide a briefdescription of the equations of plate (2 . .
2) of the problem, necessary as the firststep towards the numerical approximation. To end, we recall few definitions andtheorems in bifurcation and non-linear analysis (2 .
3) that justify the eigenproblemcoupling.In Section 3 we deal with the numerical approximation of the problem, provid-ing the pseudo-code used to construct the bifurcation diagram, based on Galerkinfinite element (3 .
1) and Newton method. The reduction strategy, or reduced basismethod is reported in (3 .
2) with its main features and its matrix formulation, com-pared with the finite element one. Finally, we show some preliminary results onthe spectral analysis (3 . Starting from the very well known theory of continuum mechanics, Von Kármánin 1910 proposed a mathematical model in order to describe all the possibleconfigurations that a plate under compression can take [36].
Buckling phenomenon is the mathematical way of explaining a well known physical event that veryfrequently happens in many contexts. As an example, an appropriate one since it3s exactly what we are trying to model, if we pick a thin rectangular plate at rest,we can use our hands to compress it until we reach a critical point, i.e. when theplate takes a deformed configuration, or it buckles.
Let us consider an elastic and rectangular plate Ω = [0 , L ] × [0 ,
1] in its undeformedstate, subject to a λ -parametrized external load acting on its edges, then thedisplacement from its flat state and the Airy stress potential , respectively u and φ ,satisfy the Von Kármán equations (cid:40) ∆ u = [ λ h + φ, u ] + f , in Ω∆ φ = − [ u , u ] , in Ω (1)where h and f are some given functions, that we can set to specify the exter-nal forces acting on our plate, while ∆ is the biharmonic operator in Cartesiancoordinates and [ u , φ ] := ∂ u ∂ x ∂ φ∂ y − ∂ u ∂ x ∂ y ∂ φ∂ x ∂ y + ∂ u ∂ y ∂ φ∂ x , is the brackets of Monge-Ampére . Thus we aim to find the displacement and thecoupled Airy stress potential that solve the system (1) which is of fourth order, dueto the presence of the biharmonic operator, non-linear due to the product of secondderivatives in the bracket, and parametric due to the buckling coefficient λ varyingin a proper range of real numbers. Moreover we are presenting, for the sake ofsimplicity, a non-dimensional model where all the physical quantities, except forthe compression parameter λ , are set to unity.In order to have a well posed system of partial differential equations we pro-vide boundary conditions for both the unknowns that match the different physicalsetting. Among all the possible choices, we just focused on the so called simplysupported boundary conditions (cid:40) u = ∆ u = 0 , in ∂ Ω φ = ∆ φ = 0 , in ∂ Ω (2)which are physically complex to reproduce, but also the most used ones for thesimulations because of their versatility and importance also in the weak formu-lation. So from now on, unless specified otherwise, we will consider the system41) with simply supported boundary conditions (2). We remark that despite thesimple boundary conditions chosen, the goal of this work is understanding thebifurcation behaviour for a complex system, regardless the numerical constraintsthat a conforming method for more involved boundary condition could impose.Thanks to the BCs chosen we can split the system of two fourth order non-linearelliptic equations, into a system of four second order non-linear elliptic equations.In order to carry out this trick, given by Ciarlet-Raviart [16], we introduce twonew unknowns, namely U = ∆ u , Φ = ∆ φ , so that we can rewrite the system (1) withhomogeneous Dirichlet boundary conditions as ∆ U = [ λ h + φ, u ] + f , in Ω∆ u = U , in Ω∆Φ = − [ u , u ] , in Ω∆ φ = Φ , in Ω with u = 0 , in ∂ Ω U = 0 , in ∂ Ω φ = 0 , in ∂ ΩΦ = 0 . in ∂ Ω (3)We know that (1) and (3) are equivalent [38] when the boundary is regular andthe solution is smooth enough, so from now on we just consider the latter.Finally, since we are interested in the behaviour of the plate under compression,we can set the external body force f = 0 and model different kind of stresses atthe boundaries through the function h . Indeed if we choose h = − y , we obtainthat [ λ h , u ] = − λ u xx where we are assuming that the compression is acting on theedges parallel to the y direction (see Figure 1). Note also that if instead we choose h = − ( x + y ) we would have the stress component given by [ λ h , u ] = − λ ∆ u , inwhich case the compression we are considering is on the whole boundary.Ω z y x λ λ Figure 1: A rectangular bi-dimensional elastic plate compressed on its edges5 .2 Weak formulation
Starting from the considerations of the previous section we are now able to set ourproblem in a more abstract mathematical framework, which will be considered byus for the following numerical investigation. So let us consider λ ∈ D ⊂ R , where D is the parameter space, Ω ⊂ R is the bi-dimensional domain, that we identifywith the plate, whereas V = V (Ω) = (cid:16) H (Ω) (cid:17) is the Hilbert space in which we willseek the solution and V (cid:48) its dual space.Furthermore we represent the non-linear PDE with the parametrized mapping G : V × D → V (cid:48) , so that the Von Kármán system (3), in abstract form reads: given λ ∈ D , find X ( λ ) . = ( u ( λ ) , U ( λ ) , φ ( λ ) , Φ( λ )) ∈ V such that G ( X ( λ ); λ ) = 0 , in V (cid:48) . (4)We present the weak formulation where all the boundary terms vanish due tothe simply supported boundary conditions (2). Then, we seek X ( λ ) ∈ V such that ( ∇ u , ∇ w ) L (Ω) + ( U , w ) L (Ω) = 0 , ∀ w ∈ H (Ω)( ∇ U , ∇ v ) L (Ω) + ( [ λ h + φ, u ] , v ) L (Ω) = 0 , ∀ v ∈ H (Ω)( ∇ φ, ∇ θ ) L (Ω) + (Φ , θ ) L (Ω) = 0 , ∀ θ ∈ H (Ω)( ∇ Φ , ∇ ψ ) L (Ω) − ( [ u , u ] , ψ ) L (Ω) = 0 , ∀ ψ ∈ H (Ω) (5)in which we embed the simply supported boundary conditions in the choice of thespace H (Ω), where the test functions Y . = ( w , v , θ, ψ ) reside. Moreover we denotewith ( · , · ) L (Ω) the usual inner product in the Hilbert space L (Ω).Then coming back to the abstract form of our problem (4) the weak formulationreads: given λ ∈ D , find X ( λ ) ∈ V such that g ( X ( λ ) , Y ; λ ) = 0 , ∀ Y ∈ V , (6)where the parametrized variational form g ( · , · ; λ ) : V × V → R is defined as g ( Z , Y ; λ ) = (cid:104) G ( Z ; λ ) , Y (cid:105) , ∀ Z , Y ∈ V , where we denoted the duality pairing between V (cid:48) and V with (cid:104)· , ·(cid:105) . In this case the parametrized variational form g ( · , · ; λ ) is defined asfollows (7) g ( X ( λ ) , Y ; λ ) = a ( u ( λ ) , w ) + b ( U ( λ ) , w ) + a ( U ( λ ) , v ) + λ c ( h , u ( λ ) , v )+ c ( φ ( λ ) , u ( λ ) , v ) + a ( φ ( λ ) , θ ) + b (Φ( λ ) , θ )+ a (Φ( λ ) , ψ ) − c ( u ( λ ) , u ( λ ) , ψ ) , ∀ Y ∈ V , ∀ λ ∈ D , a ( x , y ) = (cid:90) Ω ∇ x · ∇ y d Ω , b ( x , y ) = (cid:90) Ω x y d Ω , c ( x , y , z ) = (cid:90) Ω [ x , y ] z d Ω . The numerical treatment of the variational form including the bracket of Monge-Ampére obviously needs a non-linear method, to this end we compute the partialFréchet derivative of g ( Z , · ; λ ) with respect to X at Z ∈ V .Thus we assume the mapping G to be continuously differentiable and denoteby D X G ( Z ; λ ) : V → V (cid:48) its partial Fréchet derivative at ( Z , λ ) ∈ V × D . In this waywe can express the partial Fréchet derivative of g ( Z , · ; λ ) at Z ∈ V as d g [ Z ]( W , Y ; λ ) = (cid:104) D X G ( Z ; λ ) W , Y (cid:105) , ∀ W , Y ∈ V . (8)These computations, for the non-linear system we are considering, show theexplicit expression of the derivative of g be of the form (9) d g [ Z ]( X ( λ ) , Y ; λ ) = a ( u ( λ ) , w ) + b ( U ( λ ) , w ) + a ( U ( λ ) , v ) + λ c ( h , u ( λ ) , v )+ c ( φ ( λ ) , Z , v ) + c ( Z , u ( λ ) , v ) + a ( φ ( λ ) , θ )+ b (Φ( λ ) , θ ) + a (Φ( λ ) , ψ ) − c ( u ( λ ) , Z , ψ ) − c ( Z , u ( λ ) , ψ ) , ∀ Z , Y ∈ V , ∀ λ ∈ D , where we denoted with Z = ( Z , Z , Z , Z ) the components of the point in whichwe are computing the derivative. The focus of our work is the efficient detection of the possible multiple solutions ofthe equations (3). Following the work done in [2, 11, 9] we recall the mathematicaldefinitions of bifurcation theory, which, as we will see, will serve us in developinga tool for the efficient detection of the buckling points.Since the undeformed configuration is a trivial solution for every λ ∈ R , i.e. G (0 , λ ) = 0, we can denote with S = { ( X , λ ) ∈ V × R : X (cid:54) = 0 , G ( X , λ ) = 0 } , (10)the set of non-trivial solutions of (4), then we can finally define the bifurcationpoints. 7 efinition 1 We say that λ ∗ ∈ R is a bifurcation point for G : V × D → V (cid:48) , from thetrivial solution, if there is a sequence ( X n , λ n ) ∈ V × R with X n (cid:54) = 0 and G ( X n , λ n ) = 0 such that ( X n , λ n ) → (0 , λ ∗ ) . In order to understand where are located these bifurcation points, we couldnumerically investigate the equations for each value of the parameter λ observingwhen the buckling phenomena occur. Of course this way is computationally veryexpensive and we need a different tool for the detection.To this end we have analyzed the path pursued in [8, 2, 6] where the link be-tween the bifurcation points and the behaviour of the eigenvalues of the linearizedproblem is highlighted. If we linearize the equations (3) around the trivial solutionthe system we obtain is simply given by (cid:40) ∆ U = λ [ h , u ] , in Ω∆ u = U , in Ω , with (cid:40) u = 0 , in ∂ Ω U = 0 , in ∂ Ω . (11)This connection is not surprising, in fact from ODE’s theory we know that thestability of the solutions is linked to eigenvalues that change sign, i.e. cross theimaginary axis varying λ . Now we briefly recall the main theorems [2], basedon an application of the Implicit Function Theorem, that validate the numericalinvestigation we carried out. Theorem 1
A necessary condition for λ ∗ to be a bifurcation point for G : V × D → V (cid:48) is that the partial derivative D X G (0; λ ∗ ) is not invertible. Theorem 2
Bifurcation points of Von Kármán equations (3) with respect to thetrivial solution, i.e. X = 0 , can occur only at the eigenvalues of the linearizedproblem (11) . The former is a general result, while for Von Kármán equations we also know [8, 6]that every bifurcation point is an eigenvalue of the linearized problem. Moreover,if we assume that all the eigenvalues are real, positive and ordered in such a waythat 0 ≤ λ ≤ λ ≤ λ ≤ . . . , from [8] we can reverse the statement. Theorem 3
From each eigenvalue of the system (11) at least one branch of non-trivial solution of (1) bifurcates. In particular from a simple eigenvalue onebranch bifurcates and from a multiple eigenvalue at least two branches bifurcate.Furthermore if λ is the smallest eigenvalue then for every value λ ≤ λ the unforcedsystem has no non-trivial solutions. Numerical approximation of the problem
As we said, in this application we have to face with different kinds of difficulties,such as: non-linearity, fourth order derivatives, and parameter dependence.This led us to use a double approach, the full order and the reduced order onethat we will present soon in detail. Let us present now, through a pseudo-code,how we addressed all these issues in order to obtain the bifurcation diagram.
Algorithm 1
A pseudo-code setting while λ < λ f do (cid:46) External loop on compression parameter if || u h || ∞ < δ then (cid:46) Pre-buckling initial guess X h ( λ ) = X (0) h end if while || δ X h || V > (cid:15) do (cid:46) Newton method D X G ( X ( i ) h ( λ ); λ ) δ X h = G ( X ( i ) h ( λ ); λ ) (cid:46) Galerkin FE method X ( i +1) h ( λ ) = X ( i ) h ( λ ) − δ X h end while λ = λ + d λ end while The Algorithm 1 is the implementation result of our approach to the difficul-ties previously mentioned. We denote by X h = ( u h , U h , φ h , Φ h ) the approximatedsolution and by D = [ λ i , λ f ] the parameter domain and we start with a slightlymodified continuation method [1], which in its basic formulation consists in awhile loop over λ ∈ D . So, at each new cycle, i.e. for every λ , we are changingthe initial guess for the non-linear problem from a properly chosen guess (usuallythe solution of the linearized problem) to the solution of the last cycle, as soon asthe solution of the latter turns out to be non-trivial. Thus we are able to detect thebuckling by looking at the behaviour of the maximum norm of the displacement.The second goal we had was to overcome the issue of the non-linearity. Wechoose the very well know Newton-Kantorovich method [15], where we use thenorm of the increment δ X h in the Hilbert space V as stopping criterion.Finally we note that in the line 6 we have to find the solution of a new, linear,weak formulation. To this end we applied the Galerkin-Finite Element method which appears to be a good candidate, also in view of the numerical extensiontowards the model order reduction. Since the fundamental importance of the lasttwo methods, we briefly recall in the next section how we applied them to theVon Kármán equations. 9 .1 Galerkin finite element method
As we already said, we decided to use the Galerkin finite element method todiscretize the problem. This is a projection-like technique, where using the versa-tility of the weak formulation, we can set the model in a finite dimensional space[14, 30].So we again consider the weak problem (6) and denote with V h a family ofspaces, dependent from the h parameter, such that V h ⊂ V , dim ( V h ) = N h < ∞ forall h >
0. Then, given λ ∈ D , we seek X h ( λ ) ∈ V h that satisfies g ( X h ( λ ) , Y h ; λ ) = 0 , ∀ Y h ∈ V h . (12)Obviously, because of the non-lineartity, we can not directly apply the finiteelement method to the equation (12). So the Newton method, which in this case isalso known as Newton-Kantorovich [15, 31], reads as follows: once assigned aninitial guess X h ( λ ) ∈ V h , for every k = 0 , , . . . we seek the variation δ X h ∈ V h suchthat d g [ X kh ( λ )]( δ X h , Y h ; λ ) = g ( X kh ( λ ) , Y h ; λ ) , ∀ Y h ∈ V h , (13)and then we update the solution for the successive step as X k +1 h ( λ ) = X kh ( λ ) − δ X h ,until we reach the convergence with the stopping criteria we discussed before.From the algebraic point of view we denote with { E j } N h j =1 a base for the space V h such that we can write every element X h ( λ ) ∈ V h as X h ( λ ) = N h (cid:88) j =1 X ( j ) h ( λ ) E j , (14)so that we obtain the solution vector (cid:174) X h ( λ ) = { X ( j ) h ( λ ) } N h j =1 . We then lead back tothe study of the solution (cid:174) X h ( λ ) ∈ R N h of the system g (cid:32) N h (cid:88) j =1 X ( j ) h ( λ ) E j , E i ; λ (cid:33) = 0 , ∀ i = 1 , . . ., N h (15)which, recalling the notations we introduced in Section 1 .
2, corresponds tothe solution of G h ( (cid:174) X h ( λ ); λ ) = 0, where the residual vector G h is defined as( G h ( (cid:174) X h ( λ ); λ )) i = g ( X h ( λ ) , E i ; λ ). Finally, the Newton method combined withthe Galerkin finite element method, applied to our weak formulation reads: find δ (cid:174) X h ∈ R N h such that J ( (cid:174) X kh ( λ ); λ ) δ (cid:174) X h = G h ( (cid:174) X kh ( λ ); λ ) , (16)10here we defined the Jacobian matrix in R N h × N h as J ( (cid:174) X kh ( λ ); λ )) i j = d g [ X kh ( λ )]( E j , E i ; λ ) , for all i , j = 1 , . . ., N h . (17)Moreover we use standard Lagrange finite element so that V h = ( ˚ X rh ) where X rh = { Y h ∈ C ( ¯Ω) : Y h | K ∈ P r , ∀ K ∈ T h } (18)and ˚ X rh = { Y h ∈ X rh : Y h | ∂ Ω = 0 } (19)is the space of globally continuous functions that are polynomials of degree r on the single element of the triangulation T h of the domain, which vanish on theboundary.To provide a clear matrix representation of the application of the Galerkinmethod, we present the projected weak formulation. In this case the Newtonmethod (16) reads: given an initial guess X h = ( u h , U h , φ h , Φ h ) ∈ V h for k = 0 , , . . . until convergence we seek δ X h = ( δ u h , δ U h , δφ h , δ Φ h ) ∈ V h such that a ( δ u h , w h ) + b ( δ U h , w h ) = a ( u kh , w h ) + b ( U kh , w h ) , ∀ w h ∈ ˚ X rh a ( δ U h , v h ) + c ( δφ h , u kh , v h ) + c ( φ kh , δ u h , v h ) + λ c ( h , δ u h , v h ) = a ( U kh , v h ) + c ( φ kh , u kh , v h ) + λ c ( h , u kh , v h ) , ∀ v h ∈ ˚ X rh a ( δφ h , θ h ) + b ( δ Φ h , θ h ) = a ( φ kh , θ h ) + b (Φ kh , θ h ) , ∀ θ h ∈ ˚ X rh a ( δ Φ h , ψ h ) − c ( δ u h , u kh , ψ h ) − c ( u kh , δ u h , ψ h ) = a (Φ kh , ψ h ) − c ( u kh , u kh , ψ h ) , ∀ ψ h ∈ ˚ X rh (20)and then set X k +1 h = X kh − δ X h . We can finally present the matrix formulation thatfollows directly from (16) and (20) (cid:169)(cid:173)(cid:173)(cid:173)(cid:171) A h B h C h + λ C h A h C h
00 0 A h B h − C h − C h A h (cid:170)(cid:174)(cid:174)(cid:174)(cid:172) (cid:169)(cid:173)(cid:173)(cid:173)(cid:171) δ u h δ U h δφ h δ Φ h (cid:170)(cid:174)(cid:174)(cid:174)(cid:172) = (cid:169)(cid:173)(cid:173)(cid:173)(cid:171) A h u kh + B h U kh A h U kh + C h u kh + λ C h u kh A h φ kh + B h Φ kh A h Φ kh − C h u kh (cid:170)(cid:174)(cid:174)(cid:174)(cid:172) , (21)where we denoted the matrices as follows( A h ) i j = a ( E j , E i ) , ( B h ) i j = b ( E j , E i ) , ( C h ) i j = c ( h , E j , E i ) , ( C h ) i j = c ( E j , u kh , E i ) , ( C h ) i j = c ( φ kh , E j , E i ) , ( C h ) i j = c ( u kh , E j , E i ) . Note that because of the symmetry of the bracket of Monge-Ampére, we easilyobtain that it holds C h ≡ C h . 11 .2 Reduced Basis method Dealing with the approximation of a parametrized problem could be very difficult,sometimes we must therefore rely on some techniques, by which we can reduce thecomputational cost. With this aim in the past years many authors, to mention fewworks [27, 20, 31, 7], developed and applied the reduced order methods (ROM), acollection of methodologies used to replace the original high dimension problem,called high fidelity approximation , with a reduced problem that is easy to manage.One of these methodologies is the
Reduced Basis method (RB), that consistsin a projection of the high fidelity problem on a subspace of smaller dimension,constructed with some properly chosen basis functions.At the beginning this method was used for non-linear structural problems[25, 26], but the computational complexity of this kind of equations was too big toprovide a deep understanding of the bifurcation phenomena that we have analyzed.As we have seen, the preliminary step is the projection of the weak formulation(6) in a discretized setting, which results in the Galerkin problem (12). Findinga numerical solution to this problem is very challenging because of the potentialhigh number of degrees of freedom N h .Thus we aim at building a discrete manifold V N induced by properly chosensolution of (12) and then project over it. This is the description of the first step,namely the offline phase , in which we explore the parameter space D in orderto construct a basis for the reduced space of dimension N . On the other side,the second step, called online phase , is the efficient and reliable part where thesolutions are computed through the projection on V N . This complexity reduction isbased on two main key points: the assumption that it holds the affine decompositionand the fact that N (cid:28) N h .As in the offline phase (12), given λ ∈ D , we seek X N ( λ ) ∈ V N that satisfies g ( X N ( λ ) , Y N ; λ ) = 0 , ∀ Y N ∈ V N . (22)At this point we have again to face with the non-linearity, and so come back tothe Newton-Kantorovich method obtaining: given an initial guess X N ( λ ) ∈ V N , forevery k = 0 , , . . . we seek the variation δ X N ∈ V h such that d g [ X kN ( λ )]( δ X N , Y N ; λ ) = g ( X kN ( λ ) , Y N ; λ ) , ∀ Y N ∈ V N , (23)and then we update the solution as X k +1 N ( λ ) = X kN ( λ ) − δ X N until convergence.From the algebraic point of view and thus to do another step towards thereduced solution, we introduce the orthonormal base { Σ m } Nm =1 for the space V N X N ( λ ) = N (cid:88) m =1 X ( m ) N ( λ )Σ m , (24)and denote with (cid:174) X N ( λ ) = { X ( m ) N ( λ ) } Nm =1 ∈ R N the reduced solution vector.Choosing properly the test element Y N ∈ V N as Y N = Σ n for every 1 ≤ n ≤ N ,we obtain the algebraic system in R N given by( G N ( (cid:174) X N ( λ ); λ )) n . = g (cid:32) N (cid:88) m =1 X ( m ) N ( λ )Σ m , Σ n ; λ (cid:33) = 0 , ∀ n = 1 , . . ., N , (25)where we denoted with G N ( (cid:174) X N ( λ ); λ ) the residual reduced vector and with V thetransformation N h × N matrix whose elements ( V ) jm = Σ m ( j ) are the nodal evaluationof the m th basis function at the j th node. Moreover, we note that (25) correspondsto the solution of V T G h ( V (cid:174) X N ( λ ); λ ) = 0 . (26)Finally we can apply again the Newton method, which combined with thereduced basis method, provides the following formulation : find δ (cid:174) X N ∈ R N suchthat J N ( (cid:174) X kN ( λ ); λ ) δ (cid:174) X N = G N ( (cid:174) X kN ( λ ); λ ) , (27)where J N is the reduced Jacobian R N × N matrix defined as J N ( (cid:174) X kN ( λ ); λ ) = V T J ( V (cid:174) X kN ( λ ); λ ) V . (28)Thus we want to construct the reduced problem through the projection on asubspace V N ⊂ V h from a collection of the so called snapshots , i.e. the solutions ofthe full order problem for specific values of the parameter selected by a samplingtechnique. The most famous strategies to construct V N are the Proper orthogonaldecomposition (POD) and the
Greedy algorithm [20, 27, 31]. In this work werelied on the former, based on an ordered sampling of the interval D , since itsphysical interpretation w.r.t the energy of the problem and because of the lack ofa rigorous a posteriori error estimate for the latter.Obviously POD increases the computational cost in the offline part, but si-multaneously gives us a reliable representation of the reduced manifold and keeptrack of the energy information that we are discarding. Moreover, we can slightlymodify it in order to consider multi parameter case, as the one in Section 4 . V N = span { Σ n , n = 1 , · · · , N } where { Σ n } Nn =1
13s the basis functions set obtained through the
Gram-Schmidt orthonormalizationprocedure.Now we show how the reduced basis method reflects the projection propertiesof the Galerkin method also in the online phase. Indeed the weak formulation thatwe obtain from the application of the Newton method at the reduced level reads as(23), with the reduced Jacobian J N ( (cid:174) X kN ( λ ); λ ) ∈ R N × N having the same structure ofthe finite element one J N ( (cid:174) X kN ( λ ); λ ) = (cid:169)(cid:173)(cid:173)(cid:173)(cid:171) A N B N C N + λ C N A N C N
00 0 A N B N − C N − C N A N (cid:170)(cid:174)(cid:174)(cid:174)(cid:172) , (29)where, if we introduce the transformation matrices with respect to the differentcomponents of the solution, V u and V φ , respectively for u and φ , we can definethe reduced matrices in the following way: C N = V T C h V , C N = ˜ N (cid:88) n =1 u ( n ) N V Tu C h (Σ n ) V u , C N = ˜ N (cid:88) n =1 φ ( n ) N V T φ C h (Σ n ) V φ , C N = ˜ N (cid:88) n =1 φ ( n ) N V T φ C h (Σ n ) V φ , A N = V T A h V , B N = V T B h V . Moreover, we highlight that also the reduced residual vector has the same formof the high order one, indeed it reads G N ( (cid:174) X kN ( λ ); λ ) = (cid:169)(cid:173)(cid:173)(cid:173)(cid:171) A N u kN + B N U kN A N U kN + C N u kN + λ C N u kN A N φ kN + B N Φ kN A N Φ kN − C N u kN (cid:170)(cid:174)(cid:174)(cid:174)(cid:172) . (30)We have illustrated the online phase, that permits an efficient evaluation of thesolution and possibly related outputs for every possible choice of a different param-eter λ ∈ D . The key point of this time saving is the so called affine decomposition .Indeed we want the computations to be independent form the, usually very high,number of degrees of freedom N h of the true discrete problem . In general thereduced matrices we have just presented are λ -dependent and an affine-recoverytechnique called Empirical Interpolation Method (EIM) is needed [5].14 .3 Spectral analysis In the previous sections we discussed about the issue of the computational com-plexity of the problem itself, that we try to avoid using the ROM. It is clear thatdrawing the bifurcation diagram is yet a difficult task. Indeed how can we inves-tigate the parameter space D without having any information on the position ofthese points?Taking some inspiration from [29, 28], where the stability property is analyzedwith the help of the spectral problem, supported by the theoretical results given inSection 1 .
3, we tried in this way to locate more precisely the buckling points.Thus we construct the eigenvalue problem for the linearized parametrizedoperator (cid:40) ∆ u + λ u xx = σ λ u , in Ω = [0 , L ] × [0 , u = ∆ u = 0 , in ∂ Ω (31)where we want to find, varying the buckling parameter, the couple ( u , σ λ ) ∈ H (Ω) × R , whose components are respectively the eigenfunction and eigenvalue. We willrestrict our simulations to the square plate with L = 1 and the rectangular one with L = 2.We are interested in the behaviour of σ λ with respect to λ , in fact since thesign of the eigenvalues is strictly linked with the stability property of the solution,we aim at observing that the first eigenvalue crosses the y-axis when the plate isbuckling. This is exactly what we found, indeed in Figure 2 for L = 1 we can seethe behaviour of the first four eigenvalues σ λ for λ ∈ [ , ] and if we use a λ -stepequal to one half the crossing happens for the value λ = 39 . λ ∈ [30 , D we alsoobserved the crossing of the successive eigenvalues.Finally, if we solve the eigenvalue problem (31) for the case of L = 2, wenote that a simple eigenvalue of the square plate becomes a multiple eigenvaluewith algebraic multiplicity equal to two (see Figure 3). This fact has a relevantconsequence from the physical point of view as we will see in the bifurcationdiagram later.Figure 3: Double eigenvalue for the rectangular plate crossing for λ = 62.What we just showed is computationally heavy to perform, so to keep in mindthe efficiency as key word of the whole analysis, we tried two different ways thatvalidate the result and at the same time reduce the computational time.Thus we consider the linear problem (11) but in its original form (cid:40) ∆ u + λ u xx = 0 , in Ω = [0 , L ] × [0 , u = ∆ u = 0 , in ∂ Ω (32)that has non trivial solutions for m , n = 1 , , . . . given by u m , n = sin (cid:16) m π xL (cid:17) sin ( n π y ) if and only if λ m , n = (cid:16) π L (cid:17) (cid:20) m + n L m (cid:21) , (33)where u m , n and λ m , n can be considered as the eigenfunctions and eigenvalues forthis new generalized eigenvalue problem. Now we have a simpler problem, indeedsince λ is now the eigenvalue, there is no parametrization. This provide us also anexplicit expression for the spectra, that we can use to validate the results.16=1 h = 1.e-1 h = 6.e-2 h = 1.e-2 h = 6.e-3 Order Exact λ , λ , λ , L = 1 : λ , = 4 π , λ , = 254 π , λ , = 1009 π , λ , = 28916 π , L = 2 : λ , = 4 π , λ , = 16936 π , λ , = λ , = 254 π , indeed we obtain the value λ , (cid:39) .
47 predicted in Figure 2 for the squareplate, while we note the presence of the double eigenvalue λ , = λ , (cid:39) .
68 thatconfirms what we saw in Figure 3 for the rectangular one.Finally, using the techniques in [4, 23], is an easy task to prove the followingtheorem that provides us a tool to better understand how good is our approximation.
Theorem 4
There exists a strictly positive constant C such that | λ − λ h |≤ C h , where λ h is an approximation, dependent on the sparsity of the grid, of the trueeigenvalue λ . The theorem above is crucial when we are dealing with problems for which wedo not know an explicit expression of the eigenvalues. For the sake of completenesswe provide in Table 1 and Table 2 the order of convergence results respectively forthe square and rectangular plate, that agree with the theoretical ones.To conclude this section we briefly discuss also the second straightforwardway to reduce the computational complexity of solving multiple times a full ordereigenvalue problem. Coming back to the parametrized eigenproblem (31), we canapply again the Reduced Basis method, and thanks to the affine decomposition, weeasily obtain in a more efficient way the same behaviour of the results discussed17=2 h = 1.e-1 h = 6.e-2 h = 1.e-2 h = 6.e-3 Order Exact λ , λ , λ , λ , λ , in the full order case.Figure 5: First eigenvalue λ , in the reduced order case.before, as we can see from Figures 4 and 5. We do not discuss further this lastapproach, in fact we can embed the computation for the eigenproblem (32) in theoffline phase. In this Section we will show how the buckling phenomena, i.e. the loss ofuniqueness of the solution, appears in the equation through the bifurcation diagram,both in the square and rectangular plate cases.Thus, in order to recap, let us consider the Von Kármán plate equations, withsimply supported BCs, in the bi-dimensional domain Ω = [0 , L ] × [0 ,
1] given by18 ∆ u + λ u xx = [ φ, u ] , in Ω∆ φ = − [ u , u ] , in Ω u = ∆ u = 0 , in ∂ Ω φ = ∆ φ = 0 , in ∂ Ω (34)and we are interested in the study of the solution while varying the bucklingparameter λ , which describes the compression along the edges parallel to the y -axis.This model was previously numerically investigated by many authors [10, 13,33], but as we already said the biggest issue was the computational complexity,that we overcame by means of the Reduced Basis method. We performed all thesimulations within FEniCS [3] for the full order case and RBniCS [32] for thereduced order one.The investigation done with the eigenproblem give us the necessary informationthat the parameters responsible of the buckling live in the interval D = [35 , We present the bifurcation diagram in Figure 6 for the square plate case Ω =[0 , × [0 , λ ∈ D on the x -axis, thecorrespondent value of the full order displacement u in its point of maximummodulo. How we predicted previously, we can observe the buckling phenomenafrom the trivial solution. Moreover, we note that the first bifurcation happens for λ value near λ , (cid:39) .
47. We did not stop at the first bifurcation, in fact choosingproperly the initial guess we have been able to detect also the second bifurcationfor the square plate. This result is confirming what we predicted, since for λ near λ , = 61 .
68 we obtain other branches.The physical symmetry issue is evident in both buckling points and the sameholds for the rectangular plate. Indeed, once we chose a bifurcation point, asolution from the upper branch is the same solution of the other one, but reflectedwith respect to the plate plane. Thus, for the first bifurcation near λ , , looking atthe contour plot, we observe a one cell like displacement as in Figure 7. Whileif we look at the second branch, so the one near λ , , we find a two cells likedisplacement as in Figure 8. So for the square plate case we obtained four differentsolutions for each λ ≥ λ , . 19oreover, the RB method worked well with this problem. Indeed as we cansee in Figure 9, the reduced basis solution approximates perfectly not only thebehaviour but also the order of magnitude. The remarkable point is that in orderto obtain the solution on the right in Figure 9 we just solved a linear system ofdimension 5 instead of the one given by the Galerkin full order method of order8 · . Figure 6: Bifurcation diagram for the square plate.Figure 7: Full order one cell so-lution for the displacement u with λ = 65 (green branch). Figure 8: Full order two cells so-lution for the displacement u with λ = 65 (blue branch).20igure 9: Comparison between the full order solution (on the left) and the reducedorder one (on the right) for the displacement u with λ = 46, belonging to the greenbranch. Below the reduced basis error plot.We just saw that the Reduced Basis method provides us a useful technique,always based on a full order method (in this case is the Galerkin-Finite Element),that allows us to obtain the same results, at the cost of a small error, but with a hugeamount of time saving. To be more precise we present in Table 3 a convergenceresults: the error between the truth approximation and the reduced one as a functionof N . The error reported, E N = max λ ∈D || u h ( λ ) − u N ( λ ) || H (Ω) is the maximum ofthe approximation error over a uniformly chosen test sample.We highlight that we present here just the full order bifurcation diagrams since,also in view of Table 3 and Figure 9, the reduced order one is exactly the same.We remark that we did not implemented the Greedy algorithm here, becausea suitable extension of Brezzi-Rappaz-Raviart (BRR) theory for the a posteriorierror estimate would be needed [9, 35, 18, 12]. However, applying BRR theory atreduced level is not straightforward and we leave it for further future investigation.Finally, as regards computational times, a RB evaluation λ → u N ( λ ) requiresjust t RB = 100 ms for N = 8; while the FE solution λ → u h ( λ ) requires t FE = 8 . .
22% of the FEM computational cost.21 E N E +002 6.90 E -013 7.81 E -024 2.53 E -025 1.88 E -026 1.24 E -027 9.02 E -038 8.46 E -03Table 3: The reduced basis convergence with respect to the number of the basis Nfor the square plate case.Figure 10: Bifurcation diagram for the rectangular plate. Now we analyze the case of the rectangular plate, where the domain is Ω =[0 , × [0 , L as a second geometrical parameter along with λ .In this case, the solution start branching, as before, from the first orderedeigenvalue λ , (cid:39) .
47. Obviously the number of the cells in the contour plot thatare formed strictly depends on the length of the domain, for example here the firstbifurcation is linked with the two cell configuration as shown in Figure 11.We observed also a new bifurcation for λ value near 46 .
5, that is the one withthree cells corresponding to the eigenvalue λ , (cid:39) .
33, in Figure 12.Finally, we comment the last, qualitatively different, buckling. As we note inFigure 10, as before, we have a buckling for the λ value near to 61 .
68 but thistime the bifurcation is linked with two eigenvalues. In fact, for the rectangularplate we have a double eigenvalue λ , = λ , that is the responsible of this doublebifurcation. In practice what we obtained is a point from which start branchingtwo sets of different solutions with one and four cells, respectively in Figure 13and Figure 14.Figure 11: Full order two cells so-lution for the displacement u with λ = 65 (blue branch). Figure 12: Full order three cellssolution for the displacement u with λ = 65 (cyan branch).Figure 13: Full order one cell so-lution for the displacement u with λ = 65 (yellow branch). Figure 14: Full order four cellssolution for the displacement u with λ = 65 (red branch).Same conclusions regarding the convergence error E N and computational sav-ings can be established also in this case. Finally we show that, also for the rect-angular plate, the RB method works well approximating efficiently the solution in23igure 15.Figure 15: Comparison between the full order solution (on the left) and the reducedorder one (on the right) for the displacement u with λ = 65, belonging to the redbranch. Below the reduced basis error plot. In this last numerical result section we want to extend the previous analysis inthe case of two parameters, thus obtaining a 3-D bifurcation plot. Consideringagain the same physical phenomenon, we aim at modelling a compression alongthe shorter sides of the plate, which is no more uniform on that boundaries.Interested in some practical applications of buckling plates in naval engineering,we parametrized the shape of the compression using a new parameter ψ . Infact, the shape of the compression is determined by the function h appearing in(1) thus we can characterize the in-plane compression along ∂ Ω generalizing thecorresponding term λ u xx in (32) obtaining: (cid:40) ∆ u + λ div( σ ∇ u ) = 0 , in Ω u = ∆ u = 0 , on ∂ Ω (35)where σ : Ω → R × , σ (cid:54) = 0 is the plane stress tensor field, which is assumedto satisfy the equilibrium equation: (cid:40) σ t = σ , in Ωdiv σ = 0 , in Ω . (36)24he more general case can be studied, with the uni-axial non-uniform com-pression given by the stress tensor σ ( ψ ) = (cid:34) (cid:16) − ψ y L (cid:17)
00 0 (cid:35) , where ψ ∈ [0 ,
2] is the parameter that takes care of the linearly varying in-planeload. Moreover, we can recover the standard case analyzed before by choosing ψ = 0.Here we want to test the strategy developed in the previous sections in this morecomplex case, where two parameters are involved in the bifurcation phenomenon.Here we restrict ourself to the most physically relevant behaviour, i.e. the evolutionwith respect to ψ of the first buckling, for the generalized system ∆ u + λ div( σ ( ψ ) ∇ u ) = [ φ, u ] , in Ω∆ φ = − [ u , u ] , in Ω u = ∆ u = 0 , in ∂ Ω φ = ∆ φ = 0 , in ∂ Ω . (37)Now we can show some preliminary results on the behaviour of the firstbuckling for the system (37). First of all we can observe in the Figure 16 the3-D bifurcation plot for the square plate, in which we are describing the firstbifurcation point and the post-buckling behaviour, without loss of generality, foreach uniformly sampled ψ ∈ [0 , ψ just introduced. Moreover, we were able to capture correctly the postbuckling behaviour, with results validated by the former analysis with ψ = 0.Finally we want to remark that here the necessity of the Reduced Order Modelsis still more evident. In fact, considering only the full order problem, we have tosolve a huge linear system (as many times as the following nested iteration): foreach Newton step, for each λ in the parameter domain, for each (selected) ψ andfor each initial guess if one is interested on multiple branches.25igure 16: 3-D bifurcation diagramfor the square plate. Figure 17: 2-D projected bifurcationdiagram for the square plate. In this work we have presented a methodology to properly detect a bifurcationphenomena at different levels, with the support of strong and consolidated theoret-ical results and the help of computational reduction strategies, allowing to predictefficiently the buckling. We showed the connection of these physical phenomenawith the eigenvalue problem, which is fundamental in dealing with different ge-ometry or more complex applications. Several numerical tests were performedconfirming the strength of the reduced basis method and its reliability, also innon-linear context. The recovery of eight of the possible solutions shows that,approaching with complex non-linear problems, we need to rely on some backuptool in order to verify that the solution we found is the one we are interested in.Moreover, with this work we showed the consistency of the theoretical results withthe numerical ones, but also the necessity to investigate the reduction strategies toapply this methodology on more complex, real applications. The extension of theresults for the multi parameter application could be also more relevant since theincreasing computational cost.We plan to extend this work in different directions. The first one is towardsthe Brezzi-Rappaz-Raviart theory providing the model with an “a posteriori errorestimate". Furthermore, we want to apply this methodology to other kind ofproblems, such as in fluid structure interaction models, as well as in vibro-acousticsand fluid mechanics frameworks. From the continuum mechanics point of view,we are also interested in the study of different type of plate models, such as theSaint Venant-Kirchhoff, as well as the extension toward the three dimensional26on Kármán equations.
Acknowledgements
The authors thank Dr. F. Ballarin (SISSA) for his great help with the RBniCSsoftware and precious discussion. The authors thank Prof. A. T. Patera for theinspiring conversations and valuable time. This work was supported by Euro-pean Union Funding for Research and Innovation through the European ResearchCouncil (project H2020 ERC CoG 2015 AROMA-CFD grant 681447, P.I. Prof.G. Rozza) by the INDAM-GNCS 2017 project “Advanced numerical methodscombined with computational reduction techniques for parameterised PDEs andapplications", and by the MIT-FVG project ROM2S "Reduced Order Methods atMIT and SISSA".
References [1] E. Allgower and K. Georg.
Introduction to Numerical Continuation Methods .Society for Industrial and Applied Mathematics, 2003.[2] A. Ambrosetti and G. Prodi.
A Primer of Nonlinear Analysis . CambridgeStudies in Advanced Mathematics. Cambridge University Press, 1995.[3] L. Anders, K. Mardal, G. N. Wells, et al.
Automated Solution of DifferentialEquations by the Finite Element Method . Springer, 2012.[4] I. Babuška and J. Osborn. Eigenvalue problems.
Handbook of numericalanalysis , 2:641–787, 1991.[5] M. Barrault, N. C. Nguyen, Y. Maday, and A. T. Patera. An “empiricalinterpolation” method: Application to efficient reduced-basis discretizationof partial differential equations.
C. R. Acad. Sci. Paris, Série I. , 339:667–672,2004.[6] L. Bauer and E. L. Reiss. Nonlinear buckling of rectangular plates.
Journalof the Society for Industrial and Applied Mathematics , 13(3):603–626, 1965.[7] P. Benner, M. Ohlberger, A. Patera, G. Rozza, and K. Urban, editors.
ModelReduction of Parametrized Systems . MS&A Series Vol. 17. Springer Inter-national Publishing, 2017. 278] M. S. Berger. On Von Kármán’s equations and the buckling of a thin elasticplate, I the clamped plate.
Communications on Pure and Applied Mathemat-ics , 20(4):687–719, 1967.[9] F. Brezzi, J. Rappaz, and P. A. Raviart. Finite dimensional approximation ofnonlinear problems.
Numerische Mathematik , 36(1):1–25, 1980.[10] Brezzi, F. Finite element approximations of the von kármán equations.
RAIRO. Anal. numér. , 12(4):303–312, 1978.[11] G. Caloz and J. Rappaz. Numerical analysis for nonlinear and bifurcationproblems.
Handbook of numerical analysis , 5:487–637, 1997.[12] C. Canuto, T. Tonn, and K. Urban. A posteriori error analysis of the reducedbasis method for nonaffine parametrized nonlinear pdes.
SIAM Journal onNumerical Analysis , 47(3):2001–2022, 2009.[13] C. S. Chien and M. S. Chen. Multiple bifurcation in the von Kármán equa-tions.
SIAM J. Sci. Comput. , 18, 1997.[14] P. G. Ciarlet.
The Finite Element Method for Elliptic Problems . Classicsin Applied Mathematics. Society for Industrial and Applied Mathematics,2002.[15] P. G. Ciarlet.
Linear and Nonlinear Functional Analysis with Applications: .Other Titles in Applied Mathematics. Society for Industrial and AppliedMathematics, 2013.[16] P. G. Ciarlet and P. A. Raviart. A mixed finite element method for thebiharmonic equation. In
Proceedings of Symposium on Mathematical Aspectsof Finite Elements in PDE , pages 125–145, 1974.[17] N. Gräbner, V. Mehrmann, S. Quraishi, C. Schröder, and U. von Wagner.Numerical methods for parametric model reduction in the simulation of diskbrake squeal.
ZAMM - Journal of Applied Mathematics and Mechanics /Zeitschrift für Angewandte Mathematik und Mechanik , 96(12):1388–1405,2016.[18] M. A. Grepl, Y. Maday, N. C. Nguyen, and A. T. Patera. Efficient reduced-basis treatment of nonaffine and nonlinear partial differential equations.
ESAIM: Mathematical Modelling and Numerical Analysis , 41(3):575–605,8 2007. 2819] H. Herrero, Y. Maday, and F. Pla. RB (Reduced Basis) for RB (Rayleigh–Bénard).
Computer Methods in Applied Mechanics and Engineering ,261:132–141, 2013.[20] J. S. Hesthaven, G. Rozza, and B. Stamm.
Certified Reduced Basis Methodsfor Parametrized Partial Differential Equations . SpringerBriefs in Mathe-matics. Springer International Publishing, 2015.[21] D. B. P. Huynh, A. T. Patera, and G. Rozza. Reduced basis approximation anda posteriori error estimation for affinely parametrized elliptic coercive partialdifferential equations: Application to transport and continuum mechanics.
Archives of Computational Methods in Engineering , 15:229–275, 2008.[22] D. B. P. Huynh, F. Pichi, and G. Rozza.
Reduced Basis Approximation andA Posteriori Error Estimation: Applications to Elasticity Problems in Sev-eral Parametric Settings , pages 203–247. Springer International Publishing,Cham, 2018.[23] F. Millar and D. Mora. A finite element method for the buckling problemof simply supported kirchhoff plates.
Journal of Computational and AppliedMathematics , 286:68 – 78, 2015.[24] A. K. Noor. On making large nonlinear problems small.
Comp. Meth. Appl.Mech. Engrg. , 34:955–985, 1982.[25] A. K. Noor and J. M. Peters. Reduced basis technique for nonlinear analysisof structures.
AIAA Journal , 18(4):455–462, April 1980.[26] A. K. Noor and J. M. Peters. Multiple-parameter reduced basis techniquefor bifurcation and post-buckling analysis of composite plates.
Int. J. Num.Meth. Engrg. , 19:1783–1803, 1983.[27] A.T. Patera and G. Rozza.
Reduced basis approximation and A posteri-ori error estimation for Parametrized Partial Differential Equation . MITPappalardo Monographs in Mechanical Engineering, Copyright MIT (2007-2010). http://augustine.mit.edu.[28] G. Pitton, A. Quaini, and G. Rozza. Computational reduction strategies forthe detection of steady bifurcations in incompressible fluid-dynamics: Ap-plications to coanda effect in cardiology.
Journal of Computational Physics ,344:534 – 557, 2017. 2929] G. Pitton and G. Rozza. On the application of reduced basis methods tobifurcation problems in incompressible fluid dynamics.
Journal of ScientificComputing , 73(1):157–177, 2017.[30] A. Quarteroni.
Numerical Models for Differential Problems . MS&A SeriesVol. 16. Springer International Publishing, 2017.[31] A. Quarteroni, A. Manzoni, and F. Negri.
Reduced Basis Methods for PartialDifferential Equations: An Introduction . UNITEXT. Springer InternationalPublishing, 2015.[32] RBniCS. http://mathlab.sissa.it/rbnics .[33] L. Reinhart. On the numerical analysis of the von karman equations: Mixedfinite element approximation and continuation techniques.
Numerische Math-ematik , 39(3):371–404, Oct 1982.[34] K. Veroy.
Reduced-Basis Methods Applied to Problems in Elasticity: Analysisand Applications . PhD thesis, Massachusetts Institute of Technology, 2003.[35] K. Veroy and A. T. Patera. Certified real-time solution of the parametrizedsteady incompressible navier stokes equations: rigorous reduced-basis aposteriori error bounds.
International Journal for Numerical Methods inFluids , 47(8-9):773–788, 2005.[36] T. Von Kármán. Festigkeitsprobleme im maschinenbau.
Encyclopädie derMathematischen Wissenschaften , 4, 1910.[37] L. Zanon.
Model Order Reduction for Nonlinear Elasticity: Applications ofthe Reduced Basis Method to Geometrical Nonlinearity and Finite Deforma-tion . PhD thesis, RWTH Aachen University, 2017.[38] S. Zhang and Z. Zhang. Invalidity of decoupling a biharmonic equation totwo poisson equations on non-convex polygons.