Lei Li
email: leili@freescale.com

Krishnendu Chakrabarty

Hybrid BIST Based on Repeating Sequences and Cluster Analysis 1

We present a hybrid BIST approach that extracts the most frequently occurring sequences from deterministic test patterns; these extracted sequences are stored on-chip. We use cluster analysis for sequence extraction, and encode deterministic patterns on the basis of the stored sequences. Experimental results for the ISCAS-89 benchmark circuits show that the proposed approach often requires less on-chip storage and test data volume than other recent BIST methods.

Introduction

High test data volume and limited tester channel bandwidth are two major problems encountered in the testing of today's system-on-chip integrated circuits. In order to mitigate these problems, a number of techniques based on test data compression, built-in self-test (BIST), and a combination of the two have been proposed in the literature.

In the test data compression approach, a deterministic test set is compressed and stored in tester memory. The compressed test set is transferred through tester channels to the circuit under test (CUT), where it is decompressed by decoding hardware. Test compression techniques based on on-chip pattern decompression are presented in [START_REF] Chandra | Test data compression and test resource partitioning for system-on-a-chip using frequency-directed runlength (FDR) codes[END_REF][START_REF] Gonciari | Improving compression ratio, area overhead, and test application time for system-on-a-chip test data compression/decompression[END_REF][START_REF] Li | Test data compression using dictionaries with fixed-length indices[END_REF][START_REF] Rajski | Embedded deterministic test for low-cost manufacturing test[END_REF][START_REF] Wuertenberger | Data compression for multiple scan chains using dictionaries with corrections[END_REF]. In BIST solutions, test patterns are generated by an on-chip pseudo-random pattern generator, usually a linear-feedback shift-register (LFSR). A number of BIST techniques based on test point insertion [START_REF] Schotten | Test point insertion for an area efficient BIST[END_REF], reseeding [START_REF] Hellebrand | A mixed-mode BIST scheme based on reseeding of folding counters[END_REF][START_REF] Liang | Two-dimensional test data compression for scan-based deterministic BIST[END_REF][START_REF] Rajski | Test data decompression for multiple scan designs with boundary scan[END_REF], bit-flipping [START_REF] Wunderlich | Bit-flipping BIST[END_REF], and weighted random pattern testing [START_REF] Wang | Low hardware overhead scan based 3-weight weighted random BIST[END_REF] have been proposed. Deterministic test patterns are applied in BIST by either controlling the state of the pattern generator [START_REF] Hellebrand | A mixed-mode BIST scheme based on reseeding of folding counters[END_REF][START_REF] Krishna | Test vector encoding using partial LFSR reseeding[END_REF][START_REF] Liang | Two-dimensional test data compression for scan-based deterministic BIST[END_REF][START_REF] Rajski | Test data decompression for multiple scan designs with boundary scan[END_REF] or by altering the output of the pattern generator [START_REF] Wang | Low hardware overhead scan based 3-weight weighted random BIST[END_REF][START_REF] Wunderlich | Bit-flipping BIST[END_REF].

Techniques based on the combination of data compression and BIST have also been developed recently [START_REF] Jas | Hybrid BIST based on weighted pseudo-random testing: a new test resource partitioning scheme[END_REF][START_REF] Krishna | Reducing test data volume using LFSR reseeding with seed compression[END_REF]. The hybrid BIST scheme presented in [START_REF] Jas | Hybrid BIST based on weighted pseudo-random testing: a new test resource partitioning scheme[END_REF] applies weighted pseudo-random patterns to the circuit to achieve 100% fault coverage. The compressed weight set is stored on ATE and decompression is carried out using an on-chip look-up table. In [START_REF] Krishna | Reducing test data volume using LFSR reseeding with seed compression[END_REF], the seeds for the LFSR are compressed using statistical coding.

In this paper, we present a hybrid BIST approach that extracts the most frequently occurring sequences from deterministic test patterns; these extracted sequences are stored onchip. The test session consists of three stages. In the first stage, pseudo-random patterns generated by an LFSR are applied to the CUT to detect easy-to-detect faults. In the second stage, semi-random patterns are generated by randomly selecting some of the stored sequences and flipping some of their data bits. Since the semi-random patterns are generated based on the sequences extracted from deterministic patterns, they are more likely to detect the hard-to-detect faults than the pseudo-random patterns. Faults that are not detected by pseudo-random and semi-random patterns are detected by deterministic patterns in the third stage. The deterministic patterns are encoded on the basis of the stored sequences to reduce the test data volume, and are decoded/generated during the third stage of the test session. The overall approach is similar to [START_REF] Wuertenberger | Data compression for multiple scan chains using dictionaries with corrections[END_REF], where repeating scan slices are stored in a dictionary, and corrections (bit-flips) are applied to selected scan slices.

The rest of this paper is organized as follows. Section 2 presents the proposed hybrid BIST approach and its associated synthesis procedure. The cluster analysis algorithms used in the synthesis procedure are described in Section 3. Section 4 describes the proposed BIST architecture. Section 5 presents experimental results for the ISCAS-89 benchmark circuits. Finally, Section 6 concludes the paper.

Proposed Approach

The proposed mixed-mode BIST approach is based on the observation that identical or similar sequences often appear in many test patterns that are applied to a logic circuit. Consider the simple circuit shown in Figure 1. In order to detect the stuck-at faults d/0, d/1, e/0, e/1, g/0, and g/1, f needs to be set to 1 to propagate the faulty values to the output, which requires that the primary inputs "abc" be set to "111". Thus the sequence "111" appears at the primary inputs "abc" in the six test patterns to detect the above stuck-at faults. Based on this observation, we extract a number of more frequently-occurring sequences from the deterministic test patterns for a logic circuit and store them on the chip. By selecting sequences and flipping some of their data bits, we can either generate semi-random patterns or encode deterministic patterns, based on how we select the sequences and the bits to be flipped.

Similar to other BIST approaches, we first apply a number of pseudo-random patterns to the CUT to detect easy-todetect faults. Next we generate deterministic patterns for the remaining faults, and extract a number of frequently occurring sequences from these deterministic patterns. The extraction of repeating sequences is carried out in two steps. First, the data bits in the test vectors are reorganized such that the positions that are specified in the same test vectors are grouped together. This is illustrated in Figure 2, where positions 1, 2, 3 and 7 are grouped, and positions 4, 5 and 6 are grouped. Corresponding to the reorganization of the data bits in the test set, the scan cells in the scan chain also need to be reorganized. Repeating sequences are next identified for each group. As shown in the figure, four sequences are extracted for the first group. Finally, we need 23 bits to store the sequences (the x in the last group does not need to be stored). In other words, the test set can be viewed as a matrix in which each row is a test vector. First, columns are grouped/reorganized to bring together the specified bits in a row. Then the grouped columns are regarded as a submatrix and repeating sequences/rows are determined.

In the simple example of Figure 2, the position groups and the repeating sequences are easy to identify, and the number of groups and sequences are small. However, a test set for a real-life circuit can lead to a large number of groups and sequences. We need to reduce the number of groups such that the encoding of a deterministic pattern is beneficial, and we need to merge some of the sequences in each group to make the storage requirement manageable. We use cluster analysis for both position grouping and sequence merging; the details of these procedures are discussed in Section 3.

After grouping the positions in the scan chain and extracting the most occurring sequences in each group, we generate semi-random patterns based on the extracted sequences. Assuming that the positions in the scan chain are grouped into G groups, and R sequences are identified for each group, we use log 2 R outputs from the LFSR to generate a random number r, 0 ≤ r < R, which selects a sequence from each group. Several outputs of the LFSR are then AND-ed to generate a flip indication sequence with a small proportion of 1s. The semirandom patterns shifted into the scan chain are obtained by doing an exclusive-or between the randomly selected sequences 0 00 00 1 00 1 01 00 1 10

Not the last group

The last bit

The last bit and the flip indication sequence.

After a number of semi-random test patterns are generated, we run fault simulation to determine the faults that are not detected by both the pseudo-random patterns and the semi-random patterns. Then we use an ATPG tool to generate deterministic test cubes for the remaining undetected faults. These deterministic cubes are then encoded by indicating the selection of sequence and bits to be flipped. For example, if we need to encode the cube "x000011x1" based on the sequence obtained in Figure 2(c), we first reorder the sequence to "000x", "011", "1" and "x", then encode it as "sequence(0, 0) flip(0), sequence(1, 0) flip( 2)", where sequence(i, j) means selecting the jth sequence in the ith group, and flip(k) means flipping the kth bit. In this example, the number of groups is 4, the maximum number of sequence in a group is 4, and the maximum number of bits in a sequence is 4. So we need 2 bits each to encode the group index, sequence index and the bit index. Since multiple groups might be needed for encoding a deterministic pattern, and multiple bits might need to be flipped in one sequence, a prefix is needed before the group index and the bit index to indicate if it is the last group for this pattern or if it is the last flipping bit in the current sequence. The encoded data for the deterministic pattern "x000011x1" is shown in Figure 3.

Cluster Analysis

