Popular in Course
verified elite notetaker
Popular in Mathematics (M)
This 86 page Class Notes was uploaded by Braeden Lind on Thursday October 15, 2015. The Class Notes belongs to MA 797A at North Carolina State University taught by Staff in Fall. Since its upload, it has received 14 views. For similar materials see /class/223716/ma-797a-north-carolina-state-university in Mathematics (M) at North Carolina State University.
Reviews for SP
Report this Material
What is Karma?
Karma is the currency of StudySoup.
You can buy or earn more Karma at anytime and redeem it for class notes, study guides, flashcards, and more!
Date Created: 10/15/15
Verification Techniques for PDE Verification Strategies for PDE Strategies Method of manufactured solutions Compare multiple methods Compare multiple stepsizes Method of nearby problems Compare to analytic analytic solution or special cases Verification Strategies for PDE Example SikoStreifer population model that groups individuals with common traits utx is size density at time t M 8 E 975 im 71 6 9007951 gt 0 gulm 1W owed gulml 0 14096 flm Methods Characteristics analytic Finite differences Finite elements Construct solution Employ biological analysis Verification Strategies for PDE Example Heat model M 82u a 05 m 2 s1n37rt ut 0 ut 2 0 u0 ac sin7r2 One Strategy Consider first amp 9t a82 ut 0 ut 2 0 u0 as sin7r2 that has the analytic solution ut av sin7rm2e mr2t4 Chapter 8 Numerical Techniques The rod beam plate and shell models developed in Chapter 7 generally preclude analytic solution due to the boundary conditions and piecewise nature of material parameters and exogenous inputs is quot t the I t o A 39 tion techniques appropriate not only for simulations and transducer characterization but also for optimization and control design Whereas we focus primarily on the rst objective in this chapter the use of these numerical models for subsequent trans ducer optimization and control design dictates that attention be paid to additional criteria such as adjoint convergence that arise in the context of constrained opti mization 7 the reader is referred to 67 and references therein for details regarding approximation techniques pertaining to optimization and control formulationsi r a of the models we rst employ Galerkin approximations in space to obtain semidiscrete vectorvalued ODE that evolve in time 7 see pages 4177419 of Appendix A for a general discussion regarding the relation between Galerkin and nite element methods This provides a natural setting for linear and nonlinear nitedimensional control design and direct simulationsi Since Galerkin or nite element approximation in space typically yield moderately stiff ODE systems A stable or stiff algorithms are advised when approximating solutions in time this includes trapezoidalbased approximations or routines such as ode155m in MAT LAB In all cases we consider approximation in the context of the weak model formulations since this reduces A 39 t in a natural manner discontinuous material parameters and inputs This necessarily involves the integration of polynomial or trigonometric basis elements which we accomplish using Gaussian quadrature routines chosen to ensure exact integration for linear and cubic basis functions We summarize these numerical integration algorithms in Section Sill Approximation techniques for rods beams plates and shells are subsequently described in Section 82783 umerical approximation of distributed structural models is an extremely broad topic and includes issues such as shear locking and approximation techniques for control design which constitute active research areas Rather than attempt to 373 i 200811 K page 373 374 Chapter 8 Numerical Techniques provide a comprehensive description of numerical techniques we instead summa rize certain fundamental methods appropriate for smart structures and indicate pitfalls and directions required to extend the algorithms to more complex settings The reader is referred to 54479 for nite element theory 390462527 for nite element implementation techniques and 163383406 for the theory of splines variational methods and Galerkin techniques A detailed discussion illustrating nite element implementation in MATLAB is provided in 276 and m les for im plementing various models discussed in this chapter are provided at the website httpwwwsiamorgbooksfr32 Additional references speci c to the various structures will be indicated in relevant sections 81 Quadrature Techniques Consider integrals of the form I f where ab is either nite or in nite and the value I is nite The goal in numerical integration or quadrature is to approximate I by nite sums of the form In Z wifWi i0 where wi and 11 respectively denote quadrature weights and nodes Various ap proximation theories 7 eg ased on Taylor expansions or the theory of orthogo nal polynomials 7 yield choices for wi and I that determine the rate at which In converges to I as n A 00 811 Newton Cotes Formulae To illustrate issues associated with the choice of weights and nodes we consider rst closed and open Newton7Cotes formulae for nite a b The nodes for the two cases are I IO ih where z 0n an 7 a h IO a closed Newton7Cotes formulae 1 7 a h j IO a h open Newton7Cotes formulae n Hence nodes lie in the closed interval a b in the rst case and the open interval a b in the second In all of the following formulae 5 is a point in a b and f is assumed suf ciently smooth so that requisite derivatives exist and are continuous Details regarding the derivation of these relations can be found in 19 The trapezoid and midpoint rules are illustrated in Figure 81 Closed Newton Cotes Formulae 1 n 1 Trapezoidal rule I 3 fltzgtdz 3mm mo 7 grlt5 81 i 200811K page 374 200811 K page 375 e 9 81 Quadrature Techniques 375 a b a b a b Figure 81 a Trapezoidal rule and b midpoint rule for the interval a z 2 n 2 Simpson s rule 12 h ha 4 fltzgtdz f104f11 m2 7 W gtlt5gt 82gt 3 n 3 Simpson s threeeighths rule 7 3h 3h5 fltzgtdz gum sfltz1gt 3m was 7 EM 83gt 4 n 4 Milne s rule I 7 fzdz 7f1032fxl12f1232f137fx4 f65 84 Open Newton Cotes Formulae 1 n 0 Midpoint rule I hS mm 2hfltzogt gm 8 5 2 n 1 b h hs fltzgtdz 3f10 fltz1gt1 BTWlt5 8 6 3 n 2 b 5 fltzgtdz 2fltzogt 7 mm 2fltz2gt1 flt4gtltsgt 8 7 4 n 3 95h5 144 7 5h fltzgtdz m zo mm fr211fzs f45 8A8 200811 K page 376 9 376 Chapter 8 Numerical Techniques Accuracy of the Newton Cotes Formulae A numerical quadrature formula In is said to have degree Ufprecisizm m if I I for all polynomials f such that degf S m and In I for degf m 1 Hence the trapezoid rule has degree of precision 1 Which corroborates the observation that it integrates linear polynomials exactlyi It is observed that formulae With an even index gain an extra degree of precision When compared With odd formulae Which often makes them preferable For later comparison With Gaussian quadrature routines we note that the errors En and degrees of precision DPn for the Newton7 Cotes formulae are E 7 Klnhm n m nodd K2nhn3fltn2gt neven n n odd DPn n1 neven Where KM and K2 are constants and it is assumed that f E Cn2a b if n is even and f E Cn1ab for odd n Whereas increasing n leads to improved accuracy Newton7Cotes formulae are typically restricted to n S 8 to avoid stability issues To achieve the requisite accu racy on large intervals ab including in nite intervals alternatives are required These include Romberg integration techniques Which improve accuracy through Richardson extrapolation composite rules and Gaussian quadrature techniques We summarize next the latter two options Composite Quadrature Techniques An obvious technique to improve accuracy is to partition nite domains a 12 into subintervals and then apply the quadrature rules on each subintervali The manner through Which partitions are constructed depends on the choice of open versus closed Newton7Cotes formulae as illustrated in Figure 82 for the composite trapezoid and midpoint formulaei In both cases is is assumed that f E C2a b and 5 is a point in abi Composite Trapezoida Rule The composite trapezoid rule for n subintervals is Abfzdz g Where h 1 77 and Il a ih for i 0i i i ni It is observed from Figure 82 that if n is doubled present values of are reused thus contributing to the ef ciency of the method fa mm 1 7 b gamems 200811 K page 377 449 81 Quadrature Techniques 377 a b Figure 82 1 Composite trapezoidal rule and b composite midpoint rule with three subintervals Composite Midpoint Rule The composite midpoint formula for even it and g l subintervals is bio 6 W 5 b n2 mm 2hZfz Here h if and z a 2i lh for i 0 n The approximation obtained With n 4 and hence 3 subintervals is illustrated in Figure 82b 812 Gaussian Quadrature Techniques It is observed that for the open and closed Newton7Cotes formulae and correspond ing composite rules the quadrature points z are xed a p7io7i and quadrature weights wi are determined to achieve a speci ed level of accuracy Hence both the accuracy and degree of precision for the methods are roughly equivalent to the de grees of freedom associated With the weights Alternatively one can let both I and wi be free parameters to achieve a maximal order of accuracy This is the basis for Gaussian quadrature routines Which provides the capability for exactly integrating polynomials up to degree 2n 7 1 using n point expansions To provide intuition we initially consider the expansion jl m With n l and n 2 De ning the error as 1 n ampWfWMZ MW 1 i1 we note that for polynomials pm a0 mm amzm Enpm aoEn l alEnz amEnzm It thus follows that the integration rule has degree of precision m if Enzi0 i0m 378 Chapter 8 Numerical Techniques 1 Point Gaussian Quadrature To determine the two parameters 11 and 1017 we consider the constraints E11 07E1I 0 1 1 1dziw107 zdziwlzlOi 1 71 This yields 11 0 and wl 2 and the general quadrature rule 01 1fzdz m 2f0i Note that this is simply the midpoint formula 85 Which is illustrated for the interval a7b in Figure 811bi 2 Point Gaussian Quadrature Here there are four parameters 1171271017102 and four constraints 1 E2zl zldz 7w11114rwgz 0 7 239 07172731 71 This yields the nonlinear system of equations w1w22 wlzl 10212 0 2 2 2 wlzl 10212 g 1011 10213 0 Which has the unique solution 11 7 7 wl 1 12 i 7 112 1 It is noted that the general quadrature rule 1 71 1 819 L l l w w l gt has the same degree of precision as Simpsonls rule 812 Which required three nodes 11 Point Gaussian Quadrature For n gt 27 solVing the nonlinear systems of equations becomes prohibitive and Gaussian quadrature rules are typically formulated using interpolation theory for i 2008111K page 378 200811 K page 379 9 81 Quadrature Techniques 379 orthogonal polynomials By considering families of orthogonal polynomials de ned on the intervals 7117 0700 and 700007 in addition to weighted integrands7 this provides substantial exibility for approximating a broad range of integralsi A complete discussion of this theory is beyond the scope of this chapter and we refer the reader to 19125528 for details regarding Gaussian quadrature routines for both single and multivariate approximation Ga USS Legend re Quadrature GaussiLegendre quadrature formulae are typically de ned in terms of degree n Legendre polynomials n n 13711 Ww1271 n012m on the interval 711 7 see 123465 for a derivation of the Legendre polynomials through application of the Gram7Schmidt process to the sequence 11712 i Note that the rst ve Legendre polynomials are 1301 1 P1 1 1321 612 71 1331 613 7 31 P4z 3514 7 3on 3 The quadrature relation is 1 n mm m Emmi 8 10 1 i1 where the nodes 11 are zeroes of Pn and the weights are wi iluin7 72 7 71 1P IiPn1Ii 7 as summarized in Table Sill For f E 02 1717 the errors are given by f 2 5 EN 5 2n where 5 is a point in 7171 and en m 41 as n A 00 19 As noted previously7 this implies that the quadrature formula 810 is exact for polynomials having degree less than or equal to 2n 7 1 Hence when approximating the solution to weak formulations for structural models7 the degree n is chosen to be commensurate with nite element or spline basis and test functionsi 380 Chapter 8 Numerical Techniques n Nodes zl Weights wi l 0 2 1 2 ii 1 3 0 g 3 i E 3 4 iim A m 60mm iim A m 5187NE Table 81 Nodes and weights for GaussiLegendre quadrature on 7171 Gauss Laguerre and Gauss Hermite Quadrature The integrals in the polarization model 2114 involve the domains 0700 and 700700 One technique for approximating the integrals is to exploit decay exhibited by the integrand to truncate to nite domainsi Alternatively7 one can directly approximate the integrals using orthogonal polynomials de ned on the half line and real line This yields the GaussiLaguerre quadrature relation Amwmdz gutmi and GaussiHermite relation 00 2 n 6 7 m 00 i1 The weights and nodes for these formulae can be found in 528 Gauss Legendre Quadrature on 1117 and Composite Quadrature To evaluate integrals on arbitrary domains using the GaussiLegendre quadra ture relation 8 107 one can employ the linear Change of variables Ab mm Sid711 f d5 811 to map to the interval 7171 It follows that the nodes and weights 11 and wi for a z are related to and 771 for 7171 by the relations ab bia bia TT5i 7 wi 2 771 Z i 200811 K page 380 81 Quadrature Techniques 381 This mapping can also be used to construct composite Gaussian quadrature routines and formulae appropriate for nite element and spline meshes To illus trate7 consider the partition of ab given by zj aj 7 1 forj 0 7n as illustrated in Figure 83 The nodes and weights for the 4point GaussiLegendre rule on the subinterval 11397171139 are then 7 l 7W 7 A 1 7 I34 7 h 2 2 l 7 wl 7 1213M V 1 x1572m 7 49h I2 7 1171 h 5 2 l 7 w 7 1218 812 7 V i V1572m 7 49h I3 I14 7 h 2 7 2 l 7 wg 7 12187xE 7 v 1 WWW 7 i 14 7 111 7 h 2 7 2 l 7 W4 7 mummy We note that this will integrate exactly piecewise polynomials of order up to 7 Figure 83 Partition of ab into n subintervals and position of quadrature points for the 4point GaussiLegendre rule on 11397171139 813 2 D Quadrature Formulae Approximation of weak formulations for plate and shell models requires the numer ical evaluation of double integrals using quadrature rules of the form 12 d n M WWmzz w w 1 C i1 j1 Through the change of variables AbCdfzydydz bf lt3 Jlt agtvWgtdo i 200811 K page 381 382 Chapter 8 Numerical Techniques these integrals can be mapped to the rectangular domain 711 gtlt 711 so we consider rst this case This also provides the framework necessary for numerical integration using rectangular elementsi Finally we summarize formulae for trian gular elements as required for general nite element analysis Rectangular Domains The formulae for rectangular domains are obtained by tensoring 1D relationsi For Gaussian formulae this yields the nodal placement depicted in Figure 814 4 Point Gauss Legendre Quadrature The tensor product of the 1D formula 89 yields the 4point quadrature relation 1 11fltswgtdvds aw u b fan a fab where a b This relation is exact for polynomials of degree 2 Hence this algorithm would be employed when integrating linear quadrilateral elements 9 Point Ga USS Legend re Quadrature The tensor product of the 3point formula from Table 81 yields 1 1 fa was 4 fa7a7 me ma me 71 71 W ab ma me gm where a g 12 0 and c This relation is exact for degree 4 polynomials so it would be used to integrate quadratic elements 390 1 1 1 1 1 1 a b C Figure 84 Quadrature points in a 2D rectangular domain a 1point rule 7 4point rule and a 9point rule Triangular Domains For general nite element analysis it is also necessary to consider quadra ture formulae for triangular domainsi This is often accomplished by considering transformations between physical space I yicoordinates and computational space 200811 K page 382 9 81 Quadrature Techniques 383 0 01 y 110 a 5 b Figure 85 a Master triangular element and 7 global element in physical space X guicoordinates7 as shown in Figure 857 so we summarize rst the construction of local coordinates that are independent of orientation The local coordinates L17L2 and L3 are de ned as the ratio between the perpendicular distance 8 to a side and the altitude h of the side as depicted in Figure 86a This implies that 0 S Li S 1 To elucidate a second property of the elements7 consider the triangle T1 delineated by L1 as shown in Figure 86b7 and the complete triangle Tl From the area relations 312 hi An 3 7 AT 7 it follows that L1 l This motivates the designation of L17L27L3 as area coordinates and establislTes the relation L1 L2 L3 1 From the de nition of the local coordinate L17 it follows that it satis es the property 1 at node i L1 0 at nodes i and k with similar properties for L2 and Lgl When combined with the linearity of the def inition7 this implies that local coordinates also provide the simplex linear elements Ni Nj7Nk depicted in Figure 87 that is7 NiL17 NjL27 1le This proves crucial when de ning quadrature properties for the nite element method a b Figure 86 a Local coordinates L17L27 L3 and 7 triangles T17T27T3 having the areas AT1 AT27 i 200811K page 383 384 Chapter 8 Numerical Techniques Figure 87 Linear elements NiNj and Nk Letting J denote the Jacobian for the transformation between physical and master triangular elements it follows that T fr7ydrdy fL17L27L3lJldL1dL2 T 22 Ji 1 n i Ewi LluLthl i1 The quadrature points weights and degrees of precision for n l 2 3 are summa rized in Table 812 and higherorder formulae can be found in 205462 To illustrate the 1point quadrature relation fewst gal371m T integrates linear functions exactly Whereas quadratic polynomials are integrated exactly by the formula 1 T f57vd5dv g fa20 f0712 flt1212gt1 Note that all of the formulae yield the triangle area A With v ll n Degree of Local Coordinates Weights Geometric Precision L1 L2 L3 w Location 1 1 1 1 1 g g 3 1 a 1 1 1 3 2 5 0 5 g a 1 1 1 5 5 0 g b c 1 1 1 0 5 5 5 C 1 1 1 27 4 3 5 5 5 R a A A u g b 15 15 15 48 u l l E C 15 15 15 48 l u l E d 15 15 15 48 Table 82 Quadrature paints and weights for triangular elements i 2008111K page 384 82 Numerical Approximation of the Rod Model 385 82 Numerical Approximation of the Rod Model In this section we illustrate approximation techniques for the distributed rod models developed in Section 73 For the general model 718 we rst consider a direct Galerkin discretization in space followed by a nite difference discretization in time This yields global mass stiffness and damping matrices a semidiscrete system appropriate for control design and a fully discrete system feasible for simulations Secondly we consider an elemental analysis to demonstrate aspects of nite element assembly often employed for 2D and 3D characterization 821 Global Discretization in Space Consider the weak model formulation 7 18 1 2 I 2 I E u Bu E u dab lo WWW lYAE Amati Eda l W Aa113 a2132lOzjq5 I 813 Bu 821i dz 7 much claw My Ml M where 13 P 7 PR which must hold for 45 in the space of test functions V lt1gt 45790 E 112075 X Rl E H107 7 0 0455 so The goal when constructing Galerkin solutions to 813 is to determine approximate solutions in nite dimensional subspaces VN of Vi 0 construct VN we consider a uniform partition of the interval OJ with points zj jh 0HN and a uniform stepsize h 7 where N denotes the number of subintervals The spatial basis 1 used to construct VN is comprised of linear splines 1 17113971 zj71 zltzj j1g Ij1ir7 IjSISIj1 7 J3917 7N1 0 otherwise 7 814 1 I IN717 INASISIN 45 h 0 as depicted in Figure 88 It is observed that the basis functions satisfy the essential boundary condition 415139 0 0 for j l N Furthermore are differen tiable throughout OJ except at the countable set of gridpoints and hence they otherwise jX NX Xj 1 Xj Xj1 XNii XN Figure 88 Piecewise linear basis functions a 415139 and b i 200811 K page 385 386 Chapter 8 Numerical Techniques are elements in H107Zi Letting goj 45139 Z and 7 q5j7g0j7 the approximating subspace is de ned to be VN span 39i39j The solution to 813 is approximated by the expansion N uNt7 I Zujtq5j 815 Which satis es uN t7 0 0 and can achieve arbitrary displacements at z Zr semidiscrete secondorder matrix system is obtained by considering the approximate solution uNt7z in 813 With the basis functions 1 employed as test functions 7 this is equivalent to projecting the system 813 onto the nite dimensional subspace VNi The interchange of integration and summation yields the system N z N Z N I Aid a Agad 7 YAgVd Elmo0 P W5 1u1t0 c M 1utO 75775 I I I Midi la1ltPiPRgta2ltPePRgt2l asng 0 0 7 1mm mm 3 cmNlttgt Nltzgt mm mm 12 am 3 Which holds for i 17 i i i 7Ni This can be Written as the secondorder vectorvalued system MutQulKut ftA a1PEt 7 PR a2PEt 7 PR2 b 816 Where T 11W lu1t7A A A WNW A The global mass7 stiffness and damping matrices have the components I pA Z jdz 7 ia Norja N 0 MM 2 pA i dem 7 iNandjN 0 I YA dz 7 ia Norjy N Klij 0 817 7 z YA dzkg 7 iNandjN 0 I ch5 q 9dz 7 ia Norij 0 Qtj 7 CA dIC 7 iNandjN 0 i 200811 K page 386 82 Numerical Approximation of the Rod Model 387 Whereas the force vectors are de ned by Z Z HF0 quidz NF0 qbgdzi 818 To evaluate the integrals it is observed that When p and Y are constant the maximal degree occurs in the mass matrix Which is comprised of quadratic polynomials on each subinterval zj1 zjli In this case 2point composite Gaussian quadrature With nodes and weights 111j1glt171 gt wl 121j71glt11 gt 71112 MiD MiD Will provide exact integration on zj1 zj Which yields the tridiagonal matrices 0 0 2710 0 1 271 MpAh K 8 19 71 2 71 0 0 0 0711th andvector b0m1Ti The damping matrix has the representation Q lK When 5 is constant To formulate a rstorder semidiscrete system appropriate for nitedimen sional control design we let 2 u MT and de ne the system matrix A and vectors Ft and B by A7 0 i1 iM lK iM lQ 8 20 Ft 7 0 B 7 0 M lfa M lb The secondorder system 816 can subsequently be reformulated as Azt Ft A a1PEt 7 PR a2PEt 7 PR2 B 8 21 20 zo Where the 2N X 1 vector zo denotes the projection of the initial conditions into the approximating subspacei Temporal Discretization The system 821 must be discretized in time to permit numerical or ex perimental implementation The choice of approximation method is dictated by i 200811 K page 387 388 Chapter 8 Numerical Techniques accuracy and stability requirements storage capabilities and sample rates This can be accomplished using MATLAB routines such as ode155m Which accommodate the moderate stiffness inherent to Galerkin approximation in space Alternatively a trapezoidal method can be advantageous for experimental implementation since it is moderately accurate is Astable and requires minimal storage When imple mented as a single step method For temporal stepsizes At a standard trapezoidal discretization yields the iteration 2m W 5v PW mm al Eum a113Etk1 a2132Etk a2132Etk1 VB zo z0 Where PE 7 PR tk hAt and zk approximates ztk The matrices 71 71 Wlt17Agt lt1 Agt VAtltlliAgt need only be created once When numerically or experimentally implementing the method This yields approximate solutions having 0h2 At2 accuracy For ap plications in Which data at future times tk1 is unavailable the modi ed trapezoidal algorithm Zk1 Wzk Jr VFlttk Jr a1PEtk 7 P13 a2PEtk 7 PR2 VB zo z0 can be employed Whereas this decreases slightly the temporal accuracy for large sample rates With correspondingly small stepsizes At the accuracy is still commen surate With that of the data Remark 821 The approximation of the eigenvalue problem associated with the undamped rod model yields the generalized eigenvalue problem KC w2M4 822 where the sti ness matrix K and mass matrix M are de ned in 819 Hence can be used to approximate the natural frequencies and modes for the undamped rod Remark 822 Due to the presence of both internal damping and damping in the boundary condition at x Z the eigenvalues of the system matrix A de ned in 820 will all have negative real part This property can be used to check the validity of the signs in the boundary condition 713 and weak formulation 813 For example an incorrect formulation which added rather than subtracted the nal boundary contribution will produce eigenvalues ofA having positive real part which is inconsistent with the damping in the model i 200811K page 388 200811 K page 389 9 82 Numerical Approximation of the Rod Model 389 822 Elemental Analysis The spatial discretization technique detailed in Section 821 illustrates the gen eral philosophy of Galerkin approximation including the definition of global basis functions and their use for defining an approximating subspace VN C Vi In this section we summarize a local elemental approach to the problem Whereas the two techniques yield equivalent mass stiffness and damping matrices the latter is sig nificantly more ef cient for general 2D and 3D geometries and hence is employed for general nite element analysis of complex structures To simplify the discussion and facilitate energy analysis we consider the regime employed in Section 732 and take c f P mg cg 1 0 in the weak formulation 813 Additionally we assume that p and Y are constant This model quanti es the dynamics of an undamped and unforced rod that is xed at z 0 and free at z i The space of test functions in this case is V 19307 lt15 E H107 l 0 Ola Local Basis Elements To illustrate the construction of local mass and stiffness matrices we initially consider the approximation of rod dynamics on a local interval 0 h as depicted in Figure 89a In accordance with the assumption that u is differentiable in I at all but a countable number of points we express displacements as utz a0ta1tz 99030 where at a0t a1tT and Lpz 1 le To formulate u in terms of the nodal values ut and urt at z 0 and z Z we note that um 7 1 0 a0t wt 7 1 h a1t ut 39ll at 10 1120 Ti1hl tiutl where a Figure 89 a Rod displacements ut and urt at the left and right endpoints of the local interval 0h 7 Local linear basis functions 451 and 452 200811 K page 390 9 390 Chapter 8 Numerical Techniques By observing that at Sut7 where sT71 1 0 g 1 7 h h it follows that displacements can be represented as utz TIuti 823 N N T Here 4151 17452 contains the local basis elements was 17 g 7 aw depicted in Figure 8i9bi Note that 823 is a local version of the expansion 8 15 employed when constructing approximate solutions using the global basis functions de ned in 814 The local region 07h and global interval zj1zj are analogous to the master and physical triangles depicted in Figure 85 whereas the local basis set 451 452 is the 1D analogue of the 2D elements Ni NV Nk shown in Figure 8 7 7 J7 Hamiltonian Formulation To specify a dynamic model quantifying the displacements u7 we employ the Hamiltonian framework detailed in Section 732 for theNin nite dimensional prob lemi Here we consider displacements u E VN spanq51q 2 as dictated by our approximation framewor i We rst note that for this class of displacements7 the squared strains and velocities can be expressed as u tz uTtSTlDgtzSut um I uTtsTrzsut where Dltzgtw1ltzgtw ltzgtTig 0 FrsoreoTri 1 z The potential and kinetic energies 720 can thus be expressed as h U gumsz DamSue K uTtST 1mm sumo Application of Hamiltonls principle in the manner detailed in Section 732 yields the relation pASTthzdzSutYASTthgtzdzSut 0 824 200811 K page 391 9 82 Numerical Approximation of the Rod Model 391 Local Mass and Stiffness Matrices The weak formulation 8124 can be written as Maud Keu 0 where h h MepAST MIMIs KEYAST Dams 0 0 are the local mass and stiffness matricesi Evaluation of the integrals yields the analytic formulations 1 71 YA 7 K8 T 71 1 J 1 Global Mass and Stiffness Matrices Me pAh OilH le le OilH To motivate the techniques used to construct global mass and stiffness ma trices7 we rst partition the rod interval 07 into two subregions as shown in Figure 8110a 7 hence h 3 With the requirement that u1Tt uggt7 the nodal unknowns in this case are ut um t7 uggt7 1127 tTi Combination of the local relations subsequently yields the global system MuKu0 where the global mass and stiffness matrices are given by 0 1 71 0 MpAh g KYTA 71 2 71 0 0 71 1 We note that the second row in the mass matrix is obtained by summing the element relations AZ 1 1 p7 ltgilli gilQZgt 0 AZ 1 1 p7 lt u21 Eu2rgt 0 after enforcing 111T 1122 A similar strategy is employed when constructing the stiffness matrix and the general procedure is illustrated in Figure 8111ai 41 0 0i 4i i i i i i ii o z 0 Xri gt9 2 a b Figure 810 Partition of the rod domain 0 into a 2 subregions and b n sub regions 392 Chapter 8 Numerical Techniques Figure 811 Combination of elemental mass and sti ness mat7ices Me and K8 to construct global matrices M and K a 2 subregions and b n subregions The process for general partitions zj jh7j 01Hl7N with h is identical and leads to a similar summation process as shown in Figure 8 11bl Enforcing the boundary condition an 0 subsequently yields the global mass and stiffness matrices 2 1 E g 0 0 271 0 0 e 71 271 MpAh 7 19 8 25 1 2 1 a E a 1 2 1 0 0 0 0 71 1 A comparison between 825 and 819 obtained through global analysis illustrates that the matrices are identical when one takes mg k1 0 in the latter set The technique can be extended to incorporate linear and nonlinear inputs by employing the augmented action integral 7 25 and an extended Hamiltonls principlel Internal and boundary damping is incorporated by employing extended constitutive relations as detailed in the previous section Whereas the global discretization techniques and elemental analysis provide the same semidiscrete systems7 the latter technique facilitates implementation for general 2D and 3D geometries discretized using triangular meshesl Additional details regarding nite element techniques for rod models can be found in 36276 823 Examples and Software The performance of the discretized rod model is illustrated in Section 733 in the context of characterizing displacements generated by a stacked PZT actuator from an AFM stage This example illustrates the effects of variable drive levels and the incorporation of frequencydependent hysteresis mechanisms via the nonlinear constitutive relations developed in Chapter 2 he performance of the approximation techniques when used to characterize the hysteretic dynamics of a TerfenolD transducer is reported in 119120 In this case7 the domain wall model was used to provide the constitutive relations which provide the basis for constructing the distributed rod modell n both cases7 the mass7 stiffness and damping matrices have the general representations 8 17 7 1 200811 K page 392 83 Numerical Approximation of the Beam Model 393 M and K have the speci c representations 819 when p and Y are constant 7 and the input vectors are given by 818 The only differences in the models occur in the scalarvalued input relations used to characterize the nonlinear hysteresis MATLAB m les for implementing both the rod model and the model for the stacked PZT actuator employed in the AFM stage can be found at the website httpwwwsiamorgbooksfr32 83 Numerical Approximation of the Beam Model The strategy for approximating the beam models developed in Section 74 is analo gous to that employed for the rod models 7 spline or nite element discretizations in space are used to construct semidiscrete systems appropriate for control design or subsequent temporal discretization for simulations Unlike rod models7 quan ti cation of the bending moments or strain energy associated with bending yields second derivatives in weak beam formulations which must be accommodated by basis functions In Section 831 we summarize the use of cubic Bsplines to con struct approximating subspaces whereas a cubic Hermite basis is constructed in Section 832 through techniques analogous to those of Section 822 831 Cubic Spline Basis We illustrate beam approximation in the context of a cantilever beam7 having a fixedend at z 0 and freeend at z L with surfacemounted patches For test functions 45 in the space V H507 lt15 E H207 W0 W0 0 the weak formulation of the model for linear inputs is Z 2 Z 1 2 2 E w 8w E w d 45 A pwm EMMA w dx Z 63 d2 Z Z d2 17w dz f dzVt kp dz 0 0 0 826 C 6126 where the piecewise constants p7YI CI and kp are de ned in 742 and 744 To formulate approximate solutions based on cubic BSplines7 consider the partition zj jh7 where h and j 07N Forj 71071Nl7 it is shown in 383 that standard cubic Bsplines are de ned by I i 11723 7 I EPW727114 h3 3h2z 7 zj1 3hz 7 zj712 7 3z 7 11397137 z E zj1zj A l W h3 3h2zj1 7 I 3hzj1 7 z2 7 3zj1 7 z3 7 z E zjzj1 1H2 903 7 I E le179512 0 7 otherwise 827 i 200811 K page 393 f f 2008115 page 394 9 394 Chapter 8 Numerical Techniques as illustrated in Figure 812 To accommodate the essential boundary conditions wt7 0 giz 0 and corresponding restrictions 150 0 required of functions in V7 we employ the basis functions 1511 3501 271I af x 7 j1 A 828 jr 7 j27N1 Which are modi ed to ensure that 415139 E V for j 17 i i 7N 1 Finally7 we consider approximate solutions N 1 Wm Z wjlttgt jltzgt in the space VN spanq5j C Vi In a manner similar to that detailed in Section 8217 consideration of wN in 826 With basis functions employed as test functions yields the semidiscrete system MWQWKW fVtb 829 Where wt pm WWW 830 The global mass7 damping and stiffness matrices are de ned by Z Mlij A P i jd1 Z Qtj wta clabg aby dz 8 31 02 KM 0 YIqbQ qbgdz Whereas the force vectors have the components I 1 mi 0 fastdz 7 b 0 we 8 32 Xj 2 Xj 1 Xj Xj1 Xj2 Figure 812 Cubic B splz39ne qu 83 Numerical Approximation of the Beam Model 395 The secondorder vectorvalue system 829 has the same form as 816 which was developed for the rod model and the techniques of Section 821 can be used to construct a rstorder system appropriate for design and control Properties of Cubic Spline Approximates We summarize here properties of the cubic spline approximates 7 additional comparison between this class of approximate solutions and the cubic Hermite ap proximates developed in Section 832 can be found in Section 833 Integral Evaluation From 831 it is observed that the construction of the mass damping and stiffness matrices requires the evaluation of integrals whose integrands are piece wise polynomials of order 2 or 6 when pc an Y are constant To ensure exact integration we employ the composite 4point GaussiLegendre rule 812 which exhibits degree of precision 7 accuracy on each subinterval zj1zji For the xedfree end boundary conditions considered here and constant Youngls modulus this yields the stiffness matrix 648 96 6 712 96 96 754 0 6 6 754 96 754 0 6 712 0 754 96 754 0 6 6 0 754 96 754 0 6 6 0 754 96 754 0 6 6 0 754 84 736 0 6 0 736 48 718 6 0 718 12 For the beam with surfacemounted patches the parameters will be piecewise constant with discontinuities at the gridpoints aligned with patch endsi Partition Construction and Accuracy of the Method For constant parameters p Y1 and 01 the accuracy of the cubic spline approx imate is 0h4i This implies that errors will diminish by a factor of approximately 16 when stepsizes h are have i he same asymptotic accuracy is achieved for piecewise constant parameters associated with the surfacemounted patches if the partition is aligned with the patch edges that is gridpoints zj must correspond with the patch endpoints 11 and 12 The accuracy degradation which occurs if this condition is neglected depends in part on the magnitude of the discontinuity of p Y1 CI and kp at 11 and 12 Alignment of the partition with the patch ends is easily addressed in simu lation and control designs for which patches are bonded at prede ned locationsi i 200811 K page 395 396 Chapter 8 Numerical Techniques A more pertinent issue arises when addressing the problem of optimal patch loca tion for which 11 and 12 are parameters determined through an optimization rou tine 132155202i This necessitates consideration of variable or adaptive meshes and constitutes an active research area Projection Method t is important to note that approximation using the modi ed Bspline basis 8 28 constitutes a projection rather than interpolation method as is the case for the cubic Hermite methods summarized in Section 832 Hence the coef cients 830 do not approximate nodal values of the true solution and modi ed basis functions are required to accommodate essential boundary conditions In this sense Bspline approximation shares a projective kinship with globally de ned spectral methods 7 eg Legendre or Fourier 7 while retaining the sparsity associated with locally de ned polynomial elements As will be detailed in Section 833 the advantage of cubic Bsplines over cubic Hermite elements lies in the fact that half as many coef ficients are required in the former case The noninterpolatory nature of Bsplines constitutes the primary disadvantage which is especially pertinent when accommo dating boundary or interface conditions between components of the structure 7 eg curved and at portions of THUNDER transducersi 832 Cubic Hermite Basis To illustrate the construction of a nite element basis which interpolates displace ments and slopes at the partition points we initially consider the model 826 in the absence of internal or air damping and external forces 7 hence c 7 f V 0 Additionally we take p and Y to be constant to highlight the structure of con stituent matricesi Local Mass and Stiffness Matrices It was detailed in Section 74 that thin beam models quantify both the trans verse displacement and rotation ofthe neutral line so we begin the elemental analysis by quantifying the displacements 144114 and slopes 924 on an arbitrary interval 0 h as depicted in Figure 813 Once we have constructed local mass and stiffness matrices we will extend the analysis to partitions of the full beam interval 0 Z to construct global system matricesi Wr 9 w l x0 xh Figure 813 Displacements wbwr and slopes 924 at the ends of a cubic element on the interval 0 i 200811 K page 396 200811 K page 397 9 83 Numerical Approximation of the Beam Model 397 Because the characterization of Wt wt7t9t7 wr t7 9TtlT involves 4 degrees of freedom 7 we consider cubic representations wt7 I a0 t a1tz a2tz2 a3 tz3 4PT 130 Where at a0t7a1t7a2t7a3tT and Lpz 171712713 By noting that 9t7z t7z and en orcing the interpolation conditions at z 0 and h7 the nodal coef cients can be represented as 833 Wt39ll at Where 1 0 0 l T 1 h h2 h3 0 1 2h 3h2 Substitution of at 39lll lwt into 833 yields the expansion ME I TIWt Where 1 41317 27 437447 comprises the local cubic Hermite basis functions 1 1 was gal e z22z h 7 4531 WW 7 z h 1 1 834 4521 7 h2 7 4541 gram 7 h As shown in Figure 8147 the elements 475139 have displacement or slope values of 0 or 1 at z 07 h This ensures that the coef cients Wt w1t761t7w7 t767 tT interpolate the beam displacements and slopes at z 07hr rom the relations 1 h 82102 l h 82102 UiiYIO lt6I2gt dz 7 KiipO dz M 520 6300 54X h Figure 814 Cubic Hermite basis functions 17 i i i 7414 398 Chapter 8 Numerical Techniques for the potential energy due to bending and kinetic energy7 it follows that for the class of approximate displacements 8337 Y1 h U 7WTtST lDzdz Swt 0 p h K Ev39vTtsT MIMIWt 0 where S 39llquot1 and 0 0 0 0 0 0 0 0 DltzgteomltzgteomTltzgt 0 0 4 m 0 0 121 3612 1 z 12 13 2 3 5 T I z z z FI PI P I 2 3 4 15 13 I4 I5 15 Application of Hamiltonls principle in the manner detailed in Section 732 yields the secondorder vector system Mew Kew 0 where the local mass and stiffness matrices are 156 22h 54 713 12 6h 712 6h M 22h 4h2 13h 73h2 K E 6h 4h2 76h 2h2 8 420 54 13h 156 722 7 e h3 712 76h 12 76h 713 73h2 722 4h2 6h 2h2 76h 4h2 Global Mass and Stiffness Matrices Global mass and stiffness matrices are constructed by combining local relations subject to the constraint that displacements and slopes match at the interfaces To illustrate7 we rst subdivide the beam support OJ into two subregions as depicted in Figure 8 10ai By enforcing the interface conditions WHO WW 7 91W 922t7 i 200811 K page 398 83 Numerical Approximation of the Beam Model 399 we obtain the global matrices 156 225 54 7135 0 0 225 452 135 7352 0 0 ph 54 135 312 0 54 7135 7135 7352 0 852 135 7352 0 0 54 13h 156 722 0 0 713 73h2 722 4h2 Where h 3 Due to the xedend conditions at z 07 it follows that my 1912 01 After reordering the vector of nodal values as WO 10140710243791407924 7 the dynamics are quanti ed by the system MW KW 0 Where 312 54 0 7135 24 712 0 6h M M 54 156 135 7225 K Y 712 12 76h 76h 835 420 0 135 852 7352 7 kg 0 76h 852 252 l 713h 722h 73h2 4h2 6h 76h 2h2 4h2 The process for general partitions zj jh With h is analogous and leads to a similar summation process When constructing global system matrices 7 see Figure 8151 As detailed in previous sections7 internal damping can be incorporated by employing more general constitutive relations based on the assumption that stress is proportional to a linear combination of strain and strain rate Global Discretization We illustrated in Section 812 that either global Galerkin techniques or lo cal elemental analysis could be employed When implementing linear nite element methods The same is true With cubic Hermite elementsi Whereas the local ele mental analysis illustrates the implementation philosophy for general 2D and 3D i 2008111K page 399 400 Chapter 8 Numerical Techniques Figure 815 Combination of elemental mass and sti ness mat7ices Me and K8 to construct global matrices M and K structures7 global Galerkin techniques analogous to those in Section 82 can be more efficient to implement for beam discretizationl The global discretization also demonstrates similarities and differences between the interpolatory cubic Hermite approximates and the projective cubic spline technique discussed in Section 8l3lll For the partition zj jh7 h 7 j 07Hl7N7 the global Hermite basis functions are taken to be 1 I Ijel2l2Ij I hl 7 I E leel7Ijl w7 I g Ij1 I2l2I Ij hl 7 I E Ij7Ij1l 0 7 otherwise I Ijel2I I1 7 I E leel7 Ijl 597 I p Ij1 I2I I1 7 I E Ij7Ij1l 0 otherwise forj 17 l l l 7 Neil The de nitions of 41wa and 41ng are analogous but involves only the interval 1N7171Nll As illustrated in Figure 8167 the global basis functions 45107 and 4597 are the concatenation of the local displacement elements 433 41 and slope elements 417 de ned in 834 and shown in Figure 812 The approximating subspace is taken to be VN spanm m 51 and the approximate solution is N N wNt7 z Z wj ow Z 9139 4597 8 36 j1 j1 We note that by omitting 45100 and 4590 elements 1 E VN are guaranteed to satisfy 110 10 0 which ensures that VN C V Hg 0 The ordering of nodal co ef cients provided by 836 yields banded7 tridiagonal mass7 damping and stiffness matrices that are Toeplitz along all but the last row and column 200811 K page 400 9 83 Numerical Approximation of the Beam Model 401 x xJ XM Figure 816 Global Hermite basis functions 4510 and 4597 The projection of 826 onto VN and use of basis functions as test functions yields the secondorder system MWQWKW fVtb where the 2N X 1 vector Wt has the nodal ordering WW lw1t A 7wNt761t7 A A79NtlTA The global system matrices M Q K and vectors f b have de nitions analogous to 831 and 832 where and 45139 now represent the combined basis 101459 and 41wa 4597 As with the matrices arising from the cubic Bspline discretization discussed in Section 831 the integrals involve piecewise polynomials of degree less than or equal to 6 so exact integration is achieved using the composite 4 point GaussiLegendre rule 812 on each subinterval zj1zjl For constant stiffness Y1 this yields the stiffness matrix 24 712 0 6h 712 24 712 76h 0 6h 712 24 712 76h 0 6h 12 12 76h 76h K l 837 0 76h 8h2 2h2 6h 0 76h 2h2 8h2 2h2 6h 0 76h 2h2 8h2 2h2 6h 76h 2h2 4h2 It is observed that the stiffness matrix K in 835 which was derived through elemental analysis is a special case of 837 when N 2 For the beam with surfacemounted patches the differing values of Y1 in the patch region are simply incorporated in those regions of the partition which coincide with the patch The mass and damping matrices are constructed in an analogous manner We note that the parameter ordering Wt w1t491tl l l wN tt9NtT eliminates the block tridiagonal structure and yields block diagonal matrices where i 2008114K page 401 402 Chapter 8 Numerical Techniques each block has a support of 6 diagonals The disadvantage of this ordering scheme is that the Toeplitz nature of the matrices is destroyed which complicates matrix con struction Additional details regarding elemental and global approximation using cubic Hermite elements can be found in 36203422 833 Comparison between Cubic Spline and Cubic Hermite Approximates The cubic Bspline and cubic Hermite techniques detailed in Sections 831 and 832 illustrate two commonly employed Galerkin techniques for approximating beam models Both provide 0h4 spatial convergence rates as long as partitions are aligned with patches to accommodate discontinuities in mass damping stiffness and input parameters As noted in Section 831 the cubic spline discretization is a projection method whereas the cubic Hermite method is interpolatory in the sense that the coef cients Wt are nodal values of the displacement and slope at the partition points Hence the cubic spline technique is more closely related to general Galerkin expansions 7 eg employing Legendre or Fourier bases 7 whereas the cubic Hermite expansion employs the nite element philosophy which as detailed on pages 4177419 of Appendix A is subsumed in the Galerkin framework The primary advantage of the cubic Hermite method lies in its interpolatory nature This simpli es the enforcement of essential boundary conditions and facil itates characterization of complex structures which require nodal matching at the junction of differing geometries For example the transition from flat tabs to the curved patch region in THUNDER transducers 7 see the models 7103 7106 or 7107 in Section 79 7 is easily accommodated by matching nodal values with Hermite elements whereas it is dif cult to implement with cubic splines The disadvantage of the Hermite approximate is that it requires roughly twice as many coef cients as the spline expansion since both displacements and slopes are discretized The increased dimensionality of system matrices must be accom modated when employing the model for control designs which require realtime implementation 834 Examples and Software Attributes of the discretized beam model when used to characterize the PVDF7 polyimide unimorph depicted in Figure 713 are illustrated in Section 741 Ex perimental validation of the discretized model for a beam with surfacemounted piezoceramic patches is addressed in Chapter 5 of 33 MATLAB m les for implementing the unimorph and beam models are located at the website httpwwwsiam orgbookserZ 84 Numerical Approximation of the Plate Model In this section we summarize approximation techniques for the rectangular and circular plate models developed in Section 75 We consider regimes in which trans verse and longitudinal displacements can be decoupled and focus on approximating i 2008111K page 402 84 Numerical Approximation of the Plate Model 403 the former using Galerkin expansions employing spline or Fourier bases in space As illustrated in Section 85 when discussing thin shell approximation linear elements can be employed to discretize longitudinal displacements if warranted by the appli cation A full discussion regarding nite element methods for plates is beyond the scope of this discussion and the reader is referred to 276390422 527 for details about this topic 841 Rectangular Plate Approximation We consider a rectangular plate with z E 0 and y E 0a as depicted in Fig ure 717 s in Section 751 we assume that the plate has thickness h and has NA surfacemounted piezoceramic patches of thickness h whose edges are parallel with the z and yaxes The plate is assumed to have xededge conditions at z 0 y 0 and free boundary conditions for the remaining two edges To simplify no tation we let 9 0 X 0a denote the plate region Finally we consider linear input regimes with voltages Vh t nglt Extension to nonlinear in put regimes follows immediately when voltage inputs are replaced by the nonlinear polarization relations rom 765 the transverse displacements are quanti ed by the weak model formulation 6 62 62 59 QJW 7 ME 7 2Mmyawy 7 M5972 7 mg do 7 0 838 for 45 in the space of test functions V 1939 45 E 1929 l 45079 x07y 0 for 0 S y S a7 z0 yz 0 0 for 0 S I S E The density p is speci ed by 748 and the moments are de ned in 757 with components speci ed in 7587763 We consider the case of linear patch inputs but note that nonlinear inputs are accommodated in an identical manner The philosophy when approximating the dynamics of 838 is identical to that employed in Sections 82 and 83 for the rod and beam models The relation is projected onto a splinebased nitedimensional subspace VN C V to obtain a semi discrete system appropriate for nitedimensional control design This vectorvalued system can subsequently be simulated by employing nite difference discretizations in time or standard software for moderately stiff systems Consider a partition of Q where rm mhgc yn nhy with has Nil hy Niy and m 0N n 0Ny Using the de nition 828 we de ne modi ed cubic spline basis functions 45m and on the intervals 0 and 0 a The product space basis is then taken to be gtk17 y am 90 y 839 and the approximating subspace is VN span 3221 i 200811 K page 403 63 404 Chapter 8 Numerical Techniques Where Nw N lNyli A 39 t quot 1 t have the Nw WNOJW Zwkt k179 k1 Nz1 Ny1 Z wmnlttgt mltzgt nltygtl m1 n1 The restriction ofthe in nitedimensional model 838 to the nitedimensional subspace VN C V yields the vectorvalued system NA 2YAd31C2 Mw wlef Vitbi 840 Q M7 ltgt lt gt Where Wt w1tuiwNw tTi The mass and stiffness matrices are de ned componentWise by Mljk P Igtjlt1gtkdw Q Kljk Klljk K2ljk Ksljk K4ljk Ksljk Where NA 7 th 2YA03 82 j 82 k lKlle 79 17 V3 171 Xpel y WWdW 7 wt 23mm NA 62 62 K2le 79 1217112 T Tili XPE Ily 812 By2 dw 7 m 23am NA 62ltIgtj 62 W 9 17 W a2 a2 d 7 wt 23mm NA 62 62 K4le 79 1217 112 T 17 vi gxpe 17y By2 812 dw 7 m YACS NA 62ltIgtj a2lt1gtk K5le 79 121 11 1 yA 12de may Braydw The damping matrix Q is constructed in a manner analogous to Ki Finally7 the vectors are de ned by flj mm 9 82 j 82 j bl11v 7 Angry lt W By 63 i 200811 K page 404 200811 K page 405 9 84 Numerical Approximation of the Plate Model 405 The secondorder system 840 has the same form as 816 Which was de veloped for the rod in Section 821 Hence the techniques in that section can be used to construct a corresponding rstorder system appropriate for control design or simulation 842 Circular Plate Approximation The circular plate model 771 With space of test functions 770 represents a set ting in Which the geometry and differential operator are not the tensor product of 1D components When combined With the inherent periodicity in 97 this motivates consideration of basis elements comprised of cubic Bsplines in 7 and Fourier com ponents in 9 We provide here an outline of the approximate system construction and refer the reader to 430 for details regarding this development The circumferential component of the basis is taken to be 45m 9 eime Where m 7M7 7M The form of the radial component is motivated by the analytic Bessel behavior of the undamped plate that is devoid of patches Let 7 denote the nth cubic spline modi ed to satisfy 457 a w 0 7 eg7 see 828 7 along With the condition wise the origin The plate basis is then taken to be 079 Tl llabmr imquot A 0 7 m 0 m 7 1 7 m y 0 As detailed in 430l7 the inclusion of the weighting term Tl ll is motivated by the asymptotic behavior of the Bessel functions as 7 A 0 and ensures the uniqueness of the solution at the origin The approximating subspace is VN span k 0 Which is required to ensure differentiability at Where and approximate solutions have the representation Nw wNt7r76 Zwk k k1 M Nm A I Z wmntTlml T m9 m7M 711 Here N 7 m 0 Nm N 1 7 m f 0 Where N denotes the number of modi ed cubic splines Details regarding the construction of component matrices and vectors are pro vided in 430 The performance of the discretized model When characterizing the dynamics of a circular plate is summarized in Section 75 and detailed in 430 406 Chapter 8 Numerical Techniques 843 Examples and Software The accuracy and limitations of the discretized circular plate model for charac terizing both axisymmetric and nonaxisymmetric plate vibrations are illustrated in Section 752 In the axisymmetric case the model accurately quanti es low to moderate frequency dynamics but overdamps high frequency modes which is characteristic of the Kelvin7Voigt damping model It is illustrated that in the non axisymmetric regime which is truly 2D the model accurately characterizes the dynamics associated with 8 of the 11 measured modes The reader can obtain MATLAB mfiles for the approximation of the rectan gular plate model at the website http www siamorgbooksfr32i 85 Numerical Approximation of the Shell Model The final structure under consideration is the cylindrical shell model 792 devel oped in Section 77 To simplify the discussion we consider fixededge conditions u v w gig 0 at z OJ and hence the space of test functions is v Han x Han x Han where Q OJ gtlt 027r denotes the shell region and Han lt25 6 H16 W079 M6 0 1139 lt15 6 1929 W019 M49 510379 50449 Ol We summarize here the cubic splinePFourier approximation method developed in 130 and illustrated for control design in 131 As noted in the first citation two phenomena which plague the approximation of shell models are shear locking and membrane or shearmembrane locking Shear locking which has also been studied extensively in the context of Reissner Mindlin plate models is due to element in compatibility when enforcing the Kirchhoff Love constraint of vanishing transverse shear strains as the shell thickness h tends to zero 14 Membrane locking occurs when the total deformation energy is bendingdominated and is due to smoothness and asymptotic constraints in the shell model which are not appropriately repre sented by the approximation method 7 eg 2122290380 If these constraints are not satisfied by approximating elements the numerical solution is often overly stiff in the sense that the model exhibits bending dynamics which the approximate solution cannot match As detailed in 290 mesh sizes must be chosen significantly smaller than the shell thickness to ensure accurate approximations with highorder finite elements in such bendingdominated regimesi It is noted that even with such mesh size restrictions loworder finite element methods often fail in such regimesi T e discussion in this section is meant to provide the reader with an overview of issues associated with shell approximation and a brief summary of one discretiza tion techniquei Details and subtleties associated with this topic can found in 39192 and previously cited references i 200811 K page 406 200811 K page 407 9 85 Numerical Approximation of the Shell Model 407 Basis Construction and Approximating Subspaces To i the u 1 1 r Iand4 1 1 i uv and w it is necessary to construct bases for nitedimensional subspaces of H3 and Hg This is accomplished using linear and cubic splines modi ed to accommodate the xed boundary conditions We consider a uniform partition along the zaxis With gridpoints In nh h an n iHNi or n Hl 7 l we employ the linear splines 1 1717717 zn71Szltzn am 7 was 7 E w 7 z as g z s W 0 otherwise and modi ed cubic splines gm 7 2111 7 21I n 71 wnz 7 3541 n72mN72 mm 7 2N1z 7 2N1z n 7 N 71 Where the standard B splines are de ned in 827 For n liHN 7 l and m 7M l l l M the product space bases in complex form are taken to be Put 97 I im9 un I ltIgtvk971 im9 v 95 wk 9 I eimeqbwn To provide an equivalent real form one can employ the representation cosmt9 buk61 l unz mluiM sinmt9 With similar de nitions for vk and wki The approximating subspaces are VuN span uk 21 VUN Span 1gtvkl 1 Vliv span 309221 andthe A 39 quot are bythe 1 Nu Max 7 Zuklttgtlt1gtuklt67zgt 161 N44 vNtt9z kat vk 91 841 61 Nw wNtt9z Zwkt wk 91 161 f f 2008115 page 408 44 408 Chapter 8 Numerical Techniques For a single partition it follows that Nu N1 Nw 2N 7 l2M 1 Through the construction of 45451 and 451 the approximate displacements uN v and wN satisfy the fixedend conditions and VN VuN gtlt VN gtlt Viv C Vi We note that v consideration of real components yields the equivalent representation N71 uNtt9z Zum wmz M N71 M N71 21 2 umnt cosm9 uz 71 Z Ema sinmt9q5un z with similar expressions for UN t 9 z and wN t 9 System Matrices To determine the generalized Fourier coef cients uktvk t and wkt we employ the same approach as in previous sections and orthogonalize the residual with respect to the linearly independent test functions used to construct the approx imating subspacesi To construct the resulting vectorvalued system we consolidate the coef cients in the vectors mm mm mm 1W E 7 VO E 7 WO uNu t UN t WNw t The full set of coef cients is then represented as 19t utvtwtT where N Nu Nu N The mass stiffness and damping matrices have the form Um M Vm 7 Wm U11 U12 V11 W12 W11 K U21 U22 V21 W22 W21 7 U31 v31 ELWSk 7 611 612 V711 V712 W11 Q U21 U22 V21 W22 W21 U31 v31 ELWSk 7 and the exogenous vectors are ftfu7fv7fwT bbubv0Ti 200811rK page 409 9 85 Numerical Approximation of the Shell Model 409 The various submatrices contain individual components Which arise When the weak formulation 792 is restricted to VNr For example7 the approximation of the mass and stiffness components in the longitudinal equation of 792 yields Umljr 7 thlt1gtttltIgtutdwt Q YhR B i uk B i ui Urrljr7917V2 61 71 w th 64gt alt1gtu7 WllljkiQlilQ 66 BI dwt th alt1gtu7 willjkA Wk at d Yh 64gt alt1gtu7 Mir9m art 66 d wt 9 6 Unljk Q 2R1V 89 89 d iftij 7 Pruittth Q Rh B i ui lbuljQl71 BI dw With similar expressions for the remaining submatricesr In the usual manner7 the secondorder system MM QW KW 7 W la1PEt 7 P12 a2PEt 7 P292 b 190 190 50 191 can be reformulated as the rstorder Cauchy equation 21 7 Am W la1PEt 7 P12 a2PEt 7 P292 B Z0 lo 7 Where z 19t i9tT and 7i ell 7 3 l 7 7i i 7M K 7M Q M t M b The system in this form is anemable to simulation7 parameter estimation and control design Note that the system can be adapted to alternative boundary conditions through modi cations of the rst and last basis functions Flexibility in this regard is also a hallmark of the method Numerical Techniques for PDE A good many times I have been present at gatherings of people who by the standards of traditional culture are thought highly educated and who have with considerable gusto been expressing their incredulity at the illiteracy of scientists Once or twice I have been provoked and have asked the company how many of them could describe the Second Law of Thermodynamics The response was cold it was also negative Yet I was asking something which is about the scientific equivalent of Have you read a work of Shakespeare s C P Snow The TWO Cultures and the Scienti c Revolution Canonical Forms and Important PDE Hyperbolic xo xo Ct FirstOrder ut out 0 mm for go 0 at Solution utx fm ct t SecondOrder utt czum 0 0737 2 a ut01 9amp5 1 1 ct Solution uta 5 mg ct fact g 7 t gydy I Canonical Forms and Important PDE Parabolic eg Heat equation M0510 f ut 0 ut L 0 Solution 00 2 ut I Z aneiaknt sinn 1121 where An m and T L an 2 fa sinnmdm 0 Canonical Forms and Important PDE 124 I 0 Elliptic eg Laplace s equation 74 741015 2 hm ux yy0 x y69 u0 y0 v2 0 7 u 2 Way aw7 may 6 39 Solution 1amp0 0 I ua y Z Bn sinAnc sinhny n1 where An 77quot and B 2 af msin 3dm n asinhnb 0 2 n Physical Examples Steady state heat behavior steady state structural deformations Canonical Forms and Important PDE Rod Model 8211 62M 8311 Ma w 39 C 8x2at f Boundary Conditions fgt ut 0 0 0 6 6211 Nt kgut game mgt Weak Formulation 3 6 4 an 6 d e 2 6 a f dm mom Jigglam mg72ut W 0 8t at for all lt15 80 lt25 155 6 V Canonical Forms and Important PDE Beam Model 9271 871 FM Pwl39W E amg fta95 wt0 Ziggy 0 Mt i Afw 0 93w 2 where Mt 3 YIaca m kpxVt wd 392 Weak Formulation e 2 Z Z 2 2 9 w 8w 8 wd 02 dwz l 3 2 18 Z 2 8 u d lt25 d Sb 1 0 c amzat dmz da 0 f dO katdm2 dm for all 6 V H 0 ApproximationTechniques for the Rod Model Galerkin Linear Basis j it em Ij1 a mj mfmj I jil acjilgacltxj 0 otherwise lijil iL j jl Approximate Solution N j1 System N E N Z N g jzzliijtO PA i jd23 jtO cA 2 9dac ujt0 YA Q 9dm e 2 0 f id MUNGWNW Ce NUWNW ml Nt No 1M fori1N ApproximationTechniques for the Rod Model Galerkin 2ndOrder Vector System M1105 Qu Kut ft where ut we uNtT and l pA mdx i Norja N 0 M1ij z pA i deme iNandjN 0 Z YA dx ia Norjy N Kb OZ YA 9dke iNandjN 0 l CA 3dac ia Norj7 N Cl 1239139 e CA 9dICg iNandjN 0 Hi 2 OZ f id ApproximationTechniques for the Rod Model Galerkin Matrices Constant coefficients 0 0 2 1 0 0 1 2 1 a E a 1 2 1 MpAh KYh A 1 2 1 0 0 g m74 0 041th FirstOrder System it Azt Ft z0zo where 0 ll 0 A Ft M 1K M 1Q M lf t Note Codes available at httpwww4ncsuedursmithSmartMaterialSystemsChapter8 ApproximationTechniques for the Rod Model Finite Element Motivating Problem pAg Zf YAg 0 Local Basis Elements Take utm a0t a1tx 90Tat where at a0t a1tT and m 1 mlT Then 11405 1 0 a0t Ugt Ur t urt 1 h a1t a 0 35 h ut 39H at By observing that at Sut where 7 1 0 Sz39l1 11 1 h it follows that displacements can be represented as ut 2 Twut where Tc 1 g g ApproximationTechniques for the Rod Model Finite Element Action Integral t1 A KeUdt to 1 K Kinetic Energy t1 h Z pAut7miYAu t Hth U Potential Energy to 0 Here ut x uTtSTlDxSut ufh ac uTtSTlFxSut where Dlxgtaltwzltxgtlg F90ltPTI Thus Aw t 1 ilTtMei1t uTtKeutdt where u u 1391 and h h Me pAST 0 MIMIS KeYAST 0 MIMIS Me Local Mass Matrix K5 Local Stiffness Matrix gt M16 2 pAh mild wlr x wlgt mlH l 51 m 1K in ApproximationTechniques for the Rod Model Finite Element Hamilton s Principle See Chapters 7 8 and Appendix of Smith 2005 d 60 where T 1777 and it is assumed that 17050 71t1 0 Here 751 AW 5T i 11 5hT M501 677 u 517TK5 u 571 dt t11 d gt d EAW l ET 5 hTMe 1151391139157397TMeh t0 nTKe u an u EMT Ken dt t1 hTMEil nTlKeu dt since uTMe TMeil to 151 n Men K511 dt integration by parts to Local System Mew Ken 0 ApproximationTechniques for the Rod Model Finite Element Global Matrices Consider two subregions so h Here ut lumtU2etU2rtlT and M Ku0 where 1 1 3 6 0 1 1 0 szAh g K 1 2 1 l 1 0 g E 0 1 1 Note Second row obtained by summing 1 1 YA 39 pAh gum gum T ulg u 0 P 1 i L u 1 1 21 0 1 I O l I PAh g zz l 671 T 22 uzr 0 139 39 after enforcing uh um and 711T llgg Finite Difference Techniques for the Rod Model Consider first utt K2 Uq51 F ut 0 ut 0 14096 f a Mam 990 Grid A itj l L i I ih j k System in 2 Kg 2 1 E Earl 2 m 231M h2 W44 271211 2414 Fan Ii gt ij1 m2 Ema Hal 21 m2 ij 13H HEW Ii Initial Conditions 31 371 N gt mil 2 um 2kgxk 2 gt 131 m7fi1 fi71 1 m2f139i 9995239 ngUoim Finite Difference Techniques for the Rod Model Now consider utt Rum K2um F ut0 ut 0 u0 96 M 5 MO I W Question What is finite difference scheme Galerkin Method for Beam Model Weak Formulation See Section 83 of Smith 2005 e 2 E E 2 2 8 w 910 8 wd ab 13 3 2 e e 8wd d I d d k t d 0 C 8m25tdm2 m 0f 0 pVdm2 m for all gas 6 V H302 Basis where I Ij723 90 E j72axj71 h3 j71 j712 Ij13 I E Ij1j A 1 19 m hg 3h2j1 95 3hj1 IV 3j1 903 90 E Ijvijrn Ij2 903 7 I E 1jlamj2 0 otherwise Galerkin Method for Beam Model Approximate Solution N1 wNOZw Z wjt j j1 System MWQWKW fVtb where 2 MM 20 P i jd z QM 0 WWW 01 39 01 2 mi Azfgbidq M Aiddag Lecture 1 Introduction Essentially all models are wrong but some are useful George EP Box Industrial Statistician Lecture 1 Introduction What is validation verification and uncertainty quantification Why do we care Motivation Etrema Products Inc SolidDrive httpwwwetremausacom TerfenolD Transducer Schematic of TerfenolD transducer in the SolidDrive TerfenolD rod can be modeled as a rod with elastic and clamped boundary conditions For uniform input fields spring equation with nonlinear inputs provides good approximation of rod dynamics What is nonlinear and hysteretic behavior of fHu Magnetomechanical Material Behavior AM AM A 400 MPa 5 Material depends on both magnetic elds and stresses Consider application of periodic eld followed by application of positive and negative stresses Data from Pitman 1990 validation veri cation and uncertainty quanti cation Model Development Model 1 lv 7 H266Nm H80Nm H132Nm DA 04 DA 4 quot2 02 02 ill Ll Experimental Data 1 n 1 n 1 n 1 n 1 0 1 a Model 2 Question Which model is bette Model Development Original data Craik and Wood 1970 a 3 V V 39 A l m r s 39 l Representation of data used for Model 1 Jiles 1995 Representation of data used for Al Model 2 Smith and Dapino 2006 H266Nm HEDNm H132Nm 7 04 04 04 quot7 quot KN 02 W 02 I 02 I Ir J l l 4 c c 71 D 1 D it D 1 D 71 D l 0 Questions g r 1 Are scales accurate quot 2 ls data accurate 3 What is source of loss in upper loop Model Development Model 1 l Fl H266Nm HEDNm H132Nm DA 04 04 y H 02 02 02 4quot th l Experimental Data 7 it D 1 D it D 1 D it I 1 0 Model 2 Notes 1 Numerical implementation of Model 1 is inaccurate verification 2 Neither model is true representation of physical device validation Modeling Issues Model Qualification Reality 4 A Model Computer Validation Simulation V Computer Model Analysis V Conceptual Model A Programming Model Verification Verification Process Conceptual Model Correct Answer Computational Model Analytic solutions Highly resolved numerical Computational SOlutlonS SOIUtlon Verification Test 39Benthark solutions Verification The process of determining that a model implementation accurately represents the developer s conceptual description of the model and the solution to the model Note Verification deals with mathematics Validation Process Real World Conceptual Model Correct Answer Provided by Experimental Data Computational Model Computational Solution Validation The process of determining the degree to which a model is an accurate representation of the real world from the perspective of the intended model users Benchmark cases System analysis Statistical analysis Validation Process Note Validation deals with physics and statistics Reference Oberkampf Trucano and Hirsch 2004 Experiment Model y Q Viewgraph Norm Experiment f Model Response Input Numerical Error Validation Metrics I Experiment O Model Response Input Deterministic I Experiment O Model Response Input Nondeterministic Computation Experiment o Model Response Input Experimental Uncertainty Validation Metrics Example Metric 1 H266Nm HEDNm H132Nm 04 04 at 04 2 7 1 I quot2 w 02 z 02 I 39I l 1 c I r gquot 10 10740 An 10 Data H268Nm HEuAm H132Nm 04 an on 02 02 02 1 u 1 o 4 a mu Model Metric 2 A B tesla Validation Metrics Example uw 0 Stress MPa Difficulties with Using Existing Experimental Data for Validation 1 Incomplete documentation of experimental conditions Boundary and initial conditions Material properties Geometry of the system 2 Limited documentation of measurement errors Lack of experimental uncertainty estimates Note It is critical that modelers numerical analysts and experimentalists work together Reference Oberkampf Trucano and Hirsch 2004 Background of Validation and Verification Vaidation and verification began to be formalized during the ASCI Accelerated Strategic Computing Initiative Program to replace physical nuclear testing by computer experiments A large body of validation and verification research has been performed at Sandia National Lab Verification Techniques for IVP The man of science has learned to believe in justification not by faith but by verification Thomas Huxley Many a good newspaper story has been ruined by oververification James Gordon Bennett Conceptual Model Computational Model Correct Answer Analytic solutions Highly resolved numerical Computational L SOlUtlons SOIUtlon Verification Test 39Benthark solutions Verification The process of determining that a model implementation accurately represents the developer s conceptual description of the model and the solution to the model Verification Versus Data Fit Recall A fit to data may not verify the accuracy of the numerical solution Model 1 H266Nm HEDNm H132Nm 04 04 04 02 02 02 m m r g Experimental Data 1 n H 1 n in 1 n 10 Model 2 IVP Verification Techniques True or False Decreasing the stepsize always increases accuracy Example 3391 2y7 110 3 3 l True 25 Appr0XImate 2 C 9 L 5 1 5 g HO LIJ 1 05 00 05 i 15 2 1 395 10 3915 Time S 10 1O 10 Number of Intervals N 5 Note 6 016 02k IVP Verification Techniques Analytic Solutions Strategy Test code using analytic solution Example Solution of mu 01 ku 2 F0 coswt is F0 ut clerlt 0262 mzwg w22 02w2 coswt 6 where n and T2 are solutions of the characteristic equation tug km 6 satisfies cos6 mw3 w2A and A 7712003 0122 czwz Notes 39 This provides an initial test of codes Analytic solutions may not exist for nonlinear problems or systems IVP Verification Techniques Linear Systems Example Finite element discretization of rod model from Lecture 2 yields the semidiscrete system 205 A205 Ft z0zo where 0 I 0 A 71 71 Ft 71 M K M Q M t 2 1 5 a 0 0 2 1 0 0 1 2 1 a 5 a 1 2 1 MpAh K 1 2 1 E 5 a 1 2 1 0 0 WTZ 0 0 1 1th Details See Chapter 8 of Smith and later lectures IVP Verification Techniques Linear Systems Example 2 2 AZ F 20 2 20 This has the solution t zt eAtzo eAtTSFsds o t PeJtPleOJr PeJtTSPT1Fsd5 0 Notes Feasibility of the method depends on the accuracy and difficulty of determining P and J MATLAB eigenvalue routines may work for relatively large matrices gt 1000 This provides analytic solution for all time IVP Verification Techniques Example Consider the HIV model from Adams et al 2005 T1 A1 lel 1 ek1VT1 T2 A2 clng 1 fek2VT2 Tf 1 ek1VT1 6Tf mlETf T 1 fek2VT2 6T ngT V NT6T1 712 CV 1 E101k1T1 1 fEP2k2T2lV bET1 11 E dET1 11 EA E TfT Kb TfT Kd E 6EE Problem How does one guarantee correct numerical solutions Notes Metatheorem To model accurately the propagation of error it is necessary to do at least two integrations RD Skeel Thirteen ways to estimate global error Numer Math 48 120 1986 PE Zadunaisky On the estimation of errors propagated in numerical integration of ordinary differential equations Numer Math 27 2739 1976 IVP Verification Techniques Asymptotic Analysis Numerical Solutions T 106 1 2 1o3 5 10 101 104 10391 0 50 100 150 0 50 100 150 106 T1 104 T2 102 VW 101 M 10392 102 0 50 100 150 50 100 150 V 102 E 106 103 MM N 10 101 0 50 100 150 50 100 150 Strategy Consider values of y where ftyO IVP Verification Techniques Different Methods Note This usually does not quantify error T 10x105 1 ode45 ode113 8 Note RelErr 1e5 67 7 4 00 5b 160 150 Time 5 Warning Remember Note Use of same code for rhs will 23 yield same error for different methods Mt M3 1 l u 0 IVP Verification Techniques Different Stepsizes Techniques Use differing orders to quantify error Use differing stepsizes as basis for extrapolation techniques T 5 1 10x 10 Re Err 1e5 Re Err 1e8 a Note T1 T1 5562 6 4 27 L39 0 i i 0 50 100 150 Time s IVP Verification Techniques Extrapolation Techniques Euler s Method 0k implies yk W 1119 0092 2J2 W 1219 0092 Order p Method yk Mb 1116 00439 My W 11011017 006 Elimination of 1th yields 3116 qu r b k y yk qp1 0 IVP Verification Techniques Nearby or Pseudoproblem Reference Zadunaisky 1976 Use to test implementation of RHS boundary conditions etc Idea fit splines to numerical solution to generate analytic solution to nearby problem 6 T1 10 i Approximate Solution Spine Fit 1057 4 10 0 50 100 150 Time s IVP Verification Techniques Nearby or Pseudoproblem Original Problem M variables W ftyt We yo Error ln absence of roundof f error en 2 yn ytn Pseudo or Nearby Problem 239 t 20 Dot 2050 13050 Note Method is also applicable to BVP and some PDE where Dt Pt ftPt Error Behavior ay where wt 205 yt 1W t Ptwt D00 0 wt2
Are you sure you want to buy this material for
You're already Subscribed!
Looks like you've already subscribed to StudySoup, you won't need to purchase another subscription to get this material. To access this material simply click 'View Full Document'