INTERNATIONAL JOURNALFOR NUMERICAL METHODS IN ENGINEERING, VOL. 39, 805-828 (1996) IDENTIFYING GLOBAL/LOCAL INTERFACE BOUNDARIES USING AN OBJECTIVE SEARCH METHOD S. SRINIVASAN Department of Bioengineering, 301 Rhodes Bldg., Clemson University, Clemson, SC 29634, U.S.A. S. B. BIGGERS JR. Department of Mechanical Engineering, Clemson University, Clemson, SC 29634, U.S.A. R. A. LATOUR JR. Department of Bioengineering, Material Science and Engineering Program, Clemson University, Clemson, SC 29634, U.S.A. SUMMARY One of the key components in computational mechanics procedures using global/local methods is the specification of the global/local interface. Global/local interfaces are usually specified by visually examining some measure of response such as colour-coded contour plots of stresses or strain energy. However, when both global and local domains are modelled in three dimensions, such a specification is not as obvious, and it is presented to specify the global/ lacks objectivity and uniqueness. An Objective Search Method (OSM) local interface in three dimensions in a precise, repeatable and automated manner. The OSM performs the search incrementally in all directions in three dimensions radiating from a location of interest until certain generalized guidelines are satisfied and the global/local interface is identified. The OSM is suited to problems where localized phenomena exist but where their domains are not known a priori. The generalized guidelines for the OSM require the identification of nodes lying on the external surfaces of the model. As an important component of the OSM, a unique method to identify surface nodes has been developed and is also presented. Finally, the uniqueness, sensitivity and versatility of the OSM is illustrated using two example problems and the computational effort involved with the OSM is discussed in the context of a third example problem. KEY WORDS: global/local analysis; 3-D finite elements; global/local interface identification INTRODUCTION The need for an Objective Search Method (OSM), as described in this paper, presented itself in ongoing research on the structural performance of a composite femoral component for use in total hip joint replacement.' The complex 3-D geometry of the implant and surrounding bone, the material inhomogeneities, and the anisotropic properties of the composite material all combine to present a difficult challenge to the structural analyst. These complexities and similar ones in many other applications, such as machine and vehicle design with advanced materials, require both accuracy and efficiency from the analysis methods. Global/local (G/L) readily lend themselves to satisfying these requirements. The central theme of G/L appro ache^"^ is the use of a coarse model sufficient for accurate representation of the global structural response. Following global structural analysis, detailed CCC 0029-5981/96/050805-24 0 1996 by John Wiley & Sons, Ltd. Received I July 1994 Revised 3 January 1995 806 S. SRINIVASAN, S.B. BIGGERS JR. AND R. A. LATOUR JR. local models of specific areas of interest within the global model are analysed to obtain detailed stress states in the localized areas of concern. In some of these approaches, the critical areas of the structure may need to be accounted for in the global structure at least in an approximate way. In some structural problems, such as in the analysis of composite structures for hip prostheses applications, the critical areas are not known a priori. In such instances, G/L approaches based ~ . ~ the global and local finite on the specified boundary displacement (SBD) m e t h ~ d , where element models are uncoupled, find the best applicability. Furthermore, the uncoupled SBD method is best suited for implementation within most commercially available general purpose finite element codes. In the SBD method, the following key steps are performed sequentially. First, an adequate global analysis is performed which captures the overall response of the structure. Second, critical spots on the structure which require detailed analysis are identified from the global solution. Third, an appropriate G/L interface is specified to provide a separation, into independent domains, of the global and local behaviour. Fourth, a process for interpolation of global-independent fields onto the local domain along the specified G/L interface is formulated. Finally, an adequate local analysis is performed to capture localized phenomena of interest. One of the critical steps in the SBD method is the specification of the G/L interface. The interface should lie outside regions with strong local gradients and be located along areas of gradually changing or stable stress fields (or any other field representing structural response) across the G/L b ~ u n d a r y The . ~ G/L interface is usually defined by the user, through visual observations of colour-coded contour plots of stresses, strains, strain energy density, or some other measure of response, at locations where the fields are changing gradually. This approach of specifying the G/L interface is fairly acceptable and straight forward in 2-D problems. However, in cases where a 3-D model is required to represent adequately the structural response (e.g. composite hip prostheses or other complex thick composite structures), the interface between G/L domains is not readily apparent and is very difficult to specify. Also, the subjectivity involved with visual identification will cause such a definition of the interface to be conservative (usual scenario) and non-unique. This problem of interface specification is further complicated in situations involving 3-D structures with multiple interacting stress concentrations. Thus, a more objective approach is needed in defining the G/L interface. The purpose of this paper is to introduce an Objective Search Method (OSM) which automatically searches for and establishes a unique, reproducible G/L interface in 3-D space. The OSM performs the search incrementally in all directions in 3-D radiating from a defined ‘hot spot’ until certain generalized guidelines are satisfied and the G/L interface is located in an inclusive and complete fashion. The term ‘objective’ as used in this paper implies that the process is unbiased, based on facts and analyst independent (opposed to being subjective and analyst dependent). Traditionally, in G/L problems, the G/L interface is located by ‘eye-balling’ colour-coded contour plots of structural response and, this process is, as such, subjective and analyst dependent. In the proposed approach, the G/L interface is identified in an automated, consistent, repeatable, analyst independent manner based on logical, rational sets of criteria and the search process is thus an objective search method. Depending on the problem being analysed, the OSM could be sensitive to the selected measure of structural response, and selection of the measure of structural response is the only area for which analyst judgement is required. The overall process involved with the specification of the G/L interface using the OSM is similar to error estimation but with significant differences. For instance, the Zienkiewicz-Zhu (Z-Z) method” is a simple error estimate that provides a representation of the error in the finite element discretization in predicting structural behaviour. The Z-Z method falls within a class of techniques which are often referred to as flux-projection methods’ where projection-type GLOBAL/LOCAL INTERFACE BOUNDARIES 807 procedures are employed to extract post-processed approximations of the flux which are then used to estimate measures of error in the finite element solution. The Z-Z error estimate provides a basis for controlling the error in the finite element solution by adaptively refiningl0+l2areas with error larger than allowed (specified by the analyst) and fusing others where error is negligible. The OSM, on the other hand, provides the means to isolate regions and define independent domains where a transition between global and local behaviour occurs and is best suited for use with the more efficient G/Lmethods. The interface defining this transition between global and local behaviour is identified by the OSM based on the rate of change of smoothed fluxes with respect to distance away from a localized event such as a concentration induced by geometric/material discontinuities.The Z-Z error estimate provides a basis for deciding whether an individual global element needs further subdivision/fusing, whereas the OSM provides a consolidated account of all elements that need to be considered together in isolation and whose independent refinement will improve the solution. In fact, if the Z-Z error estimate is used as the measure of structural response with the OSM, the OSM then provides an automated decision process to the analyst by monitoring the rate of change of error (more appropriate than merely thresholding the error), thereby defining an interface between global and local behaviour based on rate of change of error and providing the barest essential zone or region which, if refined independently, will enable the minimization of discretization errors. The Z-Z error estimate provides discrete information about the error as energy norms or a number of others while the OSM provides the means to consolidate that information and to present the same in an isolated form. Therefore, the OSM is not a replacement for the Z-Z estimate but could be easily coupled with the Z-Z estimate for use with the more efficient G/L methods where details in structural response are treated independently but in the context of global structural behaviour. The OSM is suited to problems where localized phenomena exist and need to be analyzed using G/L approaches but where the locations of the phenomena are not known a priori. This approach is especially needed to address problems where both the global and local models are fully three dimensional and where traditional eyeballing methods are too conservative(therefore contributing to the G/L computational expense) and inconsistent. METHODS The first step in performing a G/Lanalysis is, of course, to set up an adequate global model and to solve for structural response. The critical spot(s) in/on the structure where the selected measure of structural response, S, reaches an extreme value and where detailed local analysis needs to be performed are identified as 'hot spots'. Once the hot spot is identified, the OSM is used to uniquely locate the G/Linterface according to a set of defined OSM parameters. The OSM seeks and objectively identifies a bounding surface around the hot spot, or multiple interacting hot spots, where the G/L interface can be satisfactorily specified. In general, the OSM is performed about each hot spot by beginning the search at the hot node and then sequentially stepping from node to node of the global model in a designated search direction while sampling the structural response of interest at that node. A response gradient is then calculated between each set of traversed adjacent nodes and compared to a user-defined structural response gradient tolerance (SGT) limit. The search proceeds incrementally from node to node until the response gradient reaches the SGT limit, at which point, the last node in the designated search direction is identified as a node on the G/L interface. This search method is performed in all directions in space radiating away from the hot spot and a unique G/Linterface is thus identified which completely encloses the hot spot node. Once this G/L interface is identified, global independent fields (e.g. displacements, velocity, temperature) are interpolated to the local domain and the boundary 808 S. SRINIVASAN, S. B. BIGGERS JR. AND R. A. LATOUR JR. conditions for the local model are established. The accuracy of the local model solutions depend on these boundary conditions. Therefore, the location of the G/L interface should be at points in the structure where the response gradients have reached an acceptably low value (defined by the user and referred to in this paper as the SGT limit). Some of the possible situations that can arise in 3-D structures which must be dealt with in an OSM are that the hot spots can either be on the surface of the structure or lie internally within the structure. In order to cover both of these possibilities, the search procedure must have the capability to traverse along the surface of the structure or proceed into the interior of the structure seeking to identify the G/L interface. These possibilities and the required capabilities of the OSM are rationalized below to provide a set of generalized guidelines for the 3-D radiating search process. These guidelines, outlined in Table I and illustrated in Figure 1, allow the identification of the complete G/L boundary in 3-D surrounding the hot spot according to a user-defined SGT limit. The procedure outlined in Table I is described in detail in the following sections. Identifying surface nodes The generalized guidelines (Table I) require that each node the search may encounter from one increment to the next be identified as a surface node or as a node that lies internally within the structure. If a node or point on the structure is completely surrounded by material, then it must be an internal node in the structure, otherwise it must be a surface node. The search method must therefore have the capability to determine whether or not a given node is completely surrounded by material. This capability is implemented in a unique way and is explained in the following text. Table I. Summary of the generalized guidelines for the OSM that searches incrementally in all directions about a designated ‘hot spot’ Hot spot is surface node (Figure l(a)) (1) Search stays on the surface in the designated direction until response gradients reach the user defined SGT limit or (2) Search goes into the structure’s interior and (i) Proceeds until response gradients reach the user defined SGT limit; or (ii) Reaches another surface of the structure a Hot spot on Surface Hot spot is internal node (Figure l(b)) (1) Search proceeds in designated direction until response gradients reach the user defined SGT limit or (2) Search proceeds until a surface node is reached b. Internal Hot Spot Figure 1. Possibilities the OSM is required to handle. Numbers in figure refer to their corresponding numbers in Table I 809 GLOBAL/LOCAL INTERFACE BOUNDARIES For every node the search may encounter (i.e. node 0 in Figure 2), a determination of the number of adjacent elements N is made. Each adjacent 3-D element (element i ) must contribute three nodes that are adjacent to node 0, i.e. nodes 1,2 and 3 for adjacent element i (Figure 2). A spherical triangle of polar radius 5 = 1 is constructed using the angles subtended at 0 by nodes 1,2 (angle a), 2,3 (angle b) and 3 , l (angle c). The surface area of the oblique spherical trianglet3 formed from adjacent element i is given by d2E A. =' 180 where - b) tan (S - C) tan (S 2 2 S= (a + b + c) 2 Then the sum of the surfaceareas AT contributed by all N adjacent elements which share corner node 0 is given by N If the node of interest is an interior node, AT will be equivalent to the surface area of a sphere of radius 6 = 1, where Asphere = 4 d 2 Therefore, a simple comparison between AT and Aspheredetermines whether the given node is in the interior or the surface by the following procedure: AT <1 If-{ Aspher. = 1 then node 0 is a surface node then node 0 is a internal node In this way, a determination of whether or not any node 0 is a surface node is made. This procedure is accurate up to 16 digits if the coding is done in double precision. By this method, %/"- 1.2.3 Adjacent nodes to ocxk 0 contributedby element i Figure 2. Identifying a surface node by constructing spherical triangles of polar radius 6 using nodes 1,2,3 and 0 contributed by element i 810 S. SRINIVASAN, S.B. BIGGERS JR. AND R. A. LATOUR JR. Figure 3. Finite element models of a cubic structure and accompanying sample calculations a set of nodes occurring on the surface of the finite element model can be constructed independent of and prior to the ‘search’, or the nodes can be checked for surface location at each increment and direction in the search procedure. As an example of this procedure to determine a given node’s location (i.e. interior or surface node), consider a finite element model of a cubic structure (Figure 3). The angles at all nodes subtended by any two adjacent nodes are 90”due to the regular geometry and mesh. Also shown are the calculations involved for nodes A, B and C and the determination of location where node A is identified as an interior node, whereas nodes B and C are determined to be located on the surface of the structure. Searching in all directions around hot spot The goal of the OSM is to determine a unique bounding surface around the hot spot where the local effects have died out and where the interface between global and local domains may be specified. In order to be all inclusive and to define completely a bounding surface in three dimensions, the procedure must search in all directions in three dimensions radiating away from the hot spot. Due to the discretization involved in the finite element models, the search is performed incrementally where ‘incrementally’ means that the OSM searches to locate the boundary node by proceeding from one node to the next adjacent node. The bounding node where the specification of the G/L interface may be made is the node identified after the nth increment in any given direction where one of the generalized guidelines (Table I) is satisfied. The search procedure starts at the node determined to be critical according to the chosen measure of structural response. The first step in the search procedure is to determine the number of directions that need to be searched in order to identify the entire G/L boundary. Searching in all directions would mean that the vector R,,, given by R,, = xf + + Z& where x=sin4cos8 y = sin 4 sin 8 and would have to sweep out a complete spherical domain around the hot spot (Figure 4) by varying the values of 4 and 8 in a fixed sequence. The degree of refinement of the sweep depends on the user defined values or quanta by which 4 and 8 are varied, i.e. A 4 and A8, respectively. The number of unique search directions (NUSD) in which the search will proceed in a radiating GLOBALILOCALINTERFACE BOUNDARIES 81 1 Figure 4. Searching in all directions A Dmtion the incremental search proceeds in 4 Nodei 3' r ine(n=O) I0 Figure 5. Incremental search in the desired direction Rgo manner away from the hot spot is given by For each of the desired directions R,,, the OSM procedure performs the search incrementally from one node to the next adjacent node (Figure 5 ) until a boundary node which defines the G/L interface is identified according to the generalized guidelines outlined in Table I. For increment n = 0, the current node inc(n = 0) is set as the node representing the hot spot. For the next increment (n = 1)in the direction represented by R,,, the position vector r from the hot spot to each adjacent node (adjacent to node identified in previous increment, inc(n - 1))is determined. This direction r is then compared to the desired direction R,, and the angle o between the position vector r and the desired direction R,, is determined. The node with the direction vector closest to the desired direction R,, and furthermost from the hot node is then set as the next current node for the nth increment, inc(n). The search procedure is continued step by step in the desired direction R6e until a governing guideline from Table I is satisfied. An example of this search procedure is next presented for clarity. Consider Figure 5 where node 2 has been identified as the current node for increment n = 2 in the direction R,, (i.e. inc(n = 2) = 2). For the next increment (n = 3), all nodes adjacent to node 2 (excluding node 1 to prevent backtracking) are taken into consideration (i.e. nodes 3', 3", 3"'). The position vectors 812 S. SRINIVASAN, S.B. BIGGERS JR. AND R. A. LATOUR JR. r',r'',rmfor each of the three adjacent nodes are next determined with respect to the hot node and w",w"')they make with the desired direction R+e are determined. The the respective angles (a', node with the position vector closest to the desired direction ROO, and secondly, furthermost from the hot node, is then set as the next current node for the third increment (i.e. inc(3) = 37. It is important to note that the OSM considers only the adjacent nodes that are also element corner nodes while searching incrementally, and the incremental search is not allowed to go in loops. In this fashion, the search continues incrementally until one of the generalized guidelines outlined in Table I is satisfied, at which point a boundary node in search direction R+e is said to be reached. Accordingly, if the hot spot is on the surface of the structure, then for a given Rbe: (Al) The search may proceed along the surface until the response gradient tolerance (SGT) limit is satisfied at a node (see next section). This node is then defined as the bounding node in direction R,, and the G/L interface is said to be located in direction R+e. OR (A2) The search may proceed into the interior of the structure where the SGT limit is satisfied at a node (see next section). This node is defined as the bounding node in direction R+, and the G/L interface is said to be located in direction R+e. OR (A3) The search proceeds along the surface until consecutive increments result in increasing deviations from Rbe and the search begins to turn back towards the hot spot, i.e. if ( ( ~ ( n>) o ( n - 1))and (Ir(n)l <(r(n - 1)l)) then the bounding node in direction R+, is the node for which the position vector r made the smallest angle with respect to Rde (i.e. minimum deviation from R+e). The boundary node locates the G/L interface in direction R + g . OR (A4) The search may go in through the interior and then reach another surface of the structure prior to satisfying the SGT limit. If this occurs, the search proceeds incrementally until a node is found which also lies on the surface and for which the position vector r makes the smallest angle with respect to R,,, and thus the least deviation from the desired direction R+o. This node is defined as the bounding node in direction Rge and the G/L interface is said to be located in direction Rbe. If instead, the hot spot is within the interior of the structure, then for a given R+e: (Bl) The search proceeds within the interior of the structure until the SGT limit is satisfied at a node (see next section). This node is defined as the bounding node in direction R,, and the G/L interface is said to be located in direction R+e. OR (B2) The search proceeds until consecutive increments result in increasing deviations from R+, and the search begins to turn back towards the hot spot, i.e. if ((w(n)> w(n - 1)) and (Ir(n)l c ( r ( n - 1)l)) then the bounding node in direction R+e is the node for which the position vector r made the smallest angle or which had the minimum deviation from R+e. The bounding node locates the G/L interface in direction R+@. OR (B3) The search reaches a surface of the structure from within the interior. If this occurs, the search proceeds incrementally along the surface until a node is found which also lies on GLOBALLOCAL INTERFACE BOUNDARIES 813 the surface and for which the position vector r has the minimum deviation from R+e. This node is defined as the bounding node in direction R,, and the G/L interface is said to be located in direction R,,. Conditions (A3) and (A4) when the hot spot is on the surface and conditions (B2) and (B3) when the hot spot lies within the interior of the structure result in the specification of the G/L interface at points in/on the structure where the response gradients have not necessarily reached an acceptably low value. However, they are required to ensure that the OSM does not go into unending incremental search loops, and to avoid misdirecting the search when structural discontinuities such as sharp corners are encountered between increments in the search. If the values AC#Jand A0 by which the search direction is changed are sufficiently fine, such structural discontinuities and sharp changes in structural contour will be adequately captured when searching in other directions. Satisfaction of structural response gradient tolerance (SGT) limit While searching along each direction R,,, the OSM procedure monitors the response gradient, or rate of change in structural response (SED, stress, etc.) with respect to increasing distance after each increment in each search direction. This enables the determination of where the local effects have died out with respect to the user-defined response gradient tolerance (SGT) limit. In order to accomplish this, the change in response for the nth increment AS(n) in direction R,, is determined as AS(n) = S(n) - S(n - 1) where S(n) is the value of the structural response measure at the node identified after the nth increment. The change in distance Ax(n) for the nth increment is defined as Ax(n) = projection of x(n) on R,, (see Figure 6) where x(n) = r(n) - r(n - 1) Finally, the gradient of the response with respect to distance for the nth increment in direction R,, is given by Figure 6. Change in distance as a projection of x(n) onto Roe 814 S. SRINIVASAN, S. B. BIGGERS JR. AND R. A. LATOUR JR. The response gradient "(n) is monitored with respect to the maximum response gradient given by Y',,, which is the maximum value of Y ( i ) ,determined from i = 1, n - 1. If (€+I), is less than or equal to the user-defined SGT limit, then the response gradient is said to have reached an acceptable level for the location of the G/L boundary. The boundary node in direction R,, is then given by inc(n) where n is the increment when @(n) d SGT Some of the possible ways in which response gradients can reach the user-defined SGT limit are presented in Figure 7 where the hypothetical measure of structural response S is plotted as a function of net distance from the hot spot X. Figure 7(a) represents the case when the response gradients decrease gradually away from the concentration or hot spot. The OSM procedure as described above has the ability to account for gradients as shown in Figure 7(a). Figure 7(b) represents the case when the gradients due to one hot spot are interacting with gradients due to another hot spot in the neighbourhood, thus causing zones of localized response perturbation, and Figure 7(c) represents the case when the response experiences some small spikes far away from the concentration such as may occur from approximations inherent in numerical solutions. To account for the cases represented in Figures 7(b) and (c), the search process is allowed to perform n, additional increments (user defined) after the first instance when the normalized response gradient 0 reaches the specified SGT limit. This will ensure that the boundary is not defined erroneously due to a very localized event. If the tolerance level is reached for n, consecutive increments, then the boundary is said to be reached in direction R,, and the boundary node is given by inc(n - n,). Closure to description of method An Objective Search Method (OSM) has been formulated to define an all inclusive and complete interface between global and local domains in 3-D space. The OSM defines the interface by performing the search incrementally in all directions radiating away from the hot spot. The G/L interface is specified in each of the directions R,, by identifying a boundary node as a node where at least one of the generalized guidelines is satisfied. Thus, with adequately defined values of A 4 and At), the OSM determines a unique boundary of nodes representing the complete G/L interface. The degree of uncertainty in defining A 4 and At) in an arbitrary manner is lessened by t t t li (a) Q SGT Limit satisfied in localized .I'.bi 0 SGT Limit satisfied I 1 (b) -X (C) Figure 7. Some possible ways the structural response changes with increasing distance away from the hot spot GLOBAL/LOCAL INTERFACE BOUNDARIES 815 performing the search a number of times and then comparing the NUSD values with the normalized number of unique boundary nodes, [, that result. Convergence in this comparison will give the values of A 4 and A8 which will give a unique boundary of nodes that is all inclusive. However, by examining the extent of the structure being modelled, its aspect ratios, degree of discretization and the aspect ratios of the elements, appropriate values for A 4 and A8 can be chosen without having to perform the above convergence study. Furthermore, if there are multiple concentrations or hot spots which are expected to interact, the OSM can be run starting at each hot spot with a subset of boundary nodes {GL,} resulting for each of the hot spots. An estimation of whether these concentrations interact can be obtained by considering the problem of intersections of subsets as follows: {GL} = {GLi}n{GLj} If {GL} = (01 then hot spots i and j do not interact and the corresponding local domains can be treated independently. Otherwise, the complete boundary of nodes for the G/L interface is given by the union of these subsets, i.e. {GL} = {GLi}u{GLj} ILLUSTRATIVE EXAMPLE PROBLEMS In order to demonstrate the capabilities of the objective search method, two examples of structures with localized stress concentrations are presented. The first example (plate problem) is a plate with a hole under tension in the x-direction. The second example (stem problem) is a 3-D structure which cannot be modelled in two dimensions. This second example problem is a cube with an extended rectangular segment under distributed pressure load. In each example, the measure of structural response being monitored is the Strain Energy Density (SED). For a user-designated SGT limit, the OSM produces a boundary of nodes which identify a unique and reproducible G/L interface. The first step in the OSM is to evaluate parametrically the effect of A 4 and A8 on the normalized number of unique boundary nodes [ that are identified. For this step, the SGT limit is set to an artificially small value of 1 x lo-'' per cent. The values of A 4 and A8 that result in convergence in for each of the two examples are used in the next part of the demonstration. Given the converged values of A 4 and A8 (values of A 4 and A8 resulting in convergence in c), the G/L interface is obtained using the OSM for different values of the SGT limit (e.g. 50,25, 10 per cent, etc.). These steps are performed for the two example problems to demonstrate the versatility of the method. In order to represent graphically the G/L boundary, the global elements which fall within the local domain as defined by the G/L interface are shown following application of the OSM for each example problem. Example I (plate problem, Figure 8 ) The first example problem is a square plate (side = 7 cm)with a hole (diameter = 1 cm)located at its centre under tension (uniform edge displacement of 0.1 cm in the x-direction). This is a classical problem that has been widely modelled"8 and for which the solution is well established. The plate is modelled here in three dimensionswith 20-node quadratic brick elements (Figure s(a)). The objective search method is performed for the plate problem with the SGT limit 816 S. SRWVASAN, S. B. BIGGERS JR. AND R. A. LATOUR JR. (c) SGT-5% (0SGT-so% (e) SGT-25% Figure 8. Graphical representation of isolated global mesh for different values of the SGT limit for the plate problem (same view) set at 1 x lo-'' per cent and the values of A 4 and A9 are varied. Figure 9 shows the effect of the number of unique search directions (NUSD) versus the normalized number of unique boundary nodes, c, identified. From Figure 9, it is clear that convergence in the search is being approached With A 4 and A9 set at So,the SGT limit is when NUSD = 2522 (coincidingwith A 4 = A9 = then varied from 50,25,10,5 to 1 per cent. The sensitivity of the OSM is graphically presented in Figure 10which displays the fraction of the global structure discretization,1,contained within the G/L interface as a function of the SGT limit. The comparison of the number of global elements identified within the G/L boundary for varying values of the SGT limit is made here. A change in the SGT limit produces a change in the size of the local domain with the domain becoming smaller as the SGT limit is increased. Evaluation of the data presented in Figure 10 could provide an indication of the SGT limit for which the complete G/L analysis should be performed. However, the full G/L procedure should be completed and tradeoffs between efficiency and accuracy established before any such indication is conclusive. Graphical representations of the SO). 817 GLOBALFOCAL INTERFACE BOUNDARIES 5 f 0.15 0.10 0.05 0 0 so0 1 o O 0 l u x ) a M ) o u M ) 0 13oM) zdooo 39000 52000 6Ma) NUSD NUSD (b) (a) Figure 9. NUSD vs. number of unique nodes 5 on G/L boundary for the plate problem 0 10 20 30 40 50 SGTLimif (%j Figure 10. Fraction of the number of global elements in local region rl for increasing values of the SGT limit for the plate problem isolated global meshes or of global elements in the local domain for some of these values of the SGT limit are given in Figures 8(a)-(f). Example 2 (stem problem, Figure 11): The second example (stem problem), shown in Figure 11, is a 3-D model of a cube (side = 5 cm) with a rectangular block projecting from the positive y-face. The rectangular block (0.5 x 2 x 1 cm) is located and attached to the face centre of the cube. The opposite face of the cube is clamped. A pressure of 10000 N/m2 is applied on one face of the rectangular projection acting in the positive x-direction. The material properties used are those of stainless steel (E = 200 GPa, v = 0.30). This fully 3-D problem is presented here because the solution and the specification of the G/Linterface is not apparent and poses a difficulty to the analyst. The OSM is performed for per cent and the values of A 4 and A0 are varied. the stem problem with SGT set at 1 x Figure 12 shows the effect of the number of unique search directions (NUSD) versus the 818 S.SRINIVASAN, S. B. BIGGERS JR. AND R. A. LATOUR JR. 9 Y 0.5 cm Figure 11. Schematic illustration of the stem structure and boundary conditions 0.25 O N Figure 12. NUSD vs. number of unique nodes [ on G/L boundary for the stem problem normalized number of unique boundary nodes, C, identified. From Figure 12, it is clear that convergencein search is being approached when NUSD = 2522 (coincidingwith A 4 = A0 = 5"). With A 4 and A0 set at 5", the response gradient tolerance SGT is vaned from 50,25,10,5 to 1 per cent as was similarly done for the plate problem. The comparison of the fraction, A, of global elements identified within the G/L boundary for each value of the SGT limit is given in Figure 13. From this figure, the sensitivity of the OSM can be appreciated. Graphical representations of the isolated global meshes for some of these values of the SGT limit are given in Figures 14(a)-(f). 819 GLOBAL/LOCAL INTERFACE BOUNDARIES 0.50 1 0.25 0 0 10 20 30 SGTlhir(%) 40 50 Figure 13. Fraction of the number of global elements in local region, 1,for increasing values of the SGT limit for the stem problem V (b) SGT-1% (C) SGT-5% Figure 14. (Continued) 820 S. SRINIVASAN. S. B. BIGGERS JR. AND R. A. LATOUR JR. (d) SGT-1096 (c) SCT-25% Figure 14. Graphical representationof isolated global mesh for different values of the SGT limit for the stem problem (two views) DISCUSSION O F COMPUTATIONAL EFFORT OF THE OSM A numerical verification problem of a [ +45/0/90], composite cantilever beam with mid-beam cutouts under pressure (Figure 15) is presented to demonstrate the efficiency and accuracy of the *'~ the OSM). The material properties considered are representaentire G/L p r o ~ e s s ' ~(including tive of a PEEK/Carbon (APC-2/AS-4) composite. Figure 16 shows the global finite element model and the local models (obtained using the OSM with a, being selected as the measure of structural response) used in this demonstration. The global model (Figure 16(a)represents a level of discretization where the global deformation patterns (as a function of model degrees of freedom)have converged to less than 0.5 per cent. The efficiency and accuracy of the G/L process is compared against a traditionally refined model where the refined model has the same number of elements and grading in the x-y plane, i.e. in the plane of the fibres, as the global model. In the out-of-plane direction, the refined model has 16,20-node brick elements representing two 821 GLOBALLOCAL INTERFACE BOUNDARIES t. a Y 1 63.0. 2ti-M.1016. &.W. L l60.8. &.U. P-0.1818 N / m 2 C Figure 15. Schematic illustration of the quasi-isotropic composite cantilever beam with mid-beam cut-outs under pressure load P t’ u , a) Global Model i) sGTIs% ii) SGT-2.5% ii) SGT-75% iv) soT-9096 b) Local Models Figure 16. Finite element meshes of the global model and local models based on different SGT values for the problem of the cantilever beam with mid-beam cut-outs under edge pressure 20-node brick elements per ply? whereas, the global model has four 20-node brick elements through its entire thickness. In order to document the validity of the refined model used as a control, a ‘super-refined’ model was generated by dividing each element belonging to the x-y plane of the refined model into four elements (around four times the number of degrees of freedom of the refined model). The degree of refinement inherent in the traditionally refined model results in the in-plane oxxstresses converging to within 2 per cent at a distance of a lamina thickness away from the free edge at the cut-out. The degree of refinement in the ‘refined’ model is deemed necessary in order to get fairly good solutions in order to facilitate a future experimental verification of the global/local solutions. While the refined model may be made less refined had a different criteria for convergence been considered, the same rules would apply to both global and local level discretizations and thus the computational cost of the G/L approach would also be reduced, including that of the OSM, from the data reported herein. Furthermore, the OSM procedure fits in the realm of ‘processing of global solutions’. No similar consideration has been 822 S. SRINIVASAN, S. B. BIGGERS JR. AND R. A. LATOUR JR given to the significant effort, both computationally and in the man hours required (just considering the size of the data files), in processing the refined model solutions and the efficiency numbers given in the manuscript are, in this sense, conservative and biased towards the traditionally obtained finite element solution. In addition, use of the global/local approach over the traditional approach results in a significant efficiency in terms of the data handling capabilities required, an important consideration while modelling large complex structures. While this may not directly impact the CPU seconds required, it does so indirectly. In fact, use of the global/local approach may allow the solution to be performed in-core and reduces the significant cost associated with swapping, a cost which may be unavoidable in traditional finite element approaches. Part of this cost is evident while comparing the cost in terms of raw CPU seconds versus the model degrees of freedom as reported in Table 11. The local models (Figure 16(b))have the same level of refinement as the traditionally refined model except that they represent isolated discretizations of parts of the beam structure. The lack of symmetry in the local model domains about the cut-out is not surprising considering the increase in moments away from applied edge pressure and towards the clamped end, and the corresponding rate of change of the response. Figures 17(a) and 17(b)show a comparison of in-plane and interlaminar stress plots obtained using the traditionally refined converged model and those obtained from the global model and the local models based on different SGT limits. While there is controversy regarding the reliability of the finite element method in predicting free-edge stresses in composites,16 only a comparison between the traditionally refined finite element model solution and the G/L solution is of interests here and, as such, the comparison is still valid (given that the levels of discretization are similar). Tables I1 and I11 give the computational summary of the G/L method with and without accounting for the OSM. The data, where OSM computational effort is not included, represent the case when the structural analyst specifies a G/L interface which is exactly the same as one specified using the OSM, an unlikely scenario in two dimensions and a nearly impossible scenario in three dimensions. Table I1 gives the degrees of freedom and the raw CPU second count for the refined converged model (used as control), the global model, and the local models for different measures of SGT, the CPU seconds required to perform the OSM at each respective SGT level Table 11. Computationaleffort in terms of degrees of freedom, raw CPU seconds, and accuracy measures (where relevant) for the refined, global, and local models and the corresponding OSM process Name Refined (R) Global Local SGT = 5% SGT = 25% SGT = 75% SGT = 90% OSM SGT = 5% SGT = 25% SGT = 75% SGT = 90% Degrees of freedom CPU (s) 106 035 30 039 120 394 3163 16 161 11 640 6573 5769 1274 510 189 149 YOerror in YOerror in Qxx (msx) Qxz(mar) 24.97 - 2.02 0.63 - 1.07 - 0.90 - - - 70.1 - 1.44 -031 + 0-86 + 1.15 5520 4200 3120 2880 Note: All solutions obtained on a SPARCstation 20, 50 MHz, 256 Mb RAM GLOBAL/LOCAL INTERFACE BOUNDARIES 823 Figure 17(a). Comparison of the global model and local model (defined based on different SGT values) solutions in terms of in-plane u,, stress with the refined model solutions at z / t = 1.0for the cantilever beam problem (note: expanded scale) (b) Figure 17(b). Comparison of the global model and local model (defined based on different SGT values) solutions in terms of out-of-plane uzz stress with the refined model solutions at y / b = 1.0 for the cantilever beam problem and the percent error in ox, and oxrstresses at their respective maxima’s. For each of the local models, the errors in both in-plane ox, and out-of-plane ox, stresses at their respective maxima are less than 2.5 per cent in comparison with the highly refined converged solution. Table 111gives the computational efficiency attained using the G/L method with and without accounting for the OSM. It can be seen from these tables that the computational expense of the OSM is of the same order of magnitude as the global finite element solution. However, the additional expense of the OSM has to be considered in the context that the OSM now allows the complete automation of the G/L process and the G/L method, even after accounting for the OSM, is more efficient than the traditionally refined model by factors of 12 (at worst) and 19 (at best) while the accuracy of stress prediction is within 2.5 per cent. Furthermore, this additional expense due to inclusion of the OSM has to be considered in the following context. (a) The computational cost provided for the OSM reflects the performance of a ‘research code’. There are areas where the code could be optimized resulting in improved computational performance. 824 S. SRINIVASAN, S. B. BIGGERS JR. AND R. A. LATOUR JR. Table 111. Summary of the computational efficiency (CPUefficiency) for different combinations of global and local models, before and after accounting for the OSM Combinations(C) Global + local SGT = 5% SGT = 25% SGT = 75% SGT = 90% Global + local SGT = 5% SGT = 25% SGT = 75% SGT = 90% CPU efficiency R/C 27.13 32.78 3592 36-35 + OSM 12-09 15-29 18.60 19.44 (b) The efficiency numbers reported are subject to the platform upon which they were analysed, which in the current case in a SUN Microsystems SPARCstation 20 (scalar processor). Advantage can be taken of the fact that the search directions in the OSM are independent of one another and the method can be parallelized with minimum engineering cost. This will result in orders of magnitude improvement in wall clock computational times required for performing the OSM. (c) The assumption inherent in the above efficiency measure is that the analyst can specify the same G/L interface as can the OSM. In addition, the expense related to analyst’s time (i.e. for visual examination of 2-D contour plots and 3-D reconstruction of the same) is not accounted for. Furthermore, in the usual scenario (without using the OSM), the G/L interface is specified conservatively by the analyst. While this may not result in a significant increase in computational effort in two dimensions, that is not so in three dimensions, especially while analysing complex heterogeneous composite structures. For instance, if the 3-D global structural discretization is such that the G/L interface is defined as the surface of a sphere of radius r, then any increase dr in r would cause a [(r dr)3 - r 3 ] increase in volume and hence an equivalent increase in the discretization. For thick composite structures, such as one being considered currently for hip prostheses applications (structure with 100-120 plies), the above theoretical consideration is exemplified. Here, an inclusion of the space occupied by one additional global element along the planform results in an additional 200-240, 20-node brick elements having to be considered in the local model. Therefore, the concept of the ‘barest essential’ specification of the G/L interface becomes much more relevant when structures with increasing complexity are considered. As an example, using the numbers provided in Table I1 for an eight-ply laminate, assume that the analyst specifies a G/L interface (without using the OSM) represented by the SGT = 5 per cent case (conservativespecification by the analyst) and the OSM specifies a G/L interface represented by the SGT = 90 per cent case. In this situation, the efficiency of the G/L process is seen to drop by only about 40 per cent when the OSM is used compared to the non-automated, analyst-dependent method of specifying the G/L interface. The accuracy of the prediction remains unaffected (errors in stresses, both in-plane and interlaminar are less than 2.5% in either case). The drop in efficiency attributed to the OSM will be minimal when problems more complex than the simple eight-ply laminated structure are + GLOBALILOCALINTERFACE BOUNDARIES 825 considered, and in all probability, greater efficiency (as opposed to a drop in efficiency) could be attained through use of the OSM for problems such as the one involved with designing a composite hip prosthesis. (d) In cases when both global and local models need to reflect design changes or changes in load envelopes, the OSM provides a consistent means to account for interaction between global and local models and, in certain cases, the OSM may need to be performed only once. (e) Lastly, and probably most importantly, implementation of the OSM technique provides the necessary connectivity between global and local models, enabling G/L analysis to be performed in a fully automated, analyst-independent manner. It is illustrative to compare the performance of G/L method along with the newly proposed OSM against traditional processes such as adaptive mesh refinement for the practical analysis of engineering structures. The concepts of adaptive mesh based on criteria such as energy norms and others, are used to provide guidance for remeshing to attain a sufficient level of discretization and to predict structural response with some desired accuracy. Adaptive refinement is usually carried out in the context of some error estimate of either independent or dependent fields in the finite element formulation so as to control and reduce to within acceptable levels, the errors associated with the finite element solution. The general theme is to obtain an ‘optimal mesh’ after i iterations of the adaptive refinement process where the contribution to the overall error in the discretization is equal for all elements, and where each of the contributions are below an analyst-specified value, i.e. the concept of equidistribution of error.10-12This is quite a stringent requirement to satisfy for the entire mesh, especially if the structure being modelled exhibits localized behaviour in the form of multiple stress concentrations which may or may not interact as is the case with the class of problems associated with the design and analysis of a composite hip prosthesis. In fact, this would result in highly complex modelling and possibly distorted elements (contributing to error), require transitions between areas of high refinement and those with a coarser grid, and if nodal compatibility is to be maintained, will result in the propagation of highly refined areas through at least a part of the rest of the structural discretization even if concepts such as 1-irregular connections and transition or ‘green’ elements’ are used. Global/ local methods (based on the specified boundary displacements principle’) do not require such nodal compatibility and allow different levels of discretization in the global and local domains. Furthermore, it is well known that in displacement-based formulations, the displacements converge at a faster rate than dependent variables like stresses with respect to mesh refinement/enrichment. This fact is taken advantage of in discretizing the global domain with only a coarse mesh being necessary before the global deformation patterns converge. The local domain@ can be discretized independently to get sufficiently accurate stresses, etc. The OSM provides a consistent basis to account for the interaction between the global and local domains in an automated, objective fashion. The coarse global mesh and the independently discretized local domains result in the G/L method being far more efficient than the solution obtained using a mesh arrived at after the adaptive refinement process for a given level of solution error. Given the fact that the G/L interface is usually specified conservatively by the analyst, use of the OSM provides a certain amount of efficiency even after considering its cost (this argument has been detailed previously). Therefore, the G/L process with the OSM has the potential to be competitive with traditional adaptive refinement techniques and depending upon the model complexity (e.g. structure with multiple stress concentrations), may be substantially more efficient. Through its ability to define the ‘barest essential’ zone, the OSM is necessary to solve the class of problems such as involved with designing and analysing thick composite structures for hip prosthesis applications. This is due to the fact that at least two fully integrated 20-node brick 826 S. SRINIVASAN, S. B. BIGGERS JR. AND R. A. LATOUR JR. elements are required per ply thickness to capture adequately the response of even a simple cross-ply composite plate subjected to in-plane loads4 Even if the same were assumed to be true for the far more complex internal design of the composite hip stem which has not only in-plane loads but also in-plane and out-of-plane bending moments and twisting moments imposed, this would require about 200-240,20-node brick elements just through the thickness of the 100-120 ply composite hip stem. It has been estimated that modelling such a system adequately using standard finite element technology (h-refinement, p-enrichment, r-relocation), and maintaining reasonable element aspect ratios and tapers in order to get an 'optimal' mesh, would involve solving a problem with around 40 million degrees of freedom. It is doubtful if the problem of this magnitude can be solved at all using the best of today's computational resources, even those offered by supercomputing platforms. In the context of analysing thick composites (100-120 plies) for hip prosthesis applications and their associated free edge problems, it will be impossible to use anything but the more efficient G/L method in conjunction with the newly introduced OSM which provides a consistent basis to account for the interaction between global and local domains and provides the means to completely automate the G/L process. CONCLUDING REMARKS An automated and objective search method (OSM) has been developed to identify the G/L interface in 3-D space. Also incorporated is an unique numerical method by which surface nodes can be identified. Given a user-defined set of search increments A 4 and A0 and a user-defined structural response gradient tolerance (SGT) limit, the OSM identifies a unique and reproducible G/L boundary. This G/L boundary is all inclusive and complete for the given level of global discretization. While the OSM has the capability to handle relatively simple problems (plate, stem and beam problems), the rationale used to develop the OSM could give it the capability to handle complex structural geometries and complex structural response states. Inherent within the OSM is the ability to identify a unique boundary of nodes surrounding the hot spot(@where the selected measure of structural response has satisfied a user-defined SGT limit. However, the density or the number of interfacial nodes identified by the OSM procedure depends on the mesh density and mesh grading at the global level. Moreover, depending on the interpolation or coupling scheme used to proceed from the global to the local models, the OSM-defined G/L interface may need to be modified to a smooth or a regular boundary. However, all of the OSM-defined interfacial nodes should be included in any such modification. The OSM has been developed for intended use with G/L methods based on the specified boundary displacement method' (SBD) where a priori knowledge of the localized areas of interest is not required. The SBD-based G/L methods are essentially uncoupled G/L methods where both global and local models and their levels of discretization are independent of one another and the OSM is best suited for use with these methods. On the other hand, coupled G/L methods in addition to requiring a priori knowledge of the localized areas of interest also require that the local effects be accounted for or embedded, at least in an approximate manner, in the global model and some even require that there be an one-to-one nodal correspondence at the G/L interface. The OSM can be used even with coupled G/L methods to provide a consistent basis to account for interaction between global and local domains. The penalty, however, is that the global model needs to be solved twice, first to obtain knowledge about the localized areas of interest and their possible interactions and, second, with the approximate local embedment included. While this method can also be easily adapted to problems in two dimensions to define GLOBALLOCAL INTERFACE BOUNDARIES 827 the G/L interface more objectively, it is particularly suited to problems where both the global and local models are in three dimensions. The most important aspect of this newly developed method is that the OSM provides the necessary link which may now allow the complete automation of the G/L process involved with the engineering analysis of complex structures in today's computing environment. Independent global and local models can now, through use of the OSM, be developed and their interactions accounted for in a consistent and automated manner. The state-of-the-art would then be to perform error analysis and adaptive refinement at the global level in an automated fashion, use the OSM to define automatically the interface between global and local domains and finally perform error analysis and adaptive refinement at the local level in an automated manner. The complete automated process with minimal analyst intervention would satisfy the above ideal scenario and provide for efficient accurate analysis of complex engineering structures. ACKNOWLEDGEMENT The authors would like to thank the Whitaker Foundation, Washington D.C., for funding this research. APPENDIX Notation n increment count S user-defined measure of structural response (e.g. strain energy density, components of stress) user-defined response gradient tolerance limits (per cent) polar radius of oblique spherical triangle surface area of ith spherical triangle number of adjacent elements total surface area of N spherical triangles unit position vector with respect to hot node in the desired search direction 4, 0 angle made by R,, with the z-axis angle made by the projection of R,, on the x-y plane with the x-axis user-defined change in 4 user-defined change in 0 number of unique search directions node identified for the nth increment of the search in any 4, 0 direction position vector of inc(n) with respect to hot node angle between r(n) and R+e (= r(n) - r(n - 1)) projection of x(n) onto R,@ absolute rate of change of response with respect to distance (lAS/Axl); response gradient response gradient convergence parameter (per cent) number of additional consecutive increments where gradient convergence must be satisfied number of unique boundary nodes normalized with respect to total number of nodes in global model 828 S.SRINIVASAN, S. B. BIGGERS JR. A N D R. A. LATOUR JR. I number of elements within/on G/L boundary as a fraction of total number of elements in global model {GL,} set of boundary nodes for ith hot spot {GL} set of boundary nodes for all interacting hot spots REFERENCES 1. J. A. Davidson, ‘The challenge and opportunity for composites in structural orthopaedic applications’, J. Compos. Technol. Res., 9, 151-161 (1987). 2. A. K. Noor, ‘Global-local methodologies and their application to nonlinear analysis’, Finite Element Anal. Des., 2, 333-346 (1986). 3. 0. H. Griffin Jr., ‘The use of computers in the evaluation of three-dimensional stress effects in composite materials products’, Keynote paper, proc. 2nd In?. Con$ Comput. Aided Des. Compos. Mater. (CADCOM), Brussels, Belgium, 25-27 April, 1990. 4. D. M. Thompson and 0. H. Griffin Jr., ‘2-D to 3-D global/local finite element analysis of cross-ply composite laminates’, J. Rein$ Plast. Compos., 9, 492-502 (1990). 5. I. Hirai, B. P. Wang and W. D. Pilkey, ‘An efficient zooming method for finite element analysis’, Int. j . numer. methods eng., 20, 1671-1683 (1984). 6. C. C. Jara-Almonte and C. E. Knight, ‘The specified boundary stiffness/force SBSF method for finite element subregion analysis’, Int. j . numer. methods eng., 26, 1567-1578 (1988). 7. J. B. Ransom and N. F. Knight Jr., ‘Global/local stress analysis of composite panels’, Comput. Smct., 37, 375-395 (1990). 8. 0. H. Griffin, Jr. and M. A. Vidussoni, ‘Global/local finite element analysis of composite materials’, Proc. Int. Con$ Comput. Aided Des. Compos. Mater. Technol., Southampton, 1988, pp. 513-523. 9. 1. B. Ransom, S . L. McCleary, M. A. Aminpour and N. F. Knight Jr., ‘Computational methods for global/local analysis’, NASA TM-107591, 1992, pp. 1-22. 10. 0.C. Zienkiewicz and J. Z. Zhu, ‘A simple error estimator and adaptive procedure for practical engineering analysis,’ Int. j . numer. methods eng., 24, 331-357 (1987). 11. T. Strouboulis and K. A. Haque, ‘Recent experiences with error estimation and adaptivity, Part I: review of error estimators for scalar elliptic problems,’ Comput. Methods Appl. Mech. Eng. ASME, 97, 399-436 (1992). 12. T. Strouboulis and K. A. Haque, ‘Recent experiences with error estimation and adaptivity, Part 11: error estimation for h-adaptive approximations on grids of triangles and quadrilaterals,’ Comput. Methods Appl. Mech. Eng. ASME, 100,359-430 (1992). 13. H. Eves, Spherical triangles, in W. H. Beyer (ed.), CRC Standard Mathematical Tables, 27th edn., CRC Press, Boca Raton, FL, 1984. 14. S. Srinivasan, S. B. Biggers Jr. and R. A. Latour Jr., ‘The 3-D global/3-D local method as a tool for efficient analysis and design of FRP composite hip implants,’ Trans. 2lst Annual Mtg. Socfor biornaterials, in Conjunction with 27th Int. Biomaterials Symp.. 342, San Francisco, CA. 18-22 March 1995. 15. S. Srinivasan, S. B. Biggers Jr. and R.A. Latour Jr., ‘Analysis of composite structures using the 3-D global/3-D local method,’ Proc. 10th Tech. Con$ Compos. Mater., Am. SOC.Compos., Santa Monica, CA, 18-20 October 1995. accepted. 16. J. D. Whitcomb, I. S. Raju and J. G . Goree, ‘Reliability of the finite element method for calculating free edge stresses in composite laminates,’ Comput. Struct., 15, 23-37 (1982).

1/--страниц