Cluster analysis is used for numerical classification in many fields [START_REF] Everitt | Cluster Analysis[END_REF]. We use an agglomerative clustering algorithm for both position grouping and sequence merging [START_REF] Everitt | Cluster Analysis[END_REF]. The test set is first represented as a graph, to which the agglomerative algorithm is applied. The basic idea is to first view each node in the graph as a cluster, then continuously merge pairs of closest clusters until a predefined threshold is reached. This threshold can be either the number of clusters in the system, referred to as maxClusterN um, or the minimum distance between a pair of clusters that can be merged, referred to as minM ergeDistance.

The details of the agglomerative algorithm are described in Figure 4. The input to the procedure is the system S containing a given number of clusters, where each cluster corresponds to a node in the graph. The distance between each pair of clusters is calculated and stored in the distance matrix D. In each loop of the procedure, the minimum distance and its associated cluster indices i and j are found. If the minimum distance is not larger than the preset threshold minM ergeDistance, the clusters i and j are merged and j is deleted from the graph and the distance matrix. All the elements in the distance matrix related to i are also updated. The procedure terminates when either maxClusterN um or minM ergeDistance is reached.

Although the agglomerative algorithm is used for both position grouping and sequence merging, the calculation of the distance between a pair of clusters is different for the two problems. In position grouping, there are L nodes in the graph, where L is the length of the scan chain, and each node represents a position in the scan chain. Each node j has an associated index set I j that records all the test vectors for which this position is specified. The distance between two nodes i and j is calculated as the number of test vectors in which one node/position is specified while the other is unspecified, i.e., it simply equals |(I i -I j ) (I j -I i )|. As described above, the initial cluster system contains L clusters, each of which consists of one node. In order to calculate the distance between two clusters, we define the centroid of a cluster as follows. The centroid of a cluster is a dummy node whose associated index set is the union of the associated index sets of all the nodes contained in this cluster. The distance between two clusters is calculated as the distance between their centroids of the two clusters. We adopt this definition of a centroid because we are trying to group the positions that are specified in the same test vectors. After several positions are merged into a group, the corresponding group appears in a test vector as long as one of its position is specified.

Figure 5 illustrates position grouping. We use the same test data as in Figure 2. However, we set the threshold maxClusterN um to 2, excluding a group that contains all the positions that are unspecified in all the test vectors. In this example, position 0 is unspecified in all the 7 test vectors; therefore it is not considered for cluster analysis. As shown in After position grouping, we reorder the columns of the test data matrix to bring the positions/columns in a group together. The basic idea was illustrated in Figure 2. For each group, i.e., a submatrix, an initial cluster system is generated and the agglomerative algorithm is used for sequence merging. The number of nodes for a group is equal to the number of test vectors in which the group appears, i.e., the number of rows in the submatrix that contain at least one specified bit. A 3-valued sequence (with elements from {0, 1, x}) is associated with each node. The distance between a pair of nodes is the number of bits that are in conflict for the two sequences that are associated with them. Two bits are in conflict with each other if both of them are specified and they are set to different values. In the initial cluster system, each cluster contains exactly one node.

We define the centroid of a cluster as a dummy node associated with a three-value sequence, which is obtained as follows. The length of the sequence is the same as the sequences associated with the nodes in this cluster. For each bit in the string, we count the number of 0s (zeroCounts) and the number of 1s (oneCounts) in this bit position for all the sequences associated with the nodes in this cluster. If both numbers are equal to 0, this bit is set to 'x'. Otherwise this bit is set to 0 if zeroCounts ≥ oneCounts, and 1 if zeroCounts < oneCounts. The distance between a pair of clusters is then calculated as the product of the distance between the centroids of the pair of clusters and the total number of nodes in these two clusters. The above method for determining the centroid of a cluster and the distance between a pair of clusters has been developed to ensure that only a small number of bits need to be flipped to generate all the sequences in a cluster from the centroid sequence of the cluster.

Figure 6 illustrates cluster analysis for sequence merg- ing. The input data in this example corresponds to the first group of sequences in Figure 5(h). We set the threshold maxClusterN um to 2. Figure 6(b) shows the initial clusters and the sequences associated with their centroids. The distance matrix is shown in Figure 6(c). Note that the distance between a pair of clusters is the product of the number of conflict bits with the number of nodes in this pair of clusters. Thus although there are 2 conflict bits between Clusters 1 and 2, their distance is 4 as shown in the distance matrix. First, Clusters 1 and 3 are merged because the distance between them (2) is the smallest entry in the distance matrix. After merging Clusters 1 and 3, the distance matrix is updated. Then Clusters 2 and 4 are merged. This concludes the cluster analysis because the threshold maxClusterN um is reached. The final sequences extracted for this group are "1001x" and "10101".

Proposed BIST Architecture

Figure 7 shows the proposed BIST architecture. We first describe the architecture for a CUT with a single scan chain. We then discuss the extension for multiple scan chains. For a single scan chain, the signals connected to the two multiplexers and the exclusive-or gate are all 1-bit wide.

In the first stage of the test session, the signal Select random is set to 0. Hence the pseudo-random patterns generated the LFSR are shifted into the scan chain through MUX I and applied to the CUT. In this stage, Bit counter A is also used to indicate the end of each test pattern.

In the second stage, the signal Select random is set to 1, and Select f lip is set to 0. The test data is obtained by doing an the exclusive-or between the F lip indication R signal and the ROM output. The signal F lip indication R is obtained by AND-ing several bits from the LFSR. Thus it contains much more 0s than 1s, a feature that is used to flip some of the bits in the data sequence obtained from the ROM. The ROM contains R × C bits organized as a matrix [M ] R×C with R rows and C columns. Let the content of the ROM for the location addressed by the pair (i, j), 0 ≤ i < R and 0 ≤ j < C, be M i,j . The number of columns C in the ROM is often less than the scan chain length L. Note that the Column select signal comes from Bit counter A and it can be larger than C. In this case, the ROM simply outputs the value 0. The data in the ROM is divided horizontally into G groups, where each row in a group corresponds to an extracted sequence. The Group end decoder always outputs 0, except at the last bit of each group, when it outputs 1 to indicate the end of the current group. The signals Load r and Select r are then set to 1 to load the random number from the LFSR into the Row select register such that a random sequence is selected from the next group.

In the third stage, both Select random and Select f lip are set to 1. Thus the test data shifted into the scan chain is obtained via an exclusive-or operation between F lip indication D signal and the ROM output. F lip indication D is generated by the finite-state machine based on the encoded data from Data in. The Group counter contains an extra buffer besides the standard counter in it. The extra buffer is used to store the data that is shifted in from Data in while the standard counter is operating. With the signal Load g set to 1, the data stored in the buffer is loaded into the standard counter, which is then ready for counting down. The same pipeline structure is used for Bit counter B and Row select register. The Group counter is used to indicate whether the selected group is reached. After the selected group is reached, the Row-select register is loaded with the selected sequence index which is stored in its buffer. The Bit counter B is also loaded with the flipping bits index stored in its buffer and it then starts to count down. When the Bit counter B decrements to 0, F lip indication D is set to 1 to flip the data bit from the ROM. The Data input counter is used to count the number of bits shifted from Data in and indicate whether the data for the group index, sequence index or bit index has been completely shifted into its buffer.

For multiple scan chains, the ROM is reorganized as shown in Figure 8, which corresponds to the case of 4 scan chains for the CUT. The number of banks in the ROM equals the number of scan chains, i.e., four in this example. In each scan cycle, one bit from each bank is shifted into the scan chain. All of the four bits are from the same position of the banks, which is indicated by the same column select and row select signals. The organization of the sequences is also shown in the figure, where b i,j denotes the jth bit of the ith sequence in a group, the groups of the sequences are divided by the solid line, and the sequences in each group are divided by the dashed line. As shown in the figure, the data bits of the ith sequence are placed in the ith row of the ROM banks, and the current columns of all the banks are filled before proceeding to the next columns. The number of bits in a group may not always divide the number of ROM banks (the number of scan chains) exactly, thus some cells in the last column of this group need to be filled randomly. The cells in the scan chains also need to be reorganized, i.e., the cells belonging to the first group are first placed at the first positions of the scan chains, then at the second positions of the scan chains, and so on. Figure 9 shows the reorganization of the scan cells corresponding to the data organization in the ROM shown in Figure 8. In Figure 9, Group 0 contains 10 scan cells. The first 4 cells are placed at the first position of each scan chain, the next 4 cells are place at the second position of each scan chain, and the remaining 2 cells are placed at the third position of scan chain 0 and scan chain 1. Scan cells that are unspecified in all the deterministic patterns are used to fill the third position of scan chain 2 and scan chain 3.

For a CUT with m scan chains (m > 1) and the above organization of the ROM, the decoding architecture operates in nearly the same fashion as a single scan chain. The differences are as follows. First, the output from the ROM, the signals F lip indication R and F lip indication D, and the input signals to MUX I are all m-bit wide. Second, in the third stage of the test session, recall that for a single scan chain, the bit in- We first implemented the BIST architectures for a single scan chain using the lsi 10k library of Synopsys Design Compiler. Using the wire load model for the lsi 10k library, we designated the normalized area for a unit-length of wire to be 0.2 (assuming that the area of an inverter is 1 unit) to take into account the additional area due to interconnects. The area overhead for the FSM, measured in Synopsys gate equivalents, is 187.97, which includes 10 flip-flips and 61 logic gates. The maximum area overhead for the Group-end decoder is 120.53, which includes 65 logic gates (for circuit s38417). The hardware overhead for the architecture shown in Figure 7, excluding the LFSR and the ROM, is 834.07 for the circuit s38417, including 47 flip-flops and 158 gates. This amounts to only 2.04% overhead.

