轉(zhuǎn)向節(jié)外文文獻(xiàn)翻譯--一種無摩擦接觸問題的有限元方法

1、<p>　　一種無摩擦接觸問題的有限元方法 </p><p><b>　　摘要</b></p><p>　　文章提出了一種新的解決包括整體可能經(jīng)歷一定的運(yùn)動和變形的接觸性問題的有限元方法。這個(gè)方法是基于兩結(jié)構(gòu)接觸性問題變成兩個(gè)同時(shí)發(fā)生的問題的分解，最終結(jié)果是在幾何學(xué)上而不是在離散的接觸表面上。一個(gè)表面接觸單元是特定設(shè)計(jì)在無條件下，允許在兩接觸表面間進(jìn)行一直的

2、牽引動力的傳輸。</p><p>　　關(guān)鍵詞：無摩擦性接觸，大的變形，有限元</p><p><b>　　1.引言</b></p><p>　　有限元方法廣泛應(yīng)用在解決接觸性問題上。從完全的計(jì)算角度來看，接觸性的檢測和隨后的強(qiáng)制性約束的實(shí)現(xiàn)是有待解決一般算法結(jié)構(gòu)的發(fā)展的兩個(gè)重要問題。從康里和瑟瑞吉，禪和圖巴的早期作品以后，很多方法論都在接觸機(jī)械

3、學(xué)的文獻(xiàn)中被提出了。一個(gè)相當(dāng)廣泛的在題目中的調(diào)查在參考書目中發(fā)現(xiàn)。</p><p>　　目前的工作是參與適合的解決在大的運(yùn)動和變形的兩機(jī)構(gòu)的接觸性問題有限元分析方法的發(fā)展。這一類問題在許多實(shí)際應(yīng)用中有特別重要的意義，如金屬成形過程和車輛事故分析。各種商業(yè)和科研計(jì)算機(jī)編碼采用的算法，以解決這種問題。拉格朗日算子方法和它們的規(guī)格化（罰數(shù)和增廣的拉格朗日方法）是典型地被用在執(zhí)行不可測知中。</p><

4、;p>　　工作性條件的集成方法的選擇在線性動力平衡方面的微弱形式下與接觸牽引動力有聯(lián)系并在接觸原理建設(shè)中起著關(guān)鍵的作用。節(jié)點(diǎn)正交的使用包括一個(gè)（二個(gè)）點(diǎn)或一個(gè)（兩個(gè)）基座或點(diǎn)和面的接觸。假設(shè)存在一些連續(xù)的接觸表面，另外的整合法則也都是合適的。</p><p>　　根據(jù)兩機(jī)構(gòu)問題接近于連續(xù)的兩同時(shí)發(fā)生問題，目前的文章的重要貢獻(xiàn)是鑒定一般的程序。像傳統(tǒng)的雙行程運(yùn)算法則，兩相互作用機(jī)構(gòu)的表面被用來分析而不需要采用

5、中間的接觸表面(任意的被選擇)。提出的方法的重要優(yōu)勢與雙程節(jié)點(diǎn)上表面的運(yùn)算法則相比，如果允許一個(gè)一體化規(guī)則的簡單解釋用于接觸面，并為可容許領(lǐng)域的適當(dāng)選擇，允許一個(gè)物體到另一個(gè)物體的連續(xù)壓力的精確傳播。根據(jù)源于鐐銬的工作中的斑貼試驗(yàn)，及其以后的概括，恒壓（在大小和方向）的精確代表性的能力被看做是一個(gè)健全的必要條件和總體接觸算法的收斂。</p><p>　　接觸力學(xué)的一個(gè)簡潔的闡述在第2節(jié)中出現(xiàn)，特別強(qiáng)調(diào)在接連發(fā)生的

6、算法的發(fā)展中常常用公式表示。在第3節(jié)中提出并分析了一個(gè)二維接觸單元，而在第4節(jié)中，提出和討論了正在使用的這個(gè)單元的選擇數(shù)值模擬的結(jié)果。在第5節(jié)中，給出了結(jié)束的評論。</p><p>　　2. 兩物體接觸問題</p><p>　　考慮物體認(rèn)為等同于線性空間的開連通集，配有標(biāo)準(zhǔn)基(,,)和歐幾里德準(zhǔn)則。至少有一個(gè)物體是假定可變的。在參考位形中的一個(gè)典型的重要點(diǎn)是在線性構(gòu)造中重要點(diǎn)用向量數(shù)學(xué)地說