To extend the BIST architecture to multiple scan chains, the FSM needs to be modified to generate m data bits in each scan clock cycle, where m is the number of scan chains. The data in the ROM needs to be reorganized as discussed above, but without any change in the storage requirements. The other components are the same as for the single scan chain architecture. The hardware overhead for the extended architecture for four scan chains, excluding the LFSR and the ROM, is 870.95 for the circuit s38417, including 49 flip-flops and only 171 gates.

Experimental Results

In this section, we present experimental results for the seven largest ISCAS-89 benchmark circuits. Table 1 presents the first set of results. We first apply 10000 pseudo-random patterns to the CUT. Next we use the ATPG tool Atalanta to generate a set of deterministic test patterns for the remaining faults. The number of extracted sequences is four for all the circuits. After sequence extraction, 10000 semi-random patterns are generated and applied to the CUT. Finally, we use Atalanta to generate deterministic patterns for the remaining faults and we encode the deterministic patterns based on the extracted sequences.

The number of bits needed for storing the extracted sequences on-chip is listed in the seventh column. For the circuits s5378 and s38584, the amount of stored data is less than the test data volume corresponding to just one test pattern. For the two worst cases, i.e., for s9234 and s38417, the on-chip storage requirements are less than the test data volume for three test patterns. For the circuits s5378 and s15850, only three and two deterministic patterns need to be encoded, respectively. The CPU times for the computation range from 7 minutes to 19 minutes on a 1.4 GHz Pentium 4 PC with 512 MB of memory.

Only a small number of deterministic test patterns need to be encoded, and only these patterns are required to be fed through the tester channel. The pseudo-random test and semi-random test can be run at higher speed since no data is required from the tester in these stages.

Table 2 compares the storage requirements of the proposed approach with test vector encoding using partial LFSR reseeding [START_REF] Krishna | Test vector encoding using partial LFSR reseeding[END_REF], BIST based on reseeding of folding counter [START_REF] Hellebrand | A mixed-mode BIST scheme based on reseeding of folding counters[END_REF], and two-dimensional test data compression [START_REF] Liang | Two-dimensional test data compression for scan-based deterministic BIST[END_REF]. The storage requirement reported for the proposed approach include the data for storing the extracted sequences and the data for encoding the deterministic patterns. The results presented in the literature for these methods also rely on 10000 initial pseudorandom patterns to eliminate the easy-to-detect faults. The proposed approach requires less storage than partial LFSR reseeding method [START_REF] Krishna | Test vector encoding using partial LFSR reseeding[END_REF]. Compared to BIST based on reseeding of folding counter, the proposed method provides better results in three out of six cases. Note however that width compression is used in [START_REF] Hellebrand | A mixed-mode BIST scheme based on reseeding of folding counters[END_REF] to reduce test data volume. While width compression can indeed reduce test data volume, it requires a special scan out procedure is needed to shift out the test responses. Compared to two-dimensional test data compression, the proposed approach provides better results in four out of six cases.

Table 3 compares the proposed approach with hybrid BIST based on weighted pseudo-random patterns [START_REF] Jas | Hybrid BIST based on weighted pseudo-random testing: a new test resource partitioning scheme[END_REF]. For both methods, the data volume is divided into two parts, on-chip storage and encoded test data. The results of [START_REF] Jas | Hybrid BIST based on weighted pseudo-random testing: a new test resource partitioning scheme[END_REF] relies on 32000 initial pseudo-random patterns to eliminate easy to detect faults. Thus we also applied 32000 pseudo-random patterns to the CUT in the first stage of the test session. The number of semi-random patterns applied to the CUT is kept at 10000. Compared to hybrid BIST based on weighted pseudo-random patterns, the proposed approach requires less on-chip storage for the larger circuits although it needs slightly more on-chip storage for the smaller circuits. The encoded data volume of the proposed approach is less than that of [START_REF] Jas | Hybrid BIST based on weighted pseudo-random testing: a new test resource partitioning scheme[END_REF] for all but one circuit.

We have implemented the BIST logic using Synopsys tools and reported the results in Section 5.3. The area overhead appears to be reasonable and, and based on limited published data, of the same magnitude as the other four methods.

Conclusion

We have presented a new hybrid BIST approach that uses cluster analysis to extract the most frequently occurring sequences deterministic test patterns. Experimental results for the ISCAS-89 benchmark circuits demonstrate that the proposed approach often requires less on-chip storage and test data volume than recently-proposed methods.

  Figure 1. A simple circuit used to motivate the clustering-based approach.

Figure 2 .

 2 Figure 2. An example to illustrate the extraction of frequently occurring patterns.

Figure 3 .

 3 Figure 3. The encoded data for the deterministic test cube "x000011x1".

Figure 4 .

 4 Figure 4. The agglomerative algorithm.

  Figure 5(b), the initial system consists of eight clusters, each of which contains one position. The test vector index set associated with the centroid of each cluster is also listed. The distances between each pair of clusters are listed in Figure 5(c). First, all pairs of clusters with mutual distance 0 are merged; the remaining clusters and their pairwise distances are shown in Figures 5(d) and 5(e). Next, the Clusters 1 and 8 are merged because the distance between them is the minimum entry in

Figure 5 .

 5 Figure 5. Clustering for position grouping.the distance matrix. After the merging of clusters 1 and 8, the maxClusterN um threshold is reached. Figures5(f) and 5(g) show the eventual outcome of cluster analysis. The positions are then grouped as shown in Figure5(h).After position grouping, we reorder the columns of the test data matrix to bring the positions/columns in a group together. The basic idea was illustrated in Figure2. For each group, i.e., a submatrix, an initial cluster system is generated and the agglomerative algorithm is used for sequence merging. The number of nodes for a group is equal to the number of test vectors in which the group appears, i.e., the number of rows in the submatrix that contain at least one specified bit. A 3-valued sequence (with elements from {0, 1, x}) is associated with each node. The distance between a pair of nodes is the number of bits that are in conflict for the two sequences that are associated with them. Two bits are in conflict with each other if both of them are specified and they are set to different values. In the initial cluster system, each cluster contains exactly one node.We define the centroid of a cluster as a dummy node associated with a three-value sequence, which is obtained as follows. The length of the sequence is the same as the sequences associated with the nodes in this cluster. For each bit in the string, we count the number of 0s (zeroCounts) and the number of 1s (oneCounts) in this bit position for all the sequences associated with the nodes in this cluster. If both numbers are equal to 0, this bit is set to 'x'. Otherwise this bit is set to 0 if zeroCounts ≥ oneCounts, and 1 if zeroCounts < oneCounts. The distance between a pair of clusters is then calculated as the product of the distance between the centroids of the pair of clusters and the total number of nodes in these two clusters. The above method for determining the centroid of a cluster and the distance between a pair of clusters has been developed to ensure that only a small number of bits need to be flipped to generate all the sequences in a cluster from the centroid sequence of the cluster.Figure6illustrates cluster analysis for sequence merg-

Figure 6 .

 6 Figure 6. Clustering for sequence merging.

Figure 7 .

 7 Figure 7. Proposed BIST architecture.

b 00 b 04 b 08 b 00 b 04 b 10 b 14 b 18 b 10 b 14 b 14 b 14 bFigure 8 .c 0 c 4 c 8 c 0 c 4 c 1 c 5 c 9 c 1 c 5 c 2 c 6 c x c 2 c 6 c 3 c 7 c x c 3 c 7 GroupFigure 9 .

 141414845679 Figure 8. Organization of the stored sequences.

  of faults left undetected by 10000 pseudo-random patterns N2: No. of test patterns generated by Atalanta for the undetected faults dex is completely shifted into the Bit counter B to determine the scan cycle in which the signal F lip indication D is set to 1. The bit index is divided into two parts. The lower-order log 2 m bits are used to determine which of the m bits in the signal F lip indication D are set to 1; the remaining bits are shifted into the Bit counter B to determine in which scan cycle to set the selected bit of F lip indication D to 1.

Table 1 . Results for ISCAS-89 benchmarks.

 1

						No. of
				No. of	No. of	remaining No. of	No. of
				groups	bits	faults after	test	bits
				for the	used to application vectors needed for
	No. of			sequence store the	of 10000	need encoding
	total			extraction extracted semi-random to be	test
	Circuits patterns N1 N2 procedure sequences patterns encoded patterns
	s5378 4563	51 38	8	136	3	3	85
	s9234 6475 735 317	8	719	81	38	1447
	s13207 9664 624 366	16	1021	43	30	610
	s15850 11336 667 248	16	962	2	2	36
	s35932 35110	0	-	-	-	-	-	-
	s38417 31015 2245 1006	32	1899	105	71	2301
	s38584 34797 448 271	32	1007	34	28	668
	N1: No.

Table 2 . Comparison of storage (in bits) required for various BIST methods.

 2

		Partial	Reseeding of	2-D
	Circuit	reseeding [7]	folding counter [5]	compression [10]	Proposed approach
	s5378	502	132	196	221
	s9234	5013	2310	3800	2166
	s13207	3008	247	1044	1631
	s15850	5204	2403	3360	998
	s38417	24513	6802	11214	4200
	s38584	2942	660	2891	1675