7、明，對于每一個(gè)t給出，移植向量是通過來規(guī)定的。</p><p>　　假定至少在內(nèi)，映射在定義域內(nèi)是連續(xù)和可轉(zhuǎn)置的。向量，在線性結(jié)構(gòu)（以t變化）中它的范圍是用和來確定的，因此，則有和 。并且的外在的單位標(biāo)準(zhǔn)記作。</p><p>　　多個(gè)物體（包括一個(gè)物體）的任何系統(tǒng)的提議都受制于物質(zhì)的不可測知性，在引文11（224頁）被Truesdell和Toypin規(guī)定。這意味著這兩機(jī)構(gòu)問題一直有<

8、;/p><p><b>　?。?）</b></p><p>　　在任何給定的時(shí)間上，所說的這兩個(gè)機(jī)構(gòu)都與它們的邊界子集C有聯(lián)系，當(dāng)且僅當(dāng)</p><p><b>　　(2)</b></p><p>　　根據(jù)上面的定義，每個(gè)物體的定義可以被唯一的分為三個(gè)相互排斥的區(qū)域，根據(jù)下面的式子</p>

9、<p>　　當(dāng)?shù)依锟巳R和諾艾曼邊界條件分別地被強(qiáng)加在 上 。雖然沒有明白的之處，但是一般依照時(shí)間來說，它應(yīng)該不包含, 和 C。</p><p>　　不連續(xù)函數(shù)，可能是多值的，每個(gè)物體的邊界按照如下的定義：對每一個(gè)，是給定的為</p><p><b>　　（2a）</b></p><p>　　當(dāng)是在這樣的式子中看做指數(shù)是1。的凸面性使

10、的值是惟一的，盡管有如此的一個(gè)幾何條件的限制也不能強(qiáng)加在開端。的一個(gè)完整的相似的定義生成下面的結(jié)果：</p><p><b>　　(2b)</b></p><p>　　又有對每一,滿足。定義方程式(2a) 和(2b)即指不連續(xù)函數(shù) 和 在C上是相等的，都等于0.也就是</p><p><b>　　(3)</b></p

11、><p>　　因此不可測知條件（1）根據(jù)上面的不連續(xù)函數(shù)被重寫作</p><p><b>　　。</b></p><p>　　在不存在慣性效應(yīng)時(shí)，方程式的局部形式控制著如下給出的每個(gè)物體的運(yùn)動，</p><p>　　div in (4a)</p><p&

12、gt;　　on （4b）</p><p>　　on (4c)</p><p>　　on (4d)</p><p>　　當(dāng)是克西的應(yīng)力張量，是質(zhì)量密度，是每單位質(zhì)量的質(zhì)量力，是法定的界壁位移，是在上的法定牽引向量。</p>&l

13、t;p>　　標(biāo)準(zhǔn)加權(quán)殘值法的應(yīng)用，與拉格朗日乘數(shù)的引入p ≥ 0對費(fèi)解的限制相結(jié)合，結(jié)果在運(yùn)動方程的若形式中指出，位移解方程（ 4 ）和拉格朗日乘數(shù)外地p滿足</p><p><b>　　(5a)</b></p><p><b>　　(5b)</b></p><p>　　和q是對所有的允許函數(shù)。在大部分的缺失不被用在

14、方程(5a)中時(shí)，則不連續(xù)函數(shù)定義在C上。位移區(qū)域?qū)儆诳臻g即</p><p>　　質(zhì)量函數(shù)也屬于空間，定義為</p><p>　　允許函數(shù)p ≥ 0是分段連續(xù)的。</p><p>　　為了證明方程式(5a)的正確性，把這項(xiàng)工作看做是允許函數(shù)和上沿C的觸點(diǎn)壓力，即為</p><p><b>　　(6)</b></p&

15、gt;<p>　　在C上，在沒有摩擦和回顧時(shí)，柯西引理應(yīng)力矢量意味著</p><p><b>　　(7)</b></p><p>　　借助于方程(7)，方程(6)可寫為</p><p>　　這表明拉格朗日乘數(shù)區(qū)域是自然等同于正常的在接觸范圍內(nèi)的牽動引力（壓力）。壓力場p一般假定僅僅是分段光滑的，因此認(rèn)為每個(gè)物體的特征材料的接觸面在

16、C的附近。而且由于</p><p>　　on (8)</p><p>　　不等式(5b)是從(4d)的形式中得到的，而且假定q是非負(fù)的。</p><p>　　為了進(jìn)一步證明拉格朗日乘數(shù)在兩物體接觸問題中的作用，為接著發(fā)生的近似數(shù)值提供了一些動機(jī)，用(2b) 和(8)改寫(5b)中在接觸面C上構(gòu)成整體所必需的，即</p><p><

17、;b>　　(9)</b></p><p>　　符號僅僅被用來強(qiáng)調(diào)接觸面C不普遍的看做的一部分，就像(3)中所表明的。整式(9)表明，在每個(gè)物體的范圍內(nèi)，適當(dāng)?shù)慕o定一個(gè)拉格朗日區(qū)域，(5b)可以在每個(gè)接觸面上分別作用。然而，很顯然區(qū)域應(yīng)該滿足在（一般地）接觸面上的線動量的平衡，即為</p><p>　　以上資料數(shù)據(jù)將被應(yīng)用在(5a) 和 (5b)的近似解決方案中。.<

18、/p><p>　　(9)式的推導(dǎo)用下面的程序，等式(5a) 和(5b)被記作</p><p><b>　　(10a)</b></p><p><b>　　(10b)</b></p><p>　　當(dāng)，對所用可容許的和都適用。</p><p>　　(5a) 和 (5b)的典型的罰數(shù)調(diào)

19、整是在的條件下獲得。 </p><p>　　當(dāng)<·> 是一個(gè)Macauley支架和時(shí)，上面的（5a）式的積分關(guān)系是一個(gè)連續(xù)的損失參數(shù)。約束條件(4d)是無約束的（因此允許受約束的滲透來替代）替代區(qū)域是確定的，則有</p><p><b>　　(11)</b></p><p>　　在約束條件下，原始的無約束問題等同與一個(gè)凸面

20、的最小化問題，序列的解決方案當(dāng)時(shí)表明收斂是約束問題的解決方案(5).實(shí)踐中，罰數(shù)方法成功的用在平衡狀態(tài)下與最小總勢能不相關(guān)的情況下。</p><p>　　由于罰數(shù)參數(shù)的有限值的普及率，接觸面C的一個(gè)獨(dú)特的定義和缺陷函數(shù)并不是容易地應(yīng)用在(11)中的。在提出的明確表達(dá)中關(guān)鍵的一步是(10a) 和 (10b),易于接受的規(guī)則化的罰數(shù)的引入，例如</p><p>　　由(12)，等式(10a)可

21、以寫為</p><p><b>　　(13)</b></p><p>　　上面的(13)中的任意兩個(gè)積分在形式上是相同的，即積分邊界來源于一個(gè)新問題的規(guī)則化罰數(shù)。盡管在目前的背景下這兩種積分在缺陷函數(shù)的兩個(gè)定義(2a) 和(2b)的基礎(chǔ)上是有聯(lián)系的，但是這種分解要設(shè)計(jì)要涉及到重要的計(jì)算，這在下一節(jié)中再討論。</p><p>　　這里提出的事態(tài)發(fā)

22、展可以很容易地靠歸納兩體問題擴(kuò)展到多體接觸問題。</p><p>　　3. 三維接觸的兩個(gè)應(yīng)用</p><p>　　這篇文章的其余部分專門討論在明確表達(dá)罰數(shù)的基礎(chǔ)上的離散二維接觸問題和兩表面上的缺失函數(shù)的鑒定，像在第二節(jié)中所述。這里待確定的兩個(gè)關(guān)鍵問題是接觸單元的幾何作圖和有限元近似值的可允許區(qū)域的選擇。</p><p>　　3.1 接觸面的離散化</p>

23、;<p>　　如果有一個(gè)接觸面的離散化是連續(xù)的，那么接觸面 和 是分別在 和的條件下是唯一確定的。撇開推行細(xì)節(jié)，接觸面是被一系列正常的從 到的投影（不定靠近點(diǎn)）所離散，如圖表2中所示。的離散化采用相同的步驟。因此每個(gè)接觸單元在接觸面上的立體空間元素與其對立面相聯(lián)系。</p><p>　　顯然地，非光滑邊界的離散會導(dǎo)致單位法線和缺失函數(shù)的不連續(xù)。這里目的不是設(shè)法回避這個(gè)問題，盡管它的結(jié)果不可忽略，特別