Table 3 .

 3 Comparison with[START_REF] Jas | Hybrid BIST based on weighted pseudo-random testing: a new test resource partitioning scheme[END_REF].

		Hybrid BIST [6]	Proposed approach
		On-chip	Encoded	On-chip	Encoded
		storage	data	storage	data
	Circuit	requirement (bits)	volume (bits)	requirement (bits)	volume (bits)
	s5378	N/A	N/A	88	13
	s9234	452	865	620	1178
	s13207	168	263	263	12
	s15850	436	1070	456	108
	s38417	2336	4680	1797	2053
	s38584	712	961	417	826

This research was supported in part by the National Science Foundation under grants CCR-9875324 and CCR-0204077.

Proceedings of the Design, Automation and Test in Europe Conference and Exhibition (DATE'05)

Proceedings of the Design, Automation and Test in Europe Conference and Exhibition (DATE'05) 1530-1591/05 $ 20.00 IEEE


INTRODUCTION

In the last decades, an increasing interest has been devoted to the development of numerical techniques such as stochastic finite elements (e.g. [START_REF] Ghanem | Stochastic finite elements: A spectral approach[END_REF][START_REF] Babuška | Solving elliptic boundary value problems with uncertain coefficients by the finite element method: the stochastic formulation[END_REF], in order to assess the impact of randomness on the response of mechanical systems. Uncertainties in material properties or loading are quite well mastered in the framework of these techniques. However, there are only a few available numerical strategies to deal with uncertainties in the geometry although it is of great interest in various applications. The effect of randomness in the geometry on the response of mechanical systems is often important. It is thus imperative to be able to correctly identify the randomness in geometry for obtaining reliable numerical predictions of the behavior of such systems.

The identification of a probabilistic model is a very critical point. It usually requires many samples and robust identification techniques of random variables or fields. Experimental campaigns or in-site measurements are often very expensive, a fact that drastically limits the number of available samples and thus the quality of identification. Among the different techniques that can be used to identify random shapes, shape recovery from simple images ("pictures" of the shape) has the following advantages: it is non-intrusive, it allows obtaining many samples at a very low cost and it can be relatively precise. This identification procedure leads to a characterization of the random geometry in a form that can be directly used for the numerical simulation of the physical problem. For the computation of the response, one generally has to solve partial differential equations defined on random domains. The number of numerical strategies proposed for this kind of problems is limited (e.g. [START_REF] Xiu | Numerical methods for differential equations in random domains[END_REF][START_REF] Canuto | A fictitious domain approach to the numerical solution of pdes in stochastic domains[END_REF]. Recently, the eXtended Stochastic Finite Element Method (X-SFEM) [START_REF] Nouy | X-SFEM, a computational technique based on X-FEM to deal with random shapes[END_REF][START_REF] Nouy | An eXtended Stochastic Finite Element Method for solving stochastic partial differential equations on random domains[END_REF] has been proposed. This approach relies on the implicit representation of complex geometries using random level-set functions and on the use of a Galerkin approximation at both stochastic and deterministic levels.

The methodology proposed in this paper is an efficient identification procedure of random geometry in a form suitable for numerical simulation within the X-SFEM. The method starts from a collection of images, representing different outcomes of the random shape to identify. The key-point is to represent the random geometry in an implicit manner using the level-set technique [START_REF] Sethian | Level-set methods and fast marching methods[END_REF]. This technique consists in representing the boundary of a shape with a level-set function, which is the signed distance function to the boundary. In our case, the random geometry will be characterized by a random level-set function, one outcome of which represents an outcome of the random boundary of the shape to be recovered. The problem of random geometry identification is thus equivalent to the identification of a random level-set function, which is a random field.

Identification of random geometry for stochastic finite element analysis G. Stefanou & A. Nouy

Research Institute in Civil Engineering and Mechanics (University of Nantes, Ecole Centrale Nantes, CNRS), Nantes, France ABSTRACT: In this paper, an efficient procedure is proposed for the identification of random geometry in a form suitable for numerical simulation within the recently proposed eXtended Stochastic Finite Element Method (X-SFEM). The method starts from a collection of images representing different outcomes of the random shape to identify. The key-point of the method is to represent the random geometry in an implicit manner using the level-set technique. In this context, the problem of random geometry identification is equivalent to the identification of a random level-set function, which is a random field. This random field is represented on a polynomial chaos basis and various efficient numerical strategies are proposed in order to identify the coefficients of its polynomial chaos decomposition. The performance of these strategies is evaluated on a "manufactured" random geometry problem and useful conclusions are provided with regard to the accuracy and efficiency of the proposed method. Some general techniques for the identification of random fields based on polynomial chaos decomposition have been proposed by [START_REF] Desceliers | Maximum likelihood estimation of stochastic chaos representations from experimental data[END_REF], [START_REF] Ghanem | On the construction and analysis of stochastic models: characterization and propagation of the errors associated with limited data[END_REF]. Here, this idea is adopted and the identification of a polynomial chaos representation of the random level-set is performed. More precisely, some efficient numerical strategies for the identification of the coefficients of the decomposition are proposed.

The paper is organized as follows: in section 2, we describe the principle of shape recovery from images with the level-set technique. In section 3, we introduce an empirical Karhunen-Loève decomposition of the samples, which allows representing the samples of the level-set on a reduced basis of deterministic modes and corresponding random variables. In sections 4 and 5, some methodologies are presented in order to identify a polynomial chaos representation of random variables. In section 4, we deal with the general case of mutually dependent random variables, while in section 5 a hypothesis of independence of random variables is introduced. Section 6 will illustrate the methodologies on a "manufactured" problem and finally some conclusions will be provided in section 7.

THE LEVEL-SET TECHNIQUE FOR SHAPE RECOVERY

The problem of shape recovery from an image is now a well mastered problem within the context of level-set techniques [START_REF] Sethian | Level-set methods and fast marching methods[END_REF]. Here, we briefly recall the basis of this method. We suppose that we have a contrasted image, defined by a mapping :

I D ∈ → x
ℝ , whose value I(x) represents the grayness intensity at location D ∈ x (where 2 D ⊂ ℝ or even 3 D ⊂ ℝ ). The aim is to detect the boundary of the underlying shape. This boundary is in fact located in the region where the intensity has the highest gradients. The aim of shape recovery consists in building a level-set function ( ) φ x whose iso-zero is located in this region with high intensity gradients. The basic idea consists in propagating a front, represented by this iso-zero of a time-dependent level-set ( , ) t φ x , which will "lock" on the desired boundary.

The equation of motion of a level-set ( , )

t φ x is a Hamilton-Jacobi partial differential equation of the form: ( , ) ( , ) ( , ) 0 t t F t t φ φ ∂ + ∇ = x x x (1) 0 ( , 0) ( ) φ φ = x x
where F is the speed of the front in the outward normal direction (from negative to positive values of φ ). In order to make the iso-zero lock in high inten- sity gradients zones, the speed has to vanish in these zones. A classical choice for F [START_REF] Sethian | Level-set methods and fast marching methods[END_REF] is the following:

1 (1 ) 1 ( ) F c G I σ εκ = - + ∇ * (2)
where κ is the curvature of the front, ε>0 a small parameter, I the mapping of grayness intensity and G σ a Gaussian smoothing filter with characteristic width σ.

( ) G I σ ∇ *
represents the gradient of the image convolved with the filter. The curvature term is a classical regularization term leading to a smooth front. Parameter c allows imposing an arbitrary small value of the speed in high intensity gradients zones.

A basic choice for the initial level-set 0 φ consists in a small circular front at the interior of the boundary to be recovered. Many algorithms have been proposed in order to solve the equation of motion (1) (see [START_REF] Sethian | Level-set methods and fast marching methods[END_REF]. After discretization and resolution, this leads to a discretized level-set N φ ∈ ℝ .

REDUCTION OF INFORMATION THROUGH KARHUNEN-LOEVE DECOMPOSITION

In our problem, the random shape to identify is a discretized random level-set represented by a random vector :

N φ Θ → ℝ defined on an abstract prob- ability space ( , , ) P Θ Β . The sample discretized level- sets ( ) 1 { } k Q k φ
= are obtained by applying Q times the shape recovery technique of section 2 to a set of sample images ( )

1 { } k Q k I = .
The probabilistic identification of the random discretized level-set φ from samples ( ) { } k φ is infeasible in practice due to the high dimensionality of the underlying probabilistic inverse problem. A reduction of information is thus unavoidable. This reduction of information consists in expressing the levelset in terms of a small set of random variables 1 ( ,... ) n X X := X. The problem of identification of the random level-set is then transformed into the problem of identification of a smaller random vector X(θ) from a collection of samples

( ) 1 { } k Q k = X .
In the general case where the shape is not known a priori, a reduction of information through discrete empirical Karhunen-Loève (K-L) decomposition of level-set samples is used.

If φ µ is the empirical mean of level-set samples

( ) k φ and ( ) ( ) k k φ φ φ µ = - ɶ
the centered samples, the unbiased empirical covariance matrix of samples writes:

( ) ( ) 1 1 1 Q k k T k Q φ φ φ = = -∑ C ɶ ɶ (3)
We denote by

) ( ) ( ) 1 ˆm k k k i i i i s X φ φ φ µ = ≈ = + ∑ U , ( ) ( ) 1 k T k i i i X s φ = U (4)
where the

( ) 1 { } k m i i X = appear as the components of ( ) k φ ɶ on the reduced basis of modes 1 {( )} m i i i s = U
associated with the m larger eigenvalues.

It is easy to show that the approximation error is:

2 ( ) ( ) 1 1 1 1 Q N k k i k i m s Q φ φ = = + - = -∑ ∑ (5)
where ⋅ denotes the classical 2-norm in N ℝ . The number of modes to be kept can then be chosen as the smallest integer m such that 2 1 1

N N i KL i i m i s s ε = + = ≤ ∑ ∑ , (6)
where KL ε is a tolerance which is fixed a priori and allows controlling the reduction error. Finally, samples of the level-set are approximated by:

( ) ( ) 1 m k k i i i X φ φ µ = ≈ + ∑ U ɶ , i i i s = U U ɶ (7)
and the X i to be identified are centered normalized uncorrelated random variables.

POLYNOMIAL CHAOS REPRESENTATION OF RANDOM VARIABLES

The problem is now to identify a random vector :

m Θ → X
ℝ defined on an abstract probability space ( , , ) P Θ Β , based on a set of independent samples

( ) 1 { } k Q k = X .
In this section, no further assumption will be made and a random vector with arbitrary probability law will be identified (a simplifying assumption of independence between the random variables will be introduced in section 5). For this purpose and following [START_REF] Desceliers | Maximum likelihood estimation of stochastic chaos representations from experimental data[END_REF], we use a polynomial chaos (PC) representation of X, identified with a maximum likelihood principle. This leads to the solution of an optimization problem on a Stiefel manifold. A numerical strategy is proposed for solving this optimization problem.

Polynomial chaos decomposition

We consider that X is a second order random vector and that a non-linear mapping :

m m g → ℝ ℝ exists such that ( ) = X g ξ , where 1 ( ,..., ) m ξ ξ = ξ
is an mdimensional vector of independent random variables with known probability law P ξ and support m ⊂ Ξ ℝ . We denote by ( , , )

B P Ξ ξ Ξ the m-dimensional probability space defined by random vector ξ. Then, random vector X admits a generalized chaos representation [START_REF] Xiu | The Wiener-Askey polynomial chaos for stochastic differential equations[END_REF][START_REF] Soize | Physical systems with random uncertainties: chaos representations with arbitrary probability measure[END_REF]):

( ) ( ( )) m H α α α θ θ ∈ℑ = ∑ X X ξ (8)
where { } m H α α∈ℑ is a Hilbertian basis of 2 ( , ) L dP ξ Ξ endowed with its natural inner product. Coefficients X α of the chaos decomposition are then simply defined as the L 2 -projection of g on the :

m H α α ∀ ∈ ℑ , ( )
2 ( , )
, ( ) ( )

L P H E H α α α = = = ξ Ξ X g g ξ ξ ( ) ( ) ( ) H dP α Ξ = ∫ ξ g y y y (9)
where E denotes the mathematical expectation. The random vector X is finally approximated by truncating the generalized PC expansion to a degree p in Equation 8:

, ( ) ( ( )) m p H α α α θ θ ∈ℑ ≈ ∑ X X ξ (10) where ,
{ ; } m m p m p α α ℑ = ∈ ℑ = ≤ ℕ .

Maximum likelihood estimation of chaos coefficients

As in [START_REF] Desceliers | Maximum likelihood estimation of stochastic chaos representations from experimental data[END_REF], we impose the second order moments of samples to be conserved after the probabilistic identification. Let us denote

Maximum likelihood principle

The likelihood function, for a given matrix of coefficients A, is defined by

( ) 1 ( ) ( ; ) Q k k L p = = ∏ X A X A , (12)
where X (k) are the samples and ( ; ) p ⋅ X A is the joint probability density function (pdf) of the chaos decomposition of X. Due to the prohibitive cost of the estimation of joint pdfs and also due to the small number of samples available in practice, it was proposed by [START_REF] Desceliers | Maximum likelihood estimation of stochastic chaos representations from experimental data[END_REF] to use a pseudo likelihood function:

( ) 1 1 ( ) ( ; ) i Q m k X i k i L p X = = = ∏∏ A A (13)
whose estimation requires only the evaluation of marginal pdfs. In practice, for numerical reasons, the opposite of the pseudo log-likelihood function is rather used:

( ) 1 1 ( ) log( ( )) log( ( ; )) i Q m k X i k i f L p X = = = - = -∑∑ A A A . ( 14
)
Finally, matrix A is defined by the following optimization problem on a compact Stiefel manifold:

( , ) arg min ( )

S P m f ∈ = A A A (15)

Solution of the optimization problem

The solution of optimization problem ( 15) is a relatively hard task due to the nature of function f and the possibly high dimensionality. In particular, f may present local minima probably leading to a bad probabilistic representation of X. Thus a global optimization procedure must be provided. In this work, a global random search algorithm is used in order to provide several initializations for a classical descent algorithm. Samples of the random search algorithm are generated with respect to the uniform probability measure on the Stiefel manifold (see e.g. [START_REF] Fang | Some methods for generating both an nt-net and the uniform distribution on a Stiefel manifold and their applications[END_REF]. In order to use a wider class of descent algorithms and to reduce the dimension of the optimization problem, the problem can be reformulated as an unconstrained optimization problem. This is possible by introducing a minimal parametrization of the compact Stiefel manifold ( , ) S P m (see [START_REF] Stefanou | Identification of random shapes from images through polynomial chaos expansion of random level set functions[END_REF]. In many cases, unconstrained optimization algorithms exhibit relatively good convergence properties for problems where constrained optimization algorithms do not converge.

POLYNOMIAL CHAOS REPRESENTATION OF MUTUALLY INDEPENDENT RANDOM VARIABLES

The resolution of problem ( 15) resulting from a general maximum likelihood approach may be infeasible in practice due to the possible large number of parameters to be identified. In this section, we will assume that random variables X i are mutually independent. This hypothesis, which will be validated in the numerical example, seems reasonable in practice due to the small number of available samples. Under this assumption, each random variable can be separately decomposed on a one-dimensional polynomial chaos basis of degree p:

, , 0 0 ( ) ( ( )) ( ( )) p i i i i i X X h X h α α α α α α θ ξ θ ξ θ ∞ = = = ≈ ∑ ∑ (16)
where ξ i are independent identically distributed random variables with support Ξ ⊂ ℝ and probability law P ξ . { } h α α∈ℕ form an orthonormal polynomial basis of 2 ( , ) L dP ξ Ξ

. Two possible strategies for the identification of the one-dimensional chaos decompositions are proposed: the maximum likelihood estimation described in the previous section (greatly simplified when applied to independent random variables) and a projection method using empirical cumulative distribution functions of samples as an approximation of mapping g appearing in Equation 9 [START_REF] Stefanou | Identification of random shapes from images through polynomial chaos expansion of random level set functions[END_REF]. The latter method leads to a very fast computation of the decomposition.

NUMERICAL EXAMPLE: RANDOM ROUGH CIRCLE RECOVERED FROM SAMPLE IMAGES

In this example, the samples of the random shape to identify are obtained using the shape recovery technique based on the level-set method described in section 2 of the paper. The starting point for the identification procedure is a set of Q=100 images, which are artificially generated from a random levelset function analytically defined in a square domain [0,1]× [0,1] as follows:

( , ) ( , ) R a

φ θ θ = -- x x c (17)
The iso-zero of level-set φ is a random "rough" cir- cle C(θ) shown in Figure 1. c is the "center" of the circle and ( , ) R a θ is a random field indexed by the polar angle a of the polar coordinate system centered at c, defined as:

1 2 1 ( , ) 0.2 0.03 ( ) 0.015[ ( ) cos( ) R a Y Y k a θ θ θ = + + + 3 1 4 2 5 2 ( )sin( ) ( ) cos( ) ( ) sin( )] Y k a Y k a Y k a θ θ θ + + (18)
where k 1 , k 2 are deterministic constants and

1 5 ( ( ),..., ( )) Y Y θ θ are independent identically distrib- uted uniform random variables: ( 3, 3) i Y U ∈ - , i=1,…,5.
With an appropriate choice of the speed of the front (see Equation 2) and by applying 100 times the shape recovery technique to the initial set of sample images, the corresponding set of 100 discretized level-sets to identify is obtained. It is worth noting that, since the recovery procedure involves calculation of image gradients whose accuracy depends on the mesh size, a sufficiently fine mesh must be used for the recovery in order to be able to capture details of shape features. In our case, a 100× 100 mesh has been used to this purpose. Figure 2 illustrates the filtered gradients of two sample images I (k) and the corresponding iso-zero of the recovered level-sets. Legendre and Hermite chaos are used for the representation of the random variables X resulting from K-L decomposition of the recovered level-sets. The pseudo log-likelihood function, the pdf of the initial (resulting from K-L expansion) and identified random variables as well as the error in the probability P in (x) to be inside the hole are used as error criteria in the comparison. The probability P in (x) is defined by: ( ) ( ( , ) 0)

in P P φ θ = < x x (19)
and can be evaluated from samples as follows:

{ } ( ) 1 ( ) {1... }; ( ) 0 k in P Card k Q Q φ ≈ ∈ < x x . ( 20
)
Using a tolerance KL ε equal to 0.05 in Equation 6leads to a K-L decomposition of the sample levelsets in 8 terms and thus 8 random variables have to be identified in this case. The K-L modes corresponding to the 8 retained terms are shown in Figure 3.

U 1 U 2 U 3 U 4 U 5 U 6 U 7 U 8
Figure 3. The K-L modes corresponding to the 8 retained terms of the sample level-sets.

In Figure 4, the values of the pseudo log-likelihood calculated for the initial and identified random variables are plotted as a function of the order p of PC decomposition using the three identification alternatives described in sections 4-5 and the two types of PC. Values of the pseudo log-likelihood up to p=3 have been calculated in the case of the maximum likelihood estimation without independence hypothesis (denoted as ML(D)) due to the very large computational cost required for higher order PC expansions. For Hermite PC, the results obtained with the maximum likelihood estimation with independence hypothesis (denoted as ML(I)) and the projection method using empirical cumulative distribution functions (denoted as EP(I)) practically coincide. For low order Legendre PC, the ML(D) technique gives better results in likelihood but the results of the two other methods substantially improve as p grows up. This fact implies that the independence assumption is suitable to this problem and has not been strongly affected by the recovery procedure. In the Legendre PC case, the convergence is somewhat slow perhaps due to the complexity of the problem. The marginal pdfs of the first 5 random variables shown in Figure 5 are close to the uniform confirming that the recovery procedure was sufficiently accurate. Similar results are obtained with the EP(I) technique and Legendre chaos of order p=1 and p=3 for the error in the probability to be inside the hole, presented in Figure 6. Before closing this section, it can be stated that the two techniques assuming independence are very competitive in terms of accuracy compared to the ML(D) approach, while the projection method EP(I) is the approach having the smallest computational cost for the identification procedure.

-2 0 2 0 1 2 X 1 x PDF(x) -2 0 2 0 1 2 X 2 x PDF(x) -2 0 2 0 1 2 X 3 x PDF(x) -2 0 2 0 1 2 X 4 x PDF(x)

CONCLUSIONS

In this paper, an efficient identification procedure of random geometry has been proposed in a form suitable for numerical simulation within the eXtended Stochastic Finite Element Method (X-SFEM). The method starts from a collection of images representing different outcomes of the random shape to identify. The key-point of the method is to represent the random geometry in an implicit manner using the level-set technique. In this context, the problem of random geometry identification is equivalent to the identification of a random level-set function, which is a random field. This random field is first decomposed using an empirical Karhunen-Loève expansion, which allows representing the samples of the level-set on a reduced basis of deterministic modes. The problem is thus transformed into the probabilistic identification of a few random variables, which are the components of the random level-set on this reduced basis of modes. The random variables are represented on a polynomial chaos basis and three efficient numerical strategies are used in order to identify the coefficients of their polynomial chaos decomposition. The performance of these strategies has been evaluated on a "manufactured" random geometry problem using various error criteria. It can be concluded that the two techniques assuming independence (ML(I) and EP(I)) are very competitive in terms of accuracy compared to the ML(D) approach, while the projection method EP(I) is the approach having the smallest computational cost for the identification procedure. This conclusion is in accordance with the assumption of independence which is very often used in practice due to the small collection of images that is usually available.

U

  the eigenvalues and eigenvectors of φ C . Vectors U i form an orthonormal basis of N ℝ . Then level-set samples can be approximated in the following way:

(

  whose rows are the remaining chaos coefficients of X. The conservation of covariance matrix imposes the following constraint on A: the coefficients is then equivalent to the identification of an orthogonal matrix P m × ∈ A ℝ . The set of orthogonal matrices in P m × ℝ is the compact Stiefel manifold denoted ( , ) S P m .

Figure 1 .

 1 Figure 1. Schematic representation of a "rough" circle.

Figure 2 .

 2 Figure 2. Gradient intensity of two filtered sample images and corresponding iso-zero of the recovered level-sets.

  Pseudo log-likelihood values of the initial and identified random variables as a function of the order p of PC decomposition: (a) Legendre chaos, (b) Hermite chaos.

Figure 5 .

 5 Figure 5. Marginal pdf of the first 5 initial and identified random variables (Legendre chaos, p=1, EP(I)).

  Error in the probability to be inside the hole (EP(I), Legendre chaos) with (a) p=1 and (b) p=3.

ACKNOWLEDGEMENT

This work has been supported by the French National Research Agency (Grant ANR-06-JCJC-0064). This support is gratefully acknowledged.


M Zidi
email: zidi@univ-paris.12.fr

Large Shearing of a Prestressed Tube

This study is devoted to a prestressed and hyperelastic tube representing a vascular graft subjected to combined deformations. The analysis is carried out for a neo-Hookean response augmented with unidirectional reinforcing that is characterized by a single additional constitutive parameter for strength of reinforcement. It is shown that the stress gradients can be reduced in presence of prestress.

Introduction

Mechanical properties are of major importance when selecting a material for the fabrication of small vascular prostheses. The operation and the handing of prostheses vessel by surgeons, on the one part, the design of such grafts, on the other, induce specific loading and particularly boundary or initial conditions. Consequently, the interest in developing a theoretical model to describe the behavior of the prostheses vessel is proved ͓͑1͔͒. In this paper, we consider a thick-walled prestressed tube, hyperelastic, transversely isotropic, and incompressible assimilated to a vessel graft. We give an exact solution of the stress distributions when the tube is subjected to the simultaneous extension, inflation, torsion, azimuthal, and telescopic shears ͓͑2-10͔͒. The first theoretical re-sults, in the case of a silicone tube, indicate that the increase of prestress minimizes the stress gradients due to the effects of the shear.

Model Formulation

Consider a nonlinearly elastic opened tube defined by the angle 0 ͑Fig. 1͒. Let us suppose that the tube undergoes two successive deformations; first, including the closure of the tube which induced residual strains ͓͑11͔͒ and second, including inflation, extension, torsion, azimuthal and telescopic shears. The mapping is described by

rϭr͑R ͒ ϭ ͩ 0 ͪ ϩ␣Zϩ⌰͑r ͒ zϭ␣Zϩ⌬͑r ͒ (1)
where (R,,Z) and (r,,z) are, respectively, the reference and the deformed positions of a material particle in a cylindrical system. is a twist angle per unloaded length, ␣ and are stretch ratios ͑respectively, for the first and the second deformation͒, ⌰ is an angle which defined the azimuthal shear, and ⌬ is an axial displacement which defined the telescopic shear.

It follows from ͑1͒ that the physical components of the deformation gradient F has the following representation in a cylindrical system:

Fϭ ͫ r ˙͑R ͒ 0 0 r͑R ͒⌰ ˙͑r ͒r ˙͑R ͒ r͑R ͒ R 0 r␣ ⌬ ˙͑r ͒r ˙͑R ͒ 0 ␣ ͬ (2)
where the dot denotes the differentiation with respect to the argument.

Incompressibility then requires that Jϵdet Fϭ1, which upon integration yields

r 2 ϭr i 2 ϩ 0 ␣ ͑ R 2 ϪR i 2 ͒ (3)
where R i and r i are, respectively, the inner surfaces of the tube in the free and in the loaded configurations ͑R e and r e are the outer surfaces͒.

The strain energy density per unit undeformed volume for an elastic material, which is locally and transversely isotropic about the t(R) direction, is given by WϭW͑I 1 ,I 2 ,I 3 ,I 4 ,I 5 ͒ (4)

where

I 1 ϭTrC, I 2 ϭ 1 2 ͓͑TrC͒ 2 ϪTrC 2 ͔, I 3 ϭ1, I 4 ϭtCt, I 5 ϭtC 2 t (5)
are the principal invariants of CϭF ¯F which is the right Cauchy-Green deformation tensor ͑F ¯is the transpose of F͒.

The corresponding response equation for the Cauchy stress for transversely isotropic incompressible is ͑see ͓12͔͒

ϭϪp1ϩ2͓W 1 BϪW 2 B Ϫ1 ϩI 4 W 4 T T ϩI 4 W 5 ͑ T B"TϩT"B T͔͒ (6)
where BϭFF ¯is the left Cauchy-Green tensor, 1 the unit tensor, and p the unknown hydrostatic pressure associated with the incompressibility constraint, W i ϭ(ץW/ץI i ) (iϭ1,2,4,5) and Tϭ(1/ͱI 4 )Ft.

From ͑6͒, the equilibrium equations in the absence of body forces are reduced to

d rr dr ϩ rr Ϫ r ϭ0 ( 7 a) d r dr ϩ 2 r r ϭ0 ( 7 b) d rz dr ϩ rz r ϭ0. ( 7c
)
Suppose that ⌰ and ⌬ satisfy the following boundary conditions: ͑a͒ ⌰ϭ⌰ i , ⌬ϭ⌬ i in rϭr i and ͑b͒ ⌰ϭ⌰ e , ⌬ϭ⌬ e in r ϭr e . Then, a simple computation by integrating ͑7b͒ and ͑7c͒ gives the expression of ⌰ and ⌬.

Integrating ͑7a͒, given the boundary conditions that rr (r i )ϭ Ϫp i and rr (r e )ϭ0, and taking t(R)ϭt (R)e ϩt Z (R)e Z and using ͑3͒ yields the pressure field p:

p͑r ͒ϭ p i ϩ2W 1 ͩ R 0 r␣ ͪ 2 Ϫ2W 2 f ͑ r ͒ϩ ͵ r i r rr Ϫ s ds (8a)
where The expressions of ⌰, ⌬, and p determine all the components of the Cauchy stress tensor .

f ͑ r ͒ϭ⌬ ˙2͑r ͒ ͫ 1 ͑ ␣͒ 2 ϩ ͩ R 0 ͪ 2 ͬ ϩ⌰ ˙2͑r ͒ ͩ R 0 ͪ 2 Ϫ2 ⌰ ˙͑r ͒⌬ ˙͑r ͒ 0 2 ␣ ϩ ͩ r␣ R 0 ͪ 2 . ( 8b
)

Results

To illustrate the response of the proposed model, we use the extended Mooney Rivlin strain energy function which represents the behavior of a prosthesis ͓͑13͔͒ constituted of a silicone matrix and textile fibers,

WϭW͑I 1 ,I 4 ͒ϭ 2 ͑ I 1 Ϫ3 ͒ϩ E f 8 ͑ I 4 Ϫ1 ͒ 2 , (9)
where is the shear modulus of the isotropic matrix at infinitesimal deformations and E f is the elastic modulus of the fibers.

The local tangent vector of the fibers is chosen here as t(R) ϭcos ␥(R)e ϩsin ␥(R)e Z that represent a helical distribution of fibers ͓͑1͔͒.

From Eqs. ͑7b͒, ͑7c͒ and using ͑3͒ it easily follows that the expressions of ⌰ and ⌬ are where kϭ␣/R i 2 0 . As an illustrative result, we focus our attention only when the tube is submitted to azimuthal shear strain. Figure 2 shows the distribution of circumferential stresses generated by applied exter-nal azimuthal strain at a given pressure when taking into account the effects of such residual stresses. We show clearly that a decrease in 0 angle helps to distribute stresses in the loaded state when the shear is important. This result does not change qualitatively when varying the pressure p i .

⌰͑r ͒ϭ͑ ⌰ e Ϫ⌰ i ͒ log ͫ r r i ͱ1ϩk͑r 2 Ϫr i 2 ͒ ͬ log ͫ r e
Furthermore, the particular effects of the presence of fibers have been examined with a linear distribution of fiber orientation within the data range ␥(R i )ϭϪ40 deg and ␥(R e )ϭ40 deg. As illustrated in Fig. 3, it is shown here that the effects of the azimuthal shear upon the distribution of the circumferential stresses within the wall become significant. When the tube is prestressed, the stresses are also distributed. Clearly these results will be able to help the design and fabrication of a small vascular prosthesis ͓͑1͔͒.

Fig. 1

 1 Fig. 1 Cross section of the tube in the stress-free "a…, unloaded "b…, and loaded configuration "c…


Mohammad Tavallaeinejad
email: mohammad.tavallaeinejad@mail.mcgill.ca

Mathias Legrand
email: mathias.legrand@mcgill.ca

Michael P Païdoussis

Nonlinear dynamics of slender inverted flags in axial flow

A nonlinear fluid-elastic continuum model of a slender cantilevered plate subjected to axial flow directed from the free end to the clamped one, also known as the inverted flag problem, is proposed. A nonlinear unsteady slender wing theory is employed to express the fluid-related forces

acting on the plate. It is based on potential flow theory for inviscid fluids. The Euler-Bernoulli beam theory with exact kinematics and inextensibility is employed to derive the nonlinear partial integro-differential equation governing the dynamics of the plate. Discretization in space is carried out via a conventional Galerkin scheme using the linear modeshapes of the cantilevered beam. A Newton solver is implemented to obtain equilibrium solutions and integration in time is conducted using Gear's backward differentiation formula. A bifurcation diagram in terms of flow velocity is constructed in order to gain insight on the stability and post-critical behaviour of the system. It is shown that the undeflected static equilibrium is stable prior to a supercritical Hopf bifurcation giving rise to a flapping motion around the undeflected static equilibrium. By further increasing the flow velocity, the flag displays flapping motions around deflected static equilibria and change to fully-deflected static states at even higher flow velocities. These predictions are in excellent agreement with existing experimental data available in the literature. Second space derivative

INTRODUCTION

An inverted flag is here understood as a cantilevered plate subjected to axial flow and clamped along its downstream edge in the flow. Applications, particularly in energy harvesting systems, have motivated recent studies in the literature. The flapping flag instability was first explored experimentally by Kim et al. [START_REF] Kim | Flapping dynamics of an inverted flag[END_REF] and later by Sader et al. [START_REF] Sader | Large-amplitude flapping of an inverted flag in a uniform steady flow-a vortex-induced vibration[END_REF] for various aspect ratios. Experimental observations have revealed that the flag aspect ratio can dramatically alter the dynamics. In particular, slender flags exhibit a buckling behaviour via a divergence instability. Wide flags, on the other hand, undergo a flapping motion prior to fully bending on one side via buckling. Moreover, the undeflected static equilibrium of inverted flags with smaller aspect ratios loses stability at higher critical flow velocities. These behaviours have been explored numerically in recent studies via three-dimensional computational simulations [START_REF] Tang | Dynamics of an inverted flexible plate in a uniform flow[END_REF][START_REF] Gurugubelli | On the mechanism of large amplitude flapping of inverted foil in a uniform flow[END_REF]. In a remarkable work, Sader et al. [START_REF] Sader | Stability of slender inverted flags and rods in uniform steady flow[END_REF] performed an experimental and theoretical investigation and provided an analytical framework to explore the stability and dynamics of inverted flags in the asymptotic limit of zero aspect ratio, employing a steady-flow model to represent the aerodynamic forces. They found the emergence of a bistable region wherein both the undeflected and the deflected static regimes are stable, which cause inverted flag to display intermittent dynamics; any contribution of unsteadiness originating from the flow was neglected in their analytical model.

The present paper aims at deriving an analytical model accounting for flow unsteadiness. The resulting nonlinear dynamics is then investigated in the limit case of slender flexible plates. The sensitivity of the dynamics to the aspect ratio is also explored. The results are summarized in the form of bifurcation diagrams to identify the transition between various regimes.

MATHEMATICAL FORMULATION

In this section, the kinetic and potential energies of the cantilevered elastic beam shown in Fig. 1 are derived on the basis of the conventional nonlinear Euler-Bernoulli beam theory. Next, the virtual work associated with fluidrelated forces is formulated. Finally, the equation of motion is obtained via Hamilton's principle.

Potential and Kinetic Energies

Figure 2 shows a generic point at a distance z from the mid-plane on the cross-section of the flag; z and x stand for the transverse and axial directions, respectively, and t is time. Denoting by w the transverse motion, and by u, the longitudinal motion of the plate, the inextensibility assumption reads .

1 C u x / 2 C w 2 x D 1.
The potential and kinetic energies are expressed as

T .t / D A 2 Z L 0 . P w 2 C P u 2 / dx C I 2 Z L 0 P 2 dx; (2) V .t / D D 2 Z L 0 2 x dx; ( 3
)
where D D EH h 3 =OE12.1 2 / and denotes the Poisson's ratio.

Slender-body Nonlinear Model

The external aerodynamic force acting along the inverted flag can be decomposed into the reactive and the resistive forces, F N and F R , respectively. The resistive component models the viscous effects such as separation and drag, while the reactive force originates from the acceleration of fluid induced by the motion of the flag. The extension of large-amplitude elongated body theory [START_REF] Lighthill | Large-amplitude elongatedbody theory of fish locomotion[END_REF][START_REF] Buchak | The clapping book: Wind-driven oscillations in a stack of elastic sheets[END_REF] combined with large displacements (that is, large ) yields

F N D 4 f H 2  P V n C OEV n V x 1 2 V 2 n x à ; (4)
where V n D . P u C U / sin C P w cos and V D . P u C U / cos C P w sin denote the normal and tangential components of the relative velocity of the element with respect to the incident flow, respectively. A closed form of the resistive force can be obtained from semi-empirical expressions or analogical work based on potential-flow theory. In this paper, the resistive force is modelled on the basis of the Polhamus leading-edge suction analogy [START_REF] Polhamus | A concept of the vortex lift of sharp-edge delta wings based on a leading-edgesuction analogy[END_REF] in the form

F R D f H U 2 2 .K p sin ˛cos ˛C K v sin 2 ˛/; (5)
with K p and K v being determined employing Bollay's nonlinear wing model [START_REF] Bollay | A non-linear wing theory and its application to rectangular wings of small aspect ratio[END_REF]. As shown in Fig. 2, the instantaneous angle of attack can be expressed as ˛D  [START_REF] Tavallaeinejad | Nonlinear response of inverted flags subjected to a steady flow[END_REF] where the relation tan  D P w=.U C P u/ holds. The virtual work done by the fluid-related forces normal to the plate is

ıW f D Z L 0 .F N F R /.ıw cos ıu sin /dx: (6)

Equation of Motion

Substituting Eqs. ( 2), ( 3) and ( 6) into the extended Hamilton's principle and using the following parameters: 7) is discretized via a Galerkin expansion .x; t / D P M rD1 r .x/q r .t /, where r .x/ denotes the rth eigenmode of the cantilevered beam and q r .t /, its corresponding generalized coordinate. In this study, six modes are employed. This is sufficient to obtain converged results. The resultant set of nonlinear ordinary differential equations (ODEs) is solved using Gear's backward differentiation formula yielding the time histories of the amplitude of oscillation. Static equilibria 0 .x/ (together with u 0 .x/ and w 0 .x/) are retrieved via a nonlinear Newton procedure solving the time-independent discretized governing equations.

x D x=L, t D t = , … D f SL 2 U 2 =D, D L 2 p A=D, D f L=. h/, ˇD I =.AL 2 /, k a D A=4, k p D K p =2, k v D K v =2, c d D c d L 4 =.D /,

RESULTS AND DISCUSSION

The nonlinear dynamics of the system is now examined by means of bifurcation diagrams. Phase portraits and (flapping) flag shapes are also provided. The parameters used in the analysis are L D 30 cm, h D 1 mm, D 1200 kg m 3 , f D 1:2 kg m 3 , D D 2454 N cm 2 and Á D 0:008. The damping coefficient c d is replaced by 2 !, where is the modal damping ratio and ! is the first dimensionless linear natural frequency. The modal damping ratio, , is set to 0:05 for the results obtained in this section.

Two configurations are considered: H D 3 cm and H D 10:5 cm corresponding to A D 0:1 and A D 0:35, respectively. Figures 3 and4 depict the bifurcation diagrams for inverted flags with A D 0:1 and A D 0:35, respectively, as a function of the dimensionless flow velocity …. Bold and dashed lines represent stable and unstable branches of the solution, respectively.

As seen in Fig. 3, no flapping motion is observed for A D 0:1. When the flow velocity is increased, the inverted flag deflects abruptly from its stable trivial equilibrium (no deformation) to a stable deflected equilibrium via a subcritical pitchfork bifurcation at … D 28:9. Also, the inverted flag may undergo a sudden large-amplitude deformation at lower flow velocities due to the emergence of a stable deflected equilibrium initiating from a saddlenode bifurcation at … D 8:6. This behaviour agrees with other existing analytical findings [START_REF] Tavallaeinejad | Nonlinear response of inverted flags subjected to a steady flow[END_REF][START_REF] Sader | Stability of slender inverted flags and rods in uniform steady flow[END_REF].

Figure 4 shows the bifurcational behaviour of the second case, that is a larger aspect ratio A D 0:35. Unlike the first case, the inverted flag undergoes a large-amplitude flapping motion via a supercritical Hopf bifurcation at … D 2:05; the amplitude of flapping increases with flow velocity. Remarkably, increasing the flow velocity causes a transition from a flapping regime around the undeflected equilibria to a deformed-flapping regime which the oscillation is around deflected equilibria. The latter motion dies out by further increasing the flow velocity. The inverted flag then undergoes a fully-defected regime via a secondary Hopf bifurcation at … D 9:95. It should be noted that the response of the system is symmetrical with respect to the horizontal axis and the inverted flag may flap around deflected equilibria (buckled position) on either side. Flapping around one side is plotted only, so as to properly distinguish flapping around the undeflected static equilibria and flapping around the deflected equilibria. Moreover, the fixed points of Eq. ( 7) are drawn as continuous (stable) and dashed (unstable) lines displaying a subcrtitical pitchfork bifurcation. This signifies that in case of introducing additional damping to the system, by touching the inverted flag using an external pole for instance [START_REF] Sader | Large-amplitude flapping of an inverted flag in a uniform steady flow-a vortex-induced vibration[END_REF], the static deflected solution turns into a stable equilibrium. As such, the inverted flag resides at this new induced equilibrium. This result is in line with the experimental treatment of Sader et al. [START_REF] Sader | Large-amplitude flapping of an inverted flag in a uniform steady flow-a vortex-induced vibration[END_REF]. From a comparison of the dynamics for A D 0:1 and for A D 0:35, it can be concluded that reducing the aspect ratio decreases the contribution of the reactive force, thus causing the divergence instability to take place at high flow veloci- ties. For sufficiently large aspect ratios, on the other hand, the effect of the reactive component of the force becomes dominant and the inverted flag loses stability via a Hopf bifurcation leading to a periodic limit cycle. These findings are in excellent qualitative agreement with experimental observations reported in the literature.

Phase portraits are shown in Fig. 5 demonstrating the static equilibrium at … D 2:0, the periodic limit cycle around undeflected static equilibrium at … D 8:0, the periodic limit cycle around a buckled position at … D 9:5, and finally, the static equilibrium at … D 20:0. Snapshots of the configurations corresponding to the foregoing flow velocities are plotted in Fig. 6.

CONCLUSION

The nonlinear response of an inverted flag is investigated in this paper. An analytical model is developed taking into account the geometric nonlinearities and the flag's inertia. Fluid-related forces are treated independently employing the extension of large-amplitude elongated body theory. The nonlinear partial integro-differential equation is solved numerically using Gear's backward differentiation formula. Equilibrium solutions are obtained via implementing a nonlinear Newton solver. The results are presented as bifurcation diagrams, phase portraits, and flapping shapes to explore various static and flapping regimes.

The various responses indicate that the change in the flag's aspect ratio substantially alters the dynamics. Although inverted flags of small aspect ratio do not flapping motion for varying flow velocities, inverted flags of sufficiently large aspect ratios undergo a flapping motion. More specifically, flags of very small aspect ratio lose stability via a subcritical pitchfork bifurcation giving rise to the birth of multiple equilibrium states. In this case, any contribution of the unsteady terms of the flow can be neglected. In contrast, the time-dependent nature of the flow is expected to increase in strength with increasing aspect ratio. This results in considerably richer and more interesting bifurcational behaviour of the system. In particular, the inverted flag exhibits various dynamical regimes: (i) it first loses stability via a supercritical Hopf bifurcation in the first mode, the amplitude of which increases with flow; (ii) this transforms to a deformed-flapping motion; and (iii) finally, the system resides at the fully deflected state.

FIGURE 1 :

 1 FIGURE 1: Inverted flag, idealized as a thin Euler-Bernoulli beam, in axial flow: free at the upstream and clamped at the downstream of the flow.

FIGURE 3 :

 3 FIGURE 3: Bifurcation diagrams for A D 0:1; (a) the tip slope and (b) the transverse/longitudinal displacement versus …. Stable [ ] and unstable [ ] static equilibrium.

1 …FIGURE 4 :

 14 FIGURE 4: Bifurcation diagrams for A D 0:35. (a) minimum/maximum values of the tip slope of the inverted flag. Stable static solution (undeflected stable equilibria) [ ], unstable static solution [ ], stable periodic solution (oscillation around undeflected equilibria) [ ], stable periodic solution (oscillations around deflected equilibria) [ ] and stable static solution (deflected static equilibria) [ ]. (b) minimum/maximum values of the tip transverse/longitudinal motion. The colour scheme in (a) is used in Figs. 5 and 6.

FIGURE 5 :FIGURE 6 :

 56 FIGURE 5: Phase portraits corresponding to Fig. 4. (a) … D 2:0, (b) … D 8:0, (c) … D 9:5, (d) … D 20:0.

  Generic infinitesimal element of the beam. The relative fluidbody velocity reads V rel D . P

		P w.x; t/k	Â.x;t / ˛.x;t / V rel .x; t/
	U i
		w.x; t /k P u.x; t/i
	.x;t/	u.x; t/i
	FIGURE 2: u.x; t / C U /i C P w.x; t /k.

As a consequence, .x; i/ .z; k/ w and u can be expressed in terms of the slope , also the angle of rotation of the cross-section, and the velocities become P w.x; t / D Z x 0 P .s; t / cos .s; t /ds; P u.x; t / D Z x 0 P .s; t / sin .s; t /ds:

  and Á s D Á= generates the following nonlinear dimensionless equation of motion (with the asterisk notation dropped for brevity):

								with ˛s D boundary conditions read .0; t / D 0 and x .L; t / D 0. atan. P w=. p u//. The clamped-free …= C P Equation (
	ˇR Á s P xx	Z x	 xx C …	cos	Z x 1 c f .s; t / cos .s; t / ds Ã
	C sin Â	1 ÄZ 1	c f .s; t / sin .s; t / ds Z s
	c d	sin				cos .Á; t / dÁds
				ÄZ 1	x Z s	0	Ã	t
	C cos		x ÄZ 1	0 Z s	sin .Á; t / dÁds	t
	sin							cos .Á; t / dÁds
				x ÄZ 1	0 Z s	t t
	C cos				x		0	sin .Á; t / dÁds	t t	D 0	(7)
	where
	c f D k a	h	…	R w cos	R u sin C 2 P P u cos C P w sin
	C x P w 2 Á C 2 P cos 1 2 P u 2 3 P u P w cos sin C r C … x P u 3 P u cos 2 Á 1 2	3 2	P w 2 P u 2 cos 2 w sin cos 3 P

x 3 cos 2 1 i C k p sin ˛s cos ˛s C k v sin ˛sj sin ˛sj;

ACKNOWLEDGEMENTS

The authors gratefully acknowledge the financial support from the Natural Sciences and Engineering Research Council of Canada.