24、是在滾動問題中。對于一個(gè)特殊的離散來說會導(dǎo)出代表兩個(gè)層面的光滑邊界，見參考資料13 。</p><p><b>　　3.2 有限元區(qū)域</b></p><p>　　在這節(jié)中假設(shè)一個(gè)具體的接觸有限元是建立在方程式(10a) 的(12)的基礎(chǔ)上。選擇有限維領(lǐng)域的指導(dǎo)方針是以前的作品特別是在青年受理領(lǐng)域位移和壓力導(dǎo)致單方面接觸單元，能夠穩(wěn)定地復(fù)制費(fèi)解的條件（即滿足的LBB條

25、件）是準(zhǔn)確的。這種收斂性分析，目前盡管并不適用于這里的兩體接觸問題的動非線性這方面，但是它為可允許區(qū)域的選擇提供了一個(gè)準(zhǔn)則。</p><p>　　數(shù)值積分典型地把所有邊界術(shù)語用在接觸面上。不考慮離散化使用空間的元素（例如三角形的三個(gè)定點(diǎn)和四邊形的四個(gè)點(diǎn)）的線性規(guī)則。在單個(gè)接觸單元中的缺失函數(shù)是關(guān)于軌跡邊界坐標(biāo)系的非標(biāo)準(zhǔn)函數(shù)。因此，一般情況下在上的數(shù)值積分是不準(zhǔn)確的，所有的引入錯(cuò)誤的積分規(guī)則直接影響到用公式的確切表

26、示。所提出的接觸單元是建立在從發(fā)源到標(biāo)準(zhǔn)Q9元素（等參的九個(gè)節(jié)點(diǎn)）的移置區(qū)域和引用的辛普森標(biāo)準(zhǔn)上的。</p><p>　　A NOVEL FINITE ELEMENT FORMULATION FOR FRICTIONLESS CONTACT PROBLEMS</p><p><b>　　SUMMARY</b></p><p>　　This ar

27、ticle advocates a new methodology for the finite element solution of contact problems involving bodies that may undergo finite motions and deformations. The analysis is based on a decomposition of the two-body contact pr

28、oblem into two simultaneous sub-problems, and results naturally in geometrically unbiased discretization of the contacting surfaces. A proposed two-dimensional contact element is specifically designed to unconditionally

29、allow for exact transmission of constant normal tract</p><p>　　KEY WORDS frictionless contact; large deformations; finite elements</p><p>　　1. INTRODUCTION</p><p>　　Finite element

30、 methods are used extensively in the solution of contact problems. From a purely computational standpoint, detection of contact and subsequent satisfaction of the impenetrability constraint are the two key issues to be a

31、ddressed in the development of a general algorithmic framework. Numerous methodologies have been proposed in the literature of compulational contact mechanics since the early works of Conry and Seircg,and Chan and Tuba.A

32、 fairly comprehensive survey on the topic is f</p><p>　　The present work is concerned with the development of finite element methods suitable for the solution of two-body contact problems in the presence of

33、large motions and deformations. This class of problems is of particular significance in numerous practical applications, such as metal forming processes and vehicular crash analyses. Various commercial and research compu

34、ter codes employ algorithms for the solution of such problems. Lagrange multiplier method,and their regularizations (penalty and </p><p>　　The choice of integration method for the work-like term associated w

35、ith the contact tractions in the weak form of linear momentum-balance plays a pivotal role in the construction of contact elements. Use of nodal quadrature involving the contacting nodes of one of (resp. both) surfaces y

36、ields the standard one-pass (resp. two-pass) node-on-surface algorithms.Other integration rules are also applicable, provided there exists a continuous discretization of the contact interface.</p><p>　　The m

37、ain contribution of the present paper is in the identification of a genneral procedure according to which the two-body problem is approached as a sequence of two simultaneous sub-problems. As in the traditional two-pass

38、algorithms, the surfaces of both interacting bodies are used in the analysis without need for introduction of an (often arbitrarily chosen) intermediate contact surface. The main advantage of the proposed approach over t

39、he two-pass node-on-surface algorithms if that it allow</p><p>　　A brief exposition to contact mechanics is presented in Section 2, with particular emphasis on formulations to be used in the ensuing algorith

40、mic developments. A two-dimensional contact element is proposed and analysed in Section 3,while the results of selected numerical simulations using this element are presented and discussed in Section 4.Concluding remarks

41、 are given in Section 5.</p><p>　　2. THE TWO-BODY CONTACT PROBLEM</p><p>　　Consider bodiesidentified with open and connected sets in linear space , equipped with canonical basis (,,) and the usu

42、al Euclidean norm. At least one of the bodies is assumed to be deformable. A typical material point of in the reference configuration is algebraically specified by vector of material point in the current configuration

43、is given byat each time t, and the displacement vector is defined according to</p><p>　　The mapping is assumed smooth throughout its domain and invertible at least on .Body and its boundary in the current co

44、nfiguration (at time t) are identified with and ,respectively, henceand .Also, the outer unit normal to is denoted by .</p><p>　　The motion of any system of bodies (including a single body) is subject to t

45、he principle of impenetrability of matter, as stated by Truesdell and Toypin in Reference11 (p.244).This implies for the two-body problem that at all times</p><p><b>　　（1）</b></p><p>

46、;　　At any given time, the two bodies are said to be in contact along a subset C of their boundaries if, and only if,</p><p><b>　　(2)</b></p><p>　　It follows from the above definetion

47、 that the boundary of each body can be uniquely decomposed into three mutually exclusive regions according to</p><p>　　Where Dirichlet and Neumann boundary conditions are enforced on and , respectively. Alt

48、hough not explicitly noted, it should be clear from the above that , and C generally depend on time.</p><p>　　Gap functions , possibly multi-valued, can be defined on the boundary of each body as follows: f

49、or each , is given by </p><p><b>　?。?a）</b></p><p>　　Where is such that see Fifure 1.Convexity of renders single-valued, although such a restrictive geometric condition will not

50、be imposed at the outset . A completely analogous definition for yields</p><p><b>　　(2b)</b></p><p>　　Where,again, for each , satisfies .Defining equations (2a) and (2b) imply that

51、 gap functions and are identically equal to zero on C, namely that</p><p><b>　　(3)</b></p><p>　　Consequently, impenetrability condition (1) can be rewritten in terms of the above g

52、ap functions as</p><p>　　At the absence of inertial effects, the local form of the equations governing the motion of each body αare given by</p><p>　　div in (4a)</p&

53、gt;<p>　　on (4b)</p><p>　　on (4c)</p><p>　　on (4d)</p><p>　　Where is the Cauchy stress tensor, the mass den

54、sity, the body force per unit mass,the prescribed boundary displacement, and the prescribed traction vector on .</p><p>　　Application of the standard weighted-residual method,in conjunction with the introdu

55、ction of Lagrange multiplier p≥0 for the impenetrability constraint,results in the weak form of the equations of motion, which states that the displacement solution of equations (4) and the Lagrange multiplier field p s

56、atisfy</p><p><b>　　(5a)</b></p><p><b>　　(5b)</b></p><p>　　For all admissible function and q. Without loss of generality is used in equation (5a) for the de

57、finition of gap function on C. Displacement fields belong to spaces with</p><p>　　and weight functions belong to spaces defined as</p><p>　　The admissible functions q≥0 are piecewise continuou

58、s.</p><p>　　To prove that equation (5a) holds, note that the work done by contact forces along C on admissible functions and is given by</p><p><b>　　(6)</b></p><p>　　A

59、t the absence of friction and recalling that on C, Cauchy’s lemma on the stress vector implies that</p><p><b>　　(7)</b></p><p>　　With the aid of equation (7), equation (6) is writte

60、n as</p><p><b>　　,</b></p><p>　　Which shows that the Lagrange multiplier field is naturally identified with normal traction (pressure) along the contact region. The pressure field p

61、 is generally assumed to be only piecewise smooth, thus allowing each body to feature material interfaces in the neighbourhood of C. Moreover, since</p><p>　　on (8)</p><p>　　inequalities (5b

62、) follow from (4d) and the assumed non-negativeness of q.</p><p>　　In order to further clarify the role of the Lagrange multipliers in the two-body contact problem and provide some motivation for the ensuing

63、 numerical approximations, use (2b) and (8) to rewrite the integral in (5b) on surface C as</p><p><b>　　(9)</b></p><p>　　The notation is employed to merely emphasize that the contac

64、t surface C can be selectively viewed as part of ,as indicated in (3).The integral expression (9) suggests that, given a properly defined Lagrange multiplier field on the boundary of each body, integration of (5b) can b

65、e performed separately on each of the contacting surfaces. However, it is clear that fields should satisfy balance of linear momentum on the (common) contact surface, namely that </p><p>　　The above observa

66、tion will be exploited in the approximate solution of (5a) and (5b). </p><p>　　Folowing the procedure used in the derivation of (9), equations (5a) and (5b) are rewritten as</p><p><b>　　(1

67、0a)</b></p><p><b>　　(10b)</b></p><p>　　for all admissible and ,where .</p><p>　　The classical penalty regularization of (5a) and (5b) is obtained by setting</p&

68、gt;<p>　　in the last integral term of (5a),where <·> is the Macauley bracket and is a constant penalty parameter. Constraint conditions (4d) are relaxed (thus permitting controlled penetration to take p

69、lace) and displacement fields are determined so that</p><p><b>　　(11)</b></p><p>　　for all .Under restrictions such as formal equivalence of the original unconstrained problem to a

70、 convex minimzation problem，the sequence of solutions as can be shown to converge to the solution of the constrained problem (5).In practice, penalty methods are used successfully even when the equilibrium state is not

71、 associated with minimization of the total potential energy.</p><p>　　Due to the occurrence of penetration for finite values of the penalty parameter ,a unique definition of the contact surface C and gap fun

72、ction g is not readily available for use in (11).A crucial step in the proposed formulation is the introduction of a well-defined penalty regularization for (10a) and (10b),such that</p><p><b>　　(12)&l

73、t;/b></p><p>　　With the aid of (12),equation (10a) becomes</p><p><b>　　(13)</b></p><p>　　Each of the last two integrals of (13) identical in form to the boundary inte

74、gral emanating from penalty regularization of a Signorini problem. Although in the present context these two integrals are clearly coupled by the definitions (2a) and (2b) of gap functions , this decomposition has import

75、ant computational implications, as will be discussed in the next section.</p><p>　　The developments presented here can be easily extended to encompass the muli-body contact problem by reducing it to a series

76、 of coupled two-body problems.</p><p>　　3. APPLICATION TO TWO-DIMENSIONAL CONTACT</p><p>　　The remainder of this article is devoted to the discretization of two-dimensional contact problems base

77、d upon the penalty formulation and the identification of gap functions on both surfaces, as suggested in Section 2.The two key issues to be addressed here are the geometric construction of contact elements and the choice

78、 of admissible fields for the finite element approximation.</p><p>　　3.1. Discretization of the contact surfaces</p><p>　　A continuous discretization of the contact surfaces is advocated. Surfac

79、es and are uniquely determined as those on which and, respectively. Setting aside implementational details, contact surface is discretized by a series of normal projections (not necessarily closest-point) from to ,a

80、s shown in Figure 2. An analogous procedure is followed for the discretization of .Consequently, each contact element relates a single spatial element edge on surface to the opposite surface.</p><p>　　The

81、apparent non-smoothness of the discrete boundaries results in discontinuity of unit normals and gap functions. No attempt is made here towards circumventing this problem, although its effect might not be negligible, espe

82、cially in problems of rolling contact. For a special discretization that results in smooth boundary representation in two-dimensions, see Reference 13.</p><p>　　3.2. Finite element fields</p><p>

83、;　　A specific contact finite element is suggested in this section, based upon equeations (10a) and (12). The choice of finite dimensional fields is guided by previous works especially on the Signorini admissible displace

84、ment and pressure fields lead to unilateral contact elements that are able to stably replicate the impenetrability condition (i.e. they satisfy the underlying LBB condition) and are accurate. Such convergence analysis, a

85、lthough not currently available for the kinematically non-linea</p><p>　　Numerical integration is typically employed for all boundary terms on the contact surface. Discounting discretizations that use straig

86、ht-edge spatial elements(e.g. three-node triangles and four-node quadrilaterals),the gap functions within a single contact element are non-poly-nomial with respect to local boundary co-ordinate systems. Thus, numerical