Jicun Wang

Dr Thomas M Michelitsch

Arne Wunderlin

Ravi Mahadeva

AGING AS A CONSEQUENCE OF MISREPAIR --A NOVEL THEORY OF AGING

Keywords: Aging (ageing), Damage, Misrepair, Longevity

It is now increasingly realized that the underlying mechanism which governs aging (ageing) is a complex interplay of genetic regulation and damage-accumulation. "Aging as a result of accumulation of "faults" on cellular and molecular levels", has been proposed in the damage (fault)-accumulation theory. However, this theory fails to explain some aging phenotypes such as fibrosis and premature aging, since terms such as "damage" and "fault" are not specified. Therefore we introduce some crucial modifications of this theory and arrive at a novel theory: aging of the body is the result of accumulation of Misrepair of tissue. It emphasizes: a) it is Misrepair, not the original damage, that accumulates and leads to aging; and b) aging can occur at different levels, however aging of the body takes place necessarily on the tissue level, but not requiring the aging of cells/molecules. The Misrepair-accumulation theory introduced in the present paper unifies the understanding of the roles of environmental damage, repair, gene regulation, and multicellular structure in the aging process. This theory gives explanations for the aging phenotypes, premature aging, the difference of longevity in different species, and it is consistent with the physical view on complex systems.

INTRODUCTION

Many theories have been proposed to answer the questions of why and how we age. These theories fall into two main groups: genetic regulation and damage-accumulation. The gene-controlling theories emphasize the genetic controlling on aging (ageing). Many genes have been found to be related to aging, and some have proven to have lifespan-extending effects on animal models by gene knock-out (or -in), however their associations with aging are debated and far from confirmed. The questions are, if there are such kinds of genes: a) how their functions are related to the variety and inhomogeneity of the phenotypes of aging? b) why a living being has such kinds of genes to make itself going to die? and c) since when (which age) they exert functions badly? Among these theories, the Mutation-accumulation Theory assumes that the weak selection by nature in late age allows a wide range of gene mutations to accumulate, with deleterious effects [START_REF] Medawar | An Unsolved Problem of Biology[END_REF][START_REF] Martin | Genetic analysis of aging: role of oxidative damage and environmental stresses[END_REF]: one interesting point which remains unanswered is how cells can survive and function well when genes mutate with deleterious effects. The Developmental Theory suggests that aging and development are coupled and regulated by the same mechanisms [START_REF] Zwaan | Linking development and aging[END_REF][START_REF] De Magalhaes | Genomes optimize reproduction: aging as a consequence of the developmental program[END_REF]. It is interesting and reasonable; however this theory does not explain how the development process influences the aging phenotypes. The Damage-accumulation theories imply that it is the extrinsic and intrinsic damage (or fault) -accumulation on cellular or molecular level that leads to aging: The causes can be free radicals, byproducts, molecular cross-linking and so on [START_REF] Kirkwood | Why do we age[END_REF][START_REF] Holliday | The multiple and irreversible causes of aging[END_REF]. However, most of these theories actually mainly pointed out the phenotypes or causes but not the underlying common mechanism of aging.

It is now increasingly realised that genetic regulation and damage-accumulation actually interplay to cause aging, and it is predicted that: "aging is a result of accumulation of "faults" at cellular and molecular level because of the limitation of maintenance and repair; the underlying driving force is damage. The genetic control of longevity comes through the regulation of the essential maintenance and repair processes that slow the build-up of faults" [START_REF] Kirkwood | Ageing: too fast by mistake[END_REF]. However this recent damage (fault)-accumulation theory does not clarify the concepts of "damage" and "fault", and is weak in explaining some aging phenomena, such as fibrosis and premature aging. Therefore we suggest here some crucial modifications of this theory and raise a novel theory of aging.

A NOVEL THEORY OF AGING

Limitations of the current damage (fault)-accumulation theory

In the above mentioned damage (fault)-accumulation theory [START_REF] Kirkwood | Ageing: too fast by mistake[END_REF], the terms of "damage" and "fault" are not clearly defined. The "damage" can mean the primary damage before reparation and/or the incomplete repair following damage. The damage can be extrinsic or intrinsic. Incomplete repair may also include partial repair and incorrect repair. However it is necessary to distinguish between them. For example, the burnt skin develops a scar. The burnt skin is the damage, and the scar is the incorrect repair. The "fault" seems to mean the original damage which is either partially repaired or not at all repaired, by mistake due to the limitation of repair system. However, if the damage is left unrepaired, it will develop a life-threatening condition, such as bleeding, infection and organ failure, by destroying the integrity of structure, no matter whether the damage is on the molecular, cellular, or tissue level. Therefore, the original damage cannot remain unrepaired, and hence the concept of accumulation of "faults" is misleading. Reparation of any damage, and reconstruction of the structure, is essential to maintain the integrity of structure and the basic function for survival: even if sometimes it is imperfect and lead to a reduction of function. In most cases, reparation is complete, so that the structure can be restored. However, incorrect repair can and even must occur in cases of serious or frequent injuries, where quick reparation is more important for the immediate survival than perfect reparation. Typical examples include the Misrepair of DNA after double-strand breaks of DNA [START_REF] Kevin | Genomic DNA is captured and amplified during double-strand break (DSB) repair in human cells[END_REF] and the scar formation. The term "Misrepair" appears applicable in a wider sense, not only to DNA, but also to other types of incorrect repair due to a similar mechanism.

The concept of Misrepair

We define the concept of Misrepair as the incorrect reconstruction of a structure, after repairing the original damage, which can lead to a change of structure (of molecule/cell/tissue) and a reduction of function. The change of structure can be the change of the kinds or amounts of functional components, and their locations or spatial relationships. The comparison between damage and Misrepair is shown in Table 1. We define the term "damage" as the structure change before any reparation has taken place, whereas "Misrepair" is the change of structure as a result of incomplete or defective reparation. Distinguishing these two concepts is necessary to understand the aging mechanism, in which damage is the cause and Misrepair is the effect. Misrepair is unavoidable when serious or frequent damage occurs. For example, the non-stopping mechanical movements of heart wall, arterial wall and skin often cause the break of muscle fibres, elastic fibres and other extracellular matrixes (ECMs). Misrepairs on different structures appear in different manners, such as change of DNA sequence, rearrangement of cytoskeleton [START_REF] Garcia | Age-Dependent Defects in TCR-Triggered Cytoskeletal Rearrangement in CD4 + T Cells[END_REF], Mallory body formation, rearrangement of ECMs [START_REF] Cook | Defective Extracellular Matrix Reorganization by Chronic Wound Fibroblasts is Associated with Alterations in TIMP-1, TIMP-2, and MMP-2 Activity[END_REF], cellular reorganization [START_REF] Eliasieh | Cellular reorganization in the human retina during normal aging[END_REF], and tissue fibrosis. Most of the degeneration changes, such as hyaline degeneration (Mallory body), mucoid degeneration, pigment particle, and calcium homeostasis are in fact the products of Misrepair. [START_REF] Bansal | Defective membrane repair in dysferlin-deficient muscular dystrophy[END_REF][START_REF] Cenacchi | Ultrastructural changes in dysferlinopathy support defective membrane repair mechanism[END_REF][START_REF] Eliasieh | Cellular reorganization in the human retina during normal aging[END_REF]. A tissue Misrepair hypothesis was also raised to explain the radiation carcinogenesis [START_REF] Kondo | Tissue misrepair hypothesis for radiation carcinogenesis[END_REF]. However the importance of Misrepair mechanism was not sufficiently emphasized and its association with aging has not yet been studied. Here we discuss mainly the importance of Misrepair in aging.

Misrepair occurs no matter how capable the gene-regulated reparation and -maintenance system is, since reparationmany steps of biochemical reactionsneeds time. Misrepair is a result of the active reparation processes beneficial for survival, but not by mistake due to the limitation of repair system. Therefore the maintenance and repair system does NOT SLOW but even PROMOTES the occurrence of Misrepair if necessary, although a more capable repairing system means lower rate of occurrence of Misrepair. Hence the appearance of Misrepair gives more chance for survival of individual, by which the individual can live at least up to the reproduction age, which is important for the survival of species. Therefore the Misrepair mechanism was selected by nature due to its evolutionary advantage.

Aging is the result of accumulation of Misrepair

A misrepaired structure cannot be "repaired" any more (restored to the pre-damage state) and persists as a "scar". So it is the Misrepair and not the original damage that accumulates. Therefore the term of Misrepairaccumulation is more appropriate than that of damage-accumulation (Figure 1). Damage drives the aging process by triggering Misrepair. The accumulation of Misrepair leads to the gradual de-organisation of structure and the gradual reduction of function. Liver cirrhosis, for instance, is a result of accumulation of misreconstruction of tissue due to the death of hepatocytes in chronic hepatises. Aging is a result of reparation, but not a result without reparation. Aging is necessary for body immediate survival and species survival. This is why evolution favours aging mechanism. Thus Misrepair might represent the mechanism by which organisms are not programmed to die but to survive (as long as possible) [START_REF] Kirkwood | Understanding the odd science of aging[END_REF], and aging is just the price to be paid. A structure (molecule, cell or tissue) containing Misrepair has less repair capacity as well as other reduced functions, therefore a misrepaired structure is more fragile to damage, and the appearance of further Misrepair is more likely. So the accumulation of Misrepair is a self-accelerating process. This is the reason why aging phenotypes are not homogenous, which we can understand from the aging dots on the face.

Aging of tissue

Aging can occur at different levels: molecular, cellular and tissue. However, an interesting point is whether the aging of tissue, which is directly related to organ function and body lifespan, should be always due to the aging of cells or molecules. The intact function of tissue not only depends on the intact functions of cells/ECMs, but also depends on the intact collaborations (spatial relationships) among them. If a defect from damage or Misrepair) occurs to the tissue structure, which affects the spatial relationships between cells/ECMs, the defect could cause directly a decline of tissue function. If, however, the defect occurs to the cells/ECMs, there can be different outcomes depending upon the location. In tissue that can regenerate, the defected cells/ECMs can be removed and replaced with new cells/ECMs, and this might not necessarily affect tissue function. In tissue that cannot regenerate, such as nerve tissue, the defective cells not only affect the cell-contributed tissue function, but also interrupt the spatial relationships amongst other functional cells/ECMs. Therefore, the residual non-functional cells/ECMs can be also regarded as a defect of tissue structure. For instance, in Alzheimer"s disease, it has been realised that the corticocortical disconnection due to the loss of pyramidal neurons, which leads to the communicating failure between different parts of brain is the pathophysiological mechanism of Alzheimer"s disease [START_REF] Delatour | Alzheimer pathology disorganizes cortico-cortical circuitry: direct evidence from a transgenic animal model[END_REF].

The change of tissue structure -the change of the spatial relationship between cells/ECMs, is sufficient for the reduction of tissue function, without requiring the change of the cells or molecules themselves. Therefore aging of body takes place necessarily on tissue level. It is well-known that the common feature of body aging is the deleterious change of tissue structure. Aging does not start only when we are very "old", but it is a lifelong gradual process. One reason why some aged cells can be found in old regenerable tissue is that, the aged cells cannot be replaced in time, or completely, due to the reduced repair function of aged tissue. In this case, the aged cells are the effect of aging of tissue, but not the cause. In this aspect our view is the same as that in the Multicellular Being Chaos Theory, which suggests that it is the failure of information transmission in multi-cellular beings between each part that leads to aging [START_REF] Mulá | Why do we age?[END_REF], however, this theory does not explain how the breakdown of information transmission occurs.

The novel theory and its consistence with the physical view on complex system

We propose that aging is the result of accumulation of Misrepair of tissue. Distinguishing between damage and Misrepair clarifies what is the cause (damage) and what is the effect (Misrepair) in the aging process.

Misrepair as an imperfect repair links the functions of gene-controlling and damage-exposure in aging. The concept of tissue Misrepair couples aging and development, which obey the same mechanism: tissue-(re)construction. Therefore this theory unifies the multi-factors in aging mechanism.

Likewise, we also predict that the aging of a cell is the result of accumulation of Misrepair of intracellular structure. Many of degeneration forms of cells, such as fatty change, Mallory body, and pigment particle, are kinds of Misrepair, in which the dead components are not degraded, but isolated to prevent their toxicity to cells. The destiny of an aged cell can be, or cleared and replaced by new cell, or isolated by fabric structure when un-degradable, or to transform later into a tumour cell if a critical mutation on DNA occurs.

Our interpretation is consistent with the physical view on complex system. In terms of the Second Law of Thermodynamics, entropy always increases in a closed system. Living system overcomes the entropy increase by metabolism and exchanging energy/substance with environment. Entropy actually characterizes the degree of disorganization. Living beings continuously need to release entropy to protect them from running into thermodynamic equilibrium (death) (Niedermueller, 1990[START_REF] Haken | Grundprinzipien der Selbstorganisation[END_REF][START_REF] Wunderlin | Evolution of Dynamical Structures in Complex Systems[END_REF]. However the entropy increases unavoidably due to the un-removable Misrepair. Therefore it is the Misrepair that accounts for the permanent increase of entropy. In addition, in the spirit of the physical theory of complex systems [START_REF] Haken | Grundprinzipien der Selbstorganisation[END_REF][START_REF] Wunderlin | Evolution of Dynamical Structures in Complex Systems[END_REF], living beings must be conceived as complex self-organised systems [START_REF] Kiss | Ageing as a price of cooperation and complexity: Self-organization of complex systems causes the ageing of constituent networks[END_REF]. The characteristics of complex systems are their "emergent" properties and functions, which are qualitatively different from those of their subsystems. Hence some properties of living beings such as longevity and aging cannot be only explained by the properties of their subsystems such as cells and molecules [START_REF] Richter | Komplexe Systeme[END_REF].

INTERPRETATION OF AGING PHENOTYPES AND LONGEVITY IN THE LIGHT OF THE NOVEL THEORY

Aging-related disorders

Fibrosis is a typical phenotype of aging, existing in many old tissues, such as lung, heart muscle and arterial wall, which is unexplainable by the damage-accumulation theory. However, in our view, fibrosis is the result of accumulation of over-produced ECMs or intracellular fibre-like components for repairing, essential for survival by limiting the damage, linking the break, isolating the un-degradable cells or by-products, and reconstructing the structure, which is more robust. For example, in Alzheimer"s disease, the neurofibrillary tangles in brain are kinds of fibrosis structures to isolate the damaged components within nerve cells and the amyloid plaques are kinds of fibrosis-like structures with amyloid to isolate the damaged nerve fibres.

Atherosclerosis is an aging-related chronic inflammation disease, and the formation of plaques is the result of the over-proliferation and accumulation of macrophages for clearing the lipid in the arterial wall and repairing the endothelium when the endothelium is frequently damaged due to the movements of arterial wall or blood flows. This could be one underlying mechanism of the effect of the caloric restriction on retarding aging [START_REF] Masoro | Overview of caloric restriction and ageing[END_REF]. Some studies showed that long-term low-dose of Aspirin is helpful to reduce heart attack [START_REF] Hung | Aspirin for cardiovascular disease prevention. Medical Issues Committee of the National Heart Foundation of Australia[END_REF][START_REF] Elwood | My health: whose responsibility? Low-dose aspirin and older people[END_REF]. The effect of Aspirin on inhibiting the Atherosclerosis might be one mechanism.

Aging-related tumours are the results of mutations due to the accumulation of Misrepair of DNA, which leads to the out-of-control of cell proliferation and finally to tissue disorganization. Actually in many cases, the mutations of DNA cause cell death, and not affect the tissue function. However, the unrepaired or misrepaired DNA damage has been implicated as a causal factor in cancer and aging [START_REF] Suh | Maintaining genetic integrity in aging: a zero sum game[END_REF]. It is now accepted that the misrepaired DSBs are the main lesions at the origin of both chromosomal abnormalities and gene mutations (Natarajan, 1993 and Bishay;[START_REF] Bishay | DNA damage-related RNA expression to assess individual sensitivity to ionizing radiation[END_REF]. Misrepair of DNA is one survival mechanism of cell under radiation [START_REF] Rothkamm | Misrepair of radiation-induced DNA double-strand breaks and its relevance for tumorigenesis and cancer treatment (review)[END_REF]. The lack of cancer in simpler organisms [START_REF] Campisi | Aging, chromatin, and food restriction--connecting the dots[END_REF] might be because their lifespans are too short to accumulate sufficient Misrepairs of DNA.

Wrinkle formation

It is well-known that the elastic structures are gradually replaced by collagen structures with aging in many organs, such as airway wall, stomach wall, blood vessel wall and skin. On skin this change leads to the formation of wrinkle. Figure 2 shows a schematic representation of the changes of elastic tissue in the process of wrinkle formation. Firstly when the elastic fibres in skin are repeatedly broken, some of them are replaced by collagen fibres by Misrepair; further when the collagen fibres are repeatedly broken during shrinking, some of them are replaced by shorter fibres (Misrepair); with time, more and more elastic fibres are replaced by collagen fibres with different lengths (Misrepair-accumulation), which gradually leads to the reduction of extension potential and the permanent folding of some longer fibres-wrinkle formation (aging).

Premature (accelerated) aging

Development and repair are both tissue/body structure formation processes. Any element that interrupts the construction process, e.g. cell division and ECMs production/formation, will affect the development/repair and lead to mis-construction (imperfect construction) and mis-reconstruction (Misrepair) of tissue/body. At this point, the aging and the development procedures are regulated by a similar genetic mechanism [START_REF] Zwaan | Linking development and aging[END_REF][START_REF] De Magalhaes | Genomes optimize reproduction: aging as a consequence of the developmental program[END_REF]. This mis-construction and mis-reconstruction could be the mechanism of premature (accelerated) aging due to gene mutation. In the pre-matured body, the imperfect tissue structure is with lower function and lower potential, and it goes to failure faster. In this case, it is the accumulation of misconstruction (abnormal development) and accumulation of Misrepair together that lead to the failure. Most such individuals cannot survive until full development. The lower repair function, due to gene mutation and abnormal tissue structure, causes more frequent Misrepair. For example, the Hutchinson-Gilford Progeria is caused by a gene mutation in the nuclear envelope protein Lamin A [START_REF] Dechat | Alterations in mitosis and cell cycle progression caused by a mutant lamin A known to accelerate human aging[END_REF]. This change of lamin A results in the nuclear blebbing, which interrupts the DNA/RNA synthesis and affects the cell division. In Fibulin-5-/-aging, mice with mutation of Fibulin-5 gene have loose and wrinkled skin, vascular abnormalities, and emphysema, all of which are thought to be due to their disorganized and fragmented elastic fibres, which interrupting the development and repair [START_REF] Hirai | Fibulin-5/DANCE has an elastogenic organizer activity that is abrogated by proteolytic cleavage in vivo[END_REF]. In Premature klotho-/-aging, Kloto gene mutation in mice causes multiple accelerated aging changes due to the disorder of calcium homeostasis, since Klotho gene encodes a type I membrane protein with homology to beta-glycosidase, which is important in the regulation of calcium homeostasis [START_REF] Lanske | Premature aging in klotho mutant mice: cause or consequence?[END_REF].

Genes restrict the maximum lifespan by shaping the body

It is unavoidable that tissue structure goes finally to breakdown due to damage and Misrepair; therefore, the potential of tissue is the key point for longevity. Death often happens from the failure of a key tissue/organ, e.g. the blockage of blood vessel in brain and heart due to Atherosclerosis. The maximum potential of tissue lies in the complexity and maintenance of the tissue structure. More complex structures often have more compensational potential in function due to their repeated and network structures. This could be the reason why species with bigger body/brain size often, although not always, has longer longevity (Aziz, 2005). It will be interesting to compare the structural complexities between different species" to explain the difference of longevity. Genes predetermine this potential by controlling the development. Therefore gene configuration restricts the maximum lifespan by shaping the body. Longer living species" often have a longer development time for their more complex functionality and need longer time for accumulating sufficient Misrepair. This may be the mechanism explaining why different species" have different maximum life-spans. An exception is the big difference in the lifespan between the queen and the worker ants (also in the case of bees): in spite of the same gene configurations. However, the queen and the workers undergo different developments, which lead to different body structures. In C. elegans, an alternative (dauer) developmental pathway results in a significant longer lifespan [START_REF] Klass | Non-ageing developmental variant of Caenorhabditis elegans[END_REF]. These examples indicate that it could be the difference of body structure that finally and directly determines the difference of maximum longevity. Some plants such as most trees, and animals such as ant queen, continue to develop after reproduction, therefore they obtain additional structure complexity and additional function. This increase of complexity enables them to live longer.

It is possible that nature develops longer-living species by increasing the body complexity; however, if the lifespan is too long, individuals would die before the reproduction age, which would lead to the extinction of the species. For example, the Ginkgo biloba tree is reported to be able to live more than 3000 years, and it can survive in extreme weather and environment; however this species is an endangered plant. The survival of species is the result of limited longevity and it is the best evidence of the essentiality of limited longevity.

REDUCING DAMAGE-EXPOSURE IS THE FEASIBLE WAY TO DELAY AGING-RELATED DISEASES

Although gene configuration determines the maximum longevity of a species, the individual lifespan is more related to the environmental damage-exposure, and it can only be extended to a certain degree. Damage triggers cell reaction and inflammation, and the repeated damage-caused chronic inflammation is important source of Misrepair. Our theory suggests that for extending lifespan all efforts need to focus on the reduction of Misrepair. This can be achieved by avoiding damage exposure as far as possible, by improving living conditions and medical care. It is especially important to prevent repeated damage and chronic inflammation.

To have a regular daily life including healthy diet, sufficient rest and sleep, is important for good reparation. At the same time, we should positively understand aging/Misrepair; without aging/Misrepair, we would have already died from once damage in young age.

An "improved" repair system by gene-manipulation would theoretically reduce the rate of Misrepair. However the repair-related genes are numerous and they need to act collectively to get their function. Up-or downregulating of one or two such genes would be insufficient for the goal to improve repairing. In this aspect, the living beings themselves could know better through evolution "what is the best". Manipulating these genes without completely understanding their functions could cause fatal danger from side effects. Those genes involved in repairing are also important in development. The extending of lifespan by Knock-out or Knock-in of some genes in animal model such as yeast and C. elegans, is probably through redirecting the development [START_REF] Wei | Life span extension by calorie restriction depends on Rim15 and transcription factors downstream of Ras/PKA, Tor, and Sch9[END_REF], and this alteration can potentially cause fatal defects of development (hypoplasia).

Anti-inflammation could be a useful intervention to judge whether Misrepair is the source of aging. From our view, anti-inflammation would reduce the chance of Misrepair, although it would leave some damage to be repaired in delay, which is dangerous. Therefore theoretically the long-term application of low-dose of antiinflammation medicine would be helpful to reduce aging-related diseases. Some investigations have shown that Aspirin, an anti-inflammation medicine, has effect of "anti-aging", including cardio-protection and reduction of cancer [START_REF] Elwood | My health: whose responsibility? Low-dose aspirin and older people[END_REF]. In addition, Bulckaen et al have observed the effect of aspirin on anti-aging on animal model [START_REF] Bulckaen | Low-dose aspirin prevents age-related endothelial dysfunction in a mouse model of physiological aging[END_REF]). However the effect is not yet fully confirmed [START_REF] Kang | Low dose aspirin and cognitive function in the women's health study cognitive cohort[END_REF]. More experiments on animal models are required. One difficulty of such research is the disturbance of the side effects of anti-inflammation to the research result. It has to be emphasized that the side effects of antiinflammation could cause fatal problems such as bleeding, infection, and organ failure due to insufficient defense and repair. Therefore the medicine, even if it is effective, should be only recommended to some older people, at least more than 50 years old, who have high risk of diseases such as vascular diseases or cancer.

CONCLUSION

In this paper we have proposed for the first time the link of Misrepair with aging. Aging is the result of accumulation of Misrepair, and aging of body is necessarily on tissue level. Misrepair and aging is active and necessary for immediate survival and species survival. This theory helps to answer the questions, such as why and how we age (for immediate survival and species survival by Misrepair), why aging mechanism has advantage in evolution, and why we have limited longevity. The following predictions can be made from this theory:

 Aging of cell is the result of accumulation of Misrepair of intracellular structure.  The longer living species" have more complex body/tissue structures, which need more time to develop.  Inhibiting chronic inflammation would help reduce aging-related diseases, but it would be accompanied with high risk of other diseases.

Figure 1 .

 1 Figure 1. Comparison of Misrepair-accumulation and damage-accumulationDamage-accumulation: due to the limitation of maintenance and repair, some damage remains partially repaired or unrepaired by mistake. This remained damage (fault) will accumulate and lead to aging. Misrepair-accumulation: serious or too frequent damage triggers Misrepair in emergency to maintain the structural integrity and protect from immediate death. However the persistent Misrepair will accumulate and lead to aging.

Figure 2 .

 2 Figure 2. The formation of wrinkleFirstly when the elastic fibres are repeatedly broken, some of them are replaced by collagen fibres by Misrepair. Further some of the collagen fibres are replaced by shorter ones (Misrepair) when the brokens take place during shrinking. With time more and more elastic fibres are replaced by collagen fibres with different lengths (Misrepair-accumulation), which lead to the reduction of extension potential and the permanent folding of some longer fibres -wrinkle.

Table 1 . Comparison of damage and Misrepair Damage Misrepair

 1

	Physical and chemical injury on molecule/cell/tissue	Imperfect reconstruction of structure of molecule/cell/tissue
	Defect before reparation	Defect after reparation
	Removable	Permanent
	e.g. Injury or death of molecule/cell/tissue	e.g. Change of the kinds or amounts of functional components and
		their locations or spatial relationships


Preconditioners based on Strong Components

List of authors: K. Kaya1 I. S. Duff 2 We propose a new method for constructing a preconditioning matrix M to accelerate the solution of the system Ax = b, when using Krylov-based iterative methods. The coefficient matrix A is large and sparse and we assume it is irreducible, that is it cannot be permuted to block triangular form. If it is reducible, then we will apply our algorithms to the irreducible blocks on the diagonal.

The proposed method is based on a hierarchical decomposition of the associated digraph into strongly connected subgraphs. It permutes the rows and columns of the original matrix A and obtains a block triangular preconditioning matrix containing a subset of the nonzeros of A where the maximum size of a diagonal block is smaller than a desired value.

Let G = (V, E) be a strongly connected digraph with n vertices and m weighted edges. A hierarchical decomposition of G into its strong subgraphs can be defined in the following way. Let σ 0 be a permutation of the edges. For 1 ≤ i ≤ m, let σ 0 (i) be the ith edge in σ 0 . Let G 0 = (V, ∅) be the graph obtained by removing all the edges from G. Consider that edges are added one by one to G 0 in the order determined by σ 0 . Let G i = (V, {σ(j) : 1 ≤ j ≤ i}) be the digraph obtained after the addition of the first i edges. Initially in G 0 , there are n strong subgraphs, one for each vertex, and, during the edge addition process, these strong components gradually coalesce until there is only one. The hierarchical decomposition of G into its strong components with respect to the edge permutation σ 0 shows which strong components are formed in this process hierarchically. Note that a strong component in the edge addition process is indeed a strong component of some digraph G i although it is only a strong subgraph of G.

More explanation by diagrams can be seen in our report [START_REF] Duff | Preconditioners based on strong components[END_REF].

Given a digraph G = (V, E) and a permutation σ 0 , the hierarchical decomposition of G can be obtained by first constructing G 0 and executing Tarjan's strong component algorithm (SCC) for each internal digraph G i obtained during the edge addition process. This would be an O(mn + m2 ) algorithm since 1 ≤ i ≤ m and the cost of SCC is O(n + m). It would thus be prohibitive for large graphs. To find the decomposition in a more efficient way, Tarjan first proposed a recursive algorithm of complexity O(mlog 2 n) that he later improved to have complexity O(m log n) [START_REF] Tarjan | An improved algorithm for hierarchical clustering using strong components[END_REF].

Note that there is a one to one correspondence between the blocks in a block triangular form of an n × n matrix with m nonzeros and the strong components of its associated digraph with n vertices and m weighted edges. The proposed method uses Tarjan's hierarchical decomposition algorithm to remove some edges with relatively small weights and obtain a reduced digraph whose associated matrix can be permuted into a block triangular form such that the dimension of each block is smaller than a desired value.

The first step in our proposed algorithm for creating M is a modified version of Tarjan's algorithm, HD. Our first modification allows us to have non-distinct edge weights, that is matrices with the same value in different positions. Another modification to the HD algorithm is that we do not want to accept any blocks larger than a predefined value, mbs, and this enables us to stop the recursion earlier than in the original algorithm. We use this modified algorithm to obtain a block triangular matrix M where the strong components correspond to the blocks on the diagonal of M. To the best of our knowledge, this is the first work that uses Tarjan's algorithm for preconditioning purposes. After we obtain the block triangular form from this modified HD we then see whether any blocks can be merged and finally use a greedy algorithm to order the blocks on the diagonal so that most of the entries are in the upper triangular part.

Figure 1 shows a sample matrix and the preconditioner obtained by our algorithm, SCPRE.

The main parameters of our algorithm are the ordering σ 0 and the largest block size mbs. Before we use algorithm HD we first scale and permute the matrix using MC64 [START_REF] Duff | On algorithms for permuting large entries to the diagonal of a sparse matrix[END_REF], which we do also for our experiments on other preconditioners. We compare our algorithm SCPRE with another block preconditioner XPABLO [START_REF] Fritzsche | Extensions of certain graph-based algorithms for preconditioning[END_REF] and a version of the industry-standard ILUT from MATLAB, on sets of matrices from circuit (left half) and device (right half) simulations from the University of Florida sparse matrix collection. We show below a very abbreviated table from the results in [START_REF] Duff | Preconditioners based on strong components[END_REF]. In the table we give the number of iterations (with the least in bold font) and in the second line for each matrix, the relative memory requirement (ratio of storage for preconditioners to that for the original matrix). Although this is a much reduced set from our technical report, the results are representative in the sense that the conclusions from this set match those from the larger set. For this set of experiments, we use mbs = 1000. For circuit simulation problems, ILUT and SCPRE converge for all matrices in this set. XPABLO fails to converge for bcircuit and ckt11752 dc 1 and so SCPRE is clearly the best block based preconditioner on this set of matrices. Although ILUT requires significantly fewer iterations on G2 circuit and ckt11752 dc 1, in both cases it requires more memory. However, for G2 circuit, if we increase mbs to 5000, the number of iterations drops to 95 and our relative memory requirement increases to 6.10 and, for ckt11752 dc 1, by increasing mbs to only 3000 we require only 11 iterations with a relative memory cost of only 1.45. Thus we feel we can we recommend using SCPRE for circuit simulation matrices especially when the amount of memory to store the preconditioner is the main concern.

For the device simulation matrices in the right-hand side of the table, the block based preconditioners are far more robust on this set with convergence for all the test matrices. We therefore feel that we can recommend SCPRE as the preconditioner for the device simulation matrices.

To balance these excellent results, we show in our paper [START_REF] Duff | Preconditioners based on strong components[END_REF] that ILUT outperforms both block approaches on matrices from CFD applications. Further research is needed to understand the effect of structure in determining the best approach.

Figure 1 :

 1 Figure 1: Matrix mult dcop 01 (on the left) and the preconditioner obtained by the algorithm SCPRE the right) with mbs = 250.

CERFACS, France

Atlas Center, RAL, England


B Hou

S Pattofatto

Y L Li

H Zhao
email: zhao@lmt.ens-cachan.fr

Impact behavior of honeycombs under combined shear-compression, PartⅡanalysis

Keywords: honeycombs, combined shear-compression, FEM, dynamic enhancement

In this paper, a numerical virtual model of honeycomb specimen as a small structure is used to simulate its combined shear-compression behavior under impact loading. With Abaqus explicit code, response of such a structure made of shell elements is calculated under the prescribed velocities as those measured in the combined shear-compression test presented in the part I of this study.

The simulated results agree well with the experimental ones in terms of overall pressure/crush curves and deformation mode. It allows for the determination of the separated normal behavior and shear behavior of honeycomb specimen under dynamic combined shear-compression. It is found that the normal strength of honeycomb decreases with the increasing shearing load. Quasi-static calculations were also performed and a significant dynamic enhancement found in experiments was validated again in the numerical work. In the end, a crushing envelope of normal strength vs. shear strength plane was obtained on the basis of these simulations.

Introduction

Dynamic multiaxial behavior of honeycombs as a basic energy-absorption design parameter is eagerly desired in order to perform numerical simulations for various industrial applications. Many previous works on this domain have been reported in the open literature and a large number of these works concern mainly the in-plane behavior, and mostly under quasi-static loading [START_REF] Gibson | Cellular material: structure and properties[END_REF][START_REF] Gibson | Failure surfaces for cellular material under multiaxial loads-I. modeling[END_REF][START_REF] Klintworth | Elasto-plastic yield limits and deformation laws for transversly crushed honeycombs[END_REF][START_REF] Yang | Elastic buckling of regular hexagonal honeycombs with plateau borders under biaxial compression[END_REF][START_REF] Papka | In-plane biaxial crushing of honeycombs-Part II: Analysis[END_REF]. However, the most interesting behavior of honeycombs for an energy absorption application is the out-of-plane crushing behavior, especially the one under combined out-of-plane shear-compression which is the most realistic loading mode for such use.

Under quasi-static loading, some testing methods for the combined shear-compression loading have been reported [START_REF] Doyoyo | Microstructural response of aluminum honeycomb to combined out-of-plane loading[END_REF][START_REF] Mohr | Experimental investigation on the plasticity of hexagonal alumininum honeycomb under multiaxial loading[END_REF][START_REF] Hong | Quasi-static crush behavior of aluminum honeycomb specimens under compression dominant combined loads[END_REF] and were used to determine the yield envelope of aluminum honeycomb under this prescribed biaxial loading states [START_REF] Hong | Quasi-static crush behavior of aluminum honeycomb specimens under compression dominant combined loads[END_REF][START_REF] Mohr | Deformation-induced folding systems in thin-walled monolithic hexagonal metallic honeycomb[END_REF][START_REF] Mohr | Large plastic deformation of metallic honeycomb: Orthotropic rate-independent constitutive model[END_REF]. For example, Hong et al. [START_REF] Hong | Quasi-static crush behavior of aluminum honeycomb specimens under compression dominant combined loads[END_REF] derived a quadratic yield criterion suitable for orthotropic material by modifying Hill's quadratic yield criterion. Mohr et al. [START_REF] Mohr | Deformation-induced folding systems in thin-walled monolithic hexagonal metallic honeycomb[END_REF][START_REF] Mohr | Large plastic deformation of metallic honeycomb: Orthotropic rate-independent constitutive model[END_REF] suggested a linear fit for the crushing envelope based on their quasi-static calculated results.

Under dynamic loading, many reported works revealed that the strength of honeycomb under uniaxial dynamic compression is higher than under quasi-static loading [START_REF] Wu | axial crush of metallic honeycombs[END_REF][START_REF] Baker | static and dynamic properties of high-density metal honeycombs[END_REF][START_REF] Zhao | On the strength enhancement of rate insensitive square tubes under impact loading[END_REF][START_REF] Harrigan | Inertia effects in impact energy absorbing materials and structures[END_REF], even if the shock wave effect is not involved [START_REF] Elnasri | Shock enhancement of cellular structures under impact loading: Part I Experiments[END_REF]. However, the behavior of honeycombs under dynamic multiaxial loading is rarely reported up to now. The main reason for such situations lies in the difficulties to achieve dynamic multiaxial experiments with accurate data measurements. Some previous works proposed the dynamic multiaxial testing methods using drop-weight or high speed machine [START_REF] Hong | Dynamic crush behaviors of aluminum honeycombs specimens under compression dominant inclined loads[END_REF][START_REF] Chung | Compressive response of circular cell polycarbonate honeycombs under inplane biaxial static and dynamic loading. Part I:experiments[END_REF], but the accuracy is not optimal at higher loading rate. In order to improve the measurement accuracy, we proposed in the Part I of this paper a new testing method using large diameter Hopkinson bars with beveled ends to perform combined shear-compression test under impact loading. It permits to obtain interesting overall force/crush curves but cannot give a separate normal and pure shear behavior. Thus, with the test presented in Part I, there is no mean to identify directly a given yield criterion.

In this part II of the work, a numerical approach is presented to study a yield criterion.

The dynamic and quasi-static combined shear-compression experiments are numerically reproduced with a detailed FEM model for the honeycomb specimen. The accuracy of these simulations is validated by comparing the numerical results with the testing ones. Such virtual tests provide a separated normal and shear behavior of honeycomb specimen, which allow for the determination of the yield envelope depicted in the macroscopic shear strength vs. compressive strength plane.

Numerical tests of honeycomb under combined shear-compression

Since the study is focused on the behavior of honeycomb under combined out-of-plane shear-compression, the modelling of the whole testing environment is not necessary. Thus, only the detailed honeycomb structure was modeled here and the loading environment was modeled by two rigid planes moving at the velocities measured during real tests. The commercial FE code of ABAQUS/explicit was employed for this numerical work.

Numerical specimen

The honeycomb structure studied here has the identical geometry to the hexagonal honeycomb used in the experiments presented in part I. It is composed of single-thickness walls (or thin walls) and double-thickness walls (or thick walls), having the following geometric parameters: single wall thickness h=76 µm, expansion angle α=30 o , and minimum cell diameter S=6.35 mm (as shown in Figure 1).

A complete model possesses the same size as the specimen employed in the experiments, which includes 39 cells on the honeycomb cross section. The dimensions of the specimen are 25×40×40 mm in the directions of T, L and W respectively (Figure 1).

The thick walls in a real honeycomb are typically made of two single-thickness thin walls which are bonded together. In this model, we ignore the rare delamination of the bonded interfaces and consider the strength of the adhesive bond as infinite. Thus, the simulations are carried out for a monolithic honeycomb, where the thick walls are also represented by a single shell element layer but with a doubled thickness value.

The model is meshed with 4-node doubly curved thick shell elements with reduced integration, finite membrane strains, active hour-glass control and 5 integration points trough the cell wall thickness. In order to determine the appropriate element size, a convergence study was performed among element sizes of 1 mm, 0.5 mm, 0.25 mm and 0.125 mm. It seems that the results converge when the element size is equal to or below 0.25 mm. With selected element size of 0.25 mm, the complete-model has totally 232600 elements.

The numerical specimen is placed between two rigid planes moving with prescribed velocities that are defined with the real input and output velocities (denoted as V input and V output in Figure 2) measured during the combined shear-compression experiments (Part I of this study). In this model, general contact with frictionless tangential behavior is defined for the whole model excluding the contact pairs of rigid planes and tested honeycomb specimen, which are redefined by surface-to-surface rough contact to make sure that no slippage occurs. As the real honeycomb is always far from perfect, it includes all kinds of imperfections which affect the initial peak value markedly, but have little influence to the crush behavior at large deformation period. These imperfections can be categorized into global ones and local ones, like irregular cell geometry, uneven or pre-buckled cell walls, wall thickness variation etc. Here in this work, we generated the imperfections by preloading the perfect specimen uniaxially by 0.1 mm before applying the prescribed experimental velocities. The value of 0.1 mm is chosen to make sure that the simulated initial peak is same as the one from experimental curve at uniaxial compression.

Quasi-static simulations were almost impossible to achieve with ABAQUS/Standard which uses Newton's method (or quasi Newton's method) as a numerical technique due to the complex nonlinear effects, e.g. the geometrical and material nonlinearity, the complex contact conditions as well as the local instability during crush. An alternative is to use also ABAQUS/Explicit for quasi-static problems. However, the explicit integration scheme of dynamic simulation codes usually leads to very small time step which in our simulation is around ten nanoseconds for the chosen element size. Thus, with the loading velocity of 0.1 mm/s, the computational duration for the quasi-static simulation (e.g. 180 s) will be too large. To overcome this difficulty, automatic mass scaling technique was employed to increase the time increment to 100 µs. The quasi-static loading conditions are guaranteed by ensuring the ratio of the kinetic energy to the strain energy as a small value with the selected time increment.

A bilinear elasto-plastic material model was employed to describe the cell wall material of this aluminum honeycomb. The loading case of dynamic uniaxial compression is taken as calculating instance. The model parameters of base material such as yield stress and hardening modulus were determined by fitting the calculation result of uniaxial compression to the result from experiment (Table 1).

Simplified models

In order to reduce the calculation cost for this complete honeycomb model, numerical models with various simplifications can be also used. For example, some researchers [START_REF] Papka | In-plane biaxial crushing of honeycombs-Part II: Analysis[END_REF][START_REF] Hönig | In-plane dynamic crushing of honeycombs. Part I: crush band initiation and wave trapping[END_REF] used one-dimensional beam elements with different microsections for the simulation of the in-plane behavior, while some others employed one layer of shell element according to the repeated behavior in cell axis direction [START_REF] Zou | Dynamic crushing of honeycombs and features of shock fronts[END_REF]. As to the out-of-plane behavior, the honeycomb specimen was usually simplified into a unit cell or a row of cells because of its periodicity [START_REF] Mohr | Deformation-induced folding systems in thin-walled monolithic hexagonal metallic honeycomb[END_REF].

These simplifications may introduce some imprecisions to the numerical model. In order to check the potential errors, two simplified models were established. By comparing the results of these three models under uniaxial out-of-plane compression, the accuracy of the simplified numerical models were insured.

The so called row-model is made up of a row of cells based on the periodicity of honeycomb specimen in L direction (as shown in Figure 4) and will be used to investigate the combined shear-compression behavior of honeycomb in TW plane. The most simplified model consists of three conjoint half walls in "Y" configuration (denoted as cell-model as in The simplified models work with symmetric boundary conditions. These displacement constraints are applied to the row-model on the two boundaries in L-direction (as shown in Figure 4). For the cell-model, symmetric boundary conditions are performed on the three non-intersecting edges of each cell wall in local y-direction (as shown in Figure 4). The same method is employed to introduce imperfections into these simplified models.

Biaxial behavior of honeycomb under combined shear-compression

In this section, the results of four loading cases under combined out-of-plane shear-compression (in TW plane) with loading angles of 30 o , 40 o , 50 o and 60 o simulated with row-model for both dynamic and quasi-static cases are presented. We first carefully examined the validity of our numerical specimens by comparing the overall pressure/crush curves and the deforming mode with experimental results at various loading angles. Finally, the normal and shear behavior of numerical honeycomb under combined shear-compression could be separated.

Validation of numerical specimen

The overall pressure/crush curves are obtained from the calculated results in order to make a comparison between the experiments and the simulations. It is worth emphasizing that the variable "crush" is defined in Part I as the relative displacement component of the two moving bevels in X 3 direcion and the "pressure" as the X 3 force component divided by specimen cross-sectional area S s . As a consequence, in the case of numerical combined shear-compression test, the pressure P(t) is calculated by dividing the contact force component in the rigid planes moving direction (X 3 direction in Figure 6) with specimen area S s . Its relationship with the directly obtained normal and shear contact forces is as follows:

s s n S t F t F t P / ) sin ) ( cos ) ( ( ) ( θ θ + = (1),
where θ is the loading angle as defined in Part I. F n (t) and F s (t) are respectively the normal and shear contact forces at the interfaces of honeycomb specimen and rigid loading planes.

The overall crush Δ(t) is derived from the relative resultant displacement of the two reference points on rigid planes (Figure 6). It has a simple relationship with the normal and shear crush (denoted as d n (t) and d s (t)), which is:

θ θ sin / ) ( cos / ) ( ) ( t d t d t s n = = Δ (2)
Figure 7 presents the experimental and calculated pressure/crush curves under dynamic uniaxial compression and combined shear-compression, for a representative loading angle θ=50 o . Each curve has two distinct stages. During stage I (from the beginning of zero crush to the position of the initial peak), the slopes of elastic segment are in good agreement for the calculations and the experiments. In addition, with the employed magnitude of imperfection which is determined in complete-model under uniaxial compression, the initial peak of the calculated curve also agrees well with the experimental one for the displayed combined loading instance equally. During stage II (defined as the following crush period after stage I to 13 mm crush), the experimental curves have more fluctuations than the calculated ones, but still can be predicted in terms of average crush strength.

A comparison between the initial peak value as well as the average crush strength for every loading angle is described in Figure 8. The average crush strength is calculated by dividing the curve area of this plateau deforming region (absorbed energy) by the corresponding crush length (same formulas as used for experimental curves), which gives:

∫ - = max * * max 1 δ δ δ δ δ pd p (3)
where * δ denotes the crush value at the point of the initial peak for each of the overall pressure/crush curve. max δ is the maximum crush of the concerned crushing duration.

Figure 6 Scheme of the decompositions of force and crush

X 2 X 3 |E F 2 F 3 F s F n RP 1 RP 2 | d n d s X 2 X 3 |E F 2 F 3 F s F n RP 1 RP 2 | d n d s
experiment is found at loading angle of 50 o . For the average crush strength, the deviation is a little more significant at larger loading angles. The quasi-static virtual testing results for uniaxial compression and combined shear-compression at loading angle θ=50 o are compared in Figure 9 together with the experimental curves. The numerical results show a good correlation with the experimental ones at the crushing stage II. The average crush strengths are calculated for all the loading angles and listed in Figure 10. The maximum error of 13.2 % is found at loading angle of 40 o . Nevertheless, during stage I of the curves, a clear difference is found for both the ascending segment slope and the initial peak value that can be attributed to a slight slippage between specimen and bevels at the beginning of the experiment.

Furthermore, the simulated deforming pattern of honeycomb under combined shear-compression during stage II is also compared with the experimental results obtained with high speed camera. Figure 11 It can be seen that the cell wall axis of all the specimens displayed incline during the crushing processes, and the inclined directions of the virtual and real specimens are in parallel to each other for the two loading angles. Besides, the phenomenon of two-side folding system is also found in the numerical results as discovered for most of the experimental shear-compression specimens. However, the numerical result can not cover at the same time the two deforming modes (rotation of cell axis or not) found during experiments (Part I). Figure 12 illustrates the rotation angle α at 40 % mean compressive strain for every loading angle. It appears that our numerical specimen have a clear preference for the deforming mode with significant cell axis rotation.

Multiaxial behavior

The validation of the simulation work in section 3.1 shows that the virtual testing results can represent well the experimental ones with the exception of quasi-static initial peak forces.

These virtual combined shear-compression tests provide more information than the real experiments and enable us to study the normal and shear behavior of honeycomb separately.

The separated normal and shear pressure/crush curves under dynamic loading are shown in Figure 13 Under quasi-static loading, the normal and shear pressure crush curves show a great similarity to the dynamic ones, i.e. the normal strength decreases with the loading angle whereas the shear strength increases. Moreover, a comparison between the quasi-static and the dynamic curves shows that the loading rate will also affect the normal and shear behavior of honeycomb under combined shear-compression. Figure 14 displays the dynamic and The average crush strength of normal and shear behavior were calculated for both dynamic and quasi-static loading at every loading angle as done in section 3.1 for the overall pressure/crush curves. All these average values are collected in Figure 15, which shows clearly the change of the normal and shear crush strength along with the loading angle as well as a strength enhancement under impact loading for every loading angle.

Macroscopic yield envelop estimation

Figure 16 shows the distribution of calculated honeycomb biaxial behavior on the normal average strength vs. shear average strength plane during the stage II. An elliptical shape is found for both the quasi-static and dynamic loading cases (Eq. 4).

1 2 0 2 0 =         +         τ τ σ σ , (4)
where 0 σ and 0 τ are respectively the normal crushing strength under uniaxial compression and the shear crushing strength under pure shear loading. By fitting the data with Levenberg-Marquardt algorithm (LMA), these two parameters are identified to be 3.98 MPa and 1.11 MPa under dynamic loading and 3.57 MPa and 1.02 MPa under quasi-static loading.

It is found in Figure 16 that the expansion of the crush envelope from quasi-static to dynamic loading is almost isotropic, even though the normal strength/shear strength ratio for one same loading angle is different under quasi-static and dynamic loading. It means that the dynamic biaxial strength for this honeycomb might be derived by using the enhancing ratio of uniaxial compression and the quasi-static crush envelope. A numerical specimen was built with simplified row model and identified material parameters and it provided simulated results in good agreement with the experiments in terms of deforming mode and the overall pressure/crush curves which were the final obtained information from the new designed combined shear-compression loading experiments.

Such numerical virtual tests enabled to separate the normal and shear behavior of honeycomb. It showed that the strength of honeycomb under compression is largely affected by the additional shear loading and exhibits a significant decrease while increasing shear loading. An obvious enhancement was also observed at dynamic loading for both the normal and shear behavior with respect to the quasi-static case at every loading angle.

In order to describe the dynamic and quasi-static biaxial behavior of honeycomb at macroscopic level, an elliptical criterion in the plane of normal strength vs. shear strength can be derived with a set of parameters obtained by fitting the data with Levenberg-Marquardt algorithm. The expansion of the crush envelope with the loading rate happened to be isotropic for the studied honeycomb in combined out-of-plane shear-compression in T and W direction.

Figure 1

 1 Figure 1 The geometry of a unit cell and the constructed honeycomb specimen

Figure 2

 2 Figure 2 Scheme of loading velocities

Figure 3 Figure 3

 33 Figure3shows the comparison between experimental and simulated pressure/crush curves, which validates the parameters of this bilinear material model.

Figure 4 )

 4 Figure 4) and can be used only in uniaxial compression to make a comparison with the other two models. Both of the two simplified models have a length of 25 mm in T direction, and the same element size of 0.25 mm as in the complete model. The numbers of elements for row-model and cell-model are 28500 and 2100 respectively.

Figure 4 Figure 5 Figure 5

 455 Figure 4 Scheme of complete and simplified models

Figure 8 Figure 7

 87 Figure 8 Comparison of the initial peak and average crush strength between dynamic calculations and experiments at various loading angles

Figure 9

 9 Figure 9 Comparison of the quasi-static pressure/crush curves from calculation and experiments

  (a) and (b) show the specimens at dynamic loading of θ=30 o and at crush of 12 mm, and Figure 11 (c) and (d) for the quasi-static loading of θ=50 o .

Figure 10 Figure 11

 1011 Figure 10 Comparison of average crush strength between quasi-static calculations and experiments at various loading angles

  (a) and (b) respectively. It is noted that the normal and shear pressures are calculated from the normal and shear contact forces (F n (t) and F s (t)) at the interfaces between rigid loading planes and honeycomb specimen. For the sake of clarity, only 0 o (not included in shear behavior), 40 o and 60 o are displayed.

Figure 12

 12 Figure 12 Comparison of the cell axis rotation at every loading angle from both experiments and simulations

Figure 13

 13 Figure 13 Normal and shear behaviors of honeycomb under dynamic combined shear-compression

Figure 14

 14 Figure 14 Comparison between dynamic and quasi-static normal and shear pressure/crush curves at loading angle of 40 o

Figure 15

 15 Figure 15 Comparison between dynamic and quasi-static normal and shear pressure vs. loading angle

Figure 16

 16 Figure 16 Crushing envelopes in normal strength vs. shear strength plane

Table 1

 1 Bilinear material parameters

Acknowledgement

The authors would like to thank 111 project of China(contract No.1307050) for funding the cooperation between NPU and LMT. B. Hou and Y. L. Li would also like to thank the supports ff the National Science Foundation of China(contract No.10932008).


Pascal Molin

Multi-precision computation of the complex error function

We give a very simple algorithm to compute the error and complementary error functions of complex argument to any given accuracy.

Introduction

With standard normalization, the error function is defined as erf(x) = 2 √ π

x 0 e -t 2 dt and the complementary error function as

erfc(x) = 2 √ π ∞ x e -t 2 dt.
Thanks to the relation erf(x) = 1erfc(x) and the symmetry erfc(-x) = 2erfc(x), it is sufficient to be able to compute erfc on the right half complex plane.

Many formula exist for this calculation, obtained either via integration or with asymptotic expansions near zero or infinity. In this note, we study an integration-based formula which is simple, efficient and easy to apply rigorously.

Although this formula is not new (we found in particular [MR71]), it seems that it has not been introduced in numerical systems so far. J. Weideman [Wei94] justifies this by the difficulties to find a correct stepsize, and some numerical instabilities. The rigorous treatment we give here aims at removing these prejudices. We also introduce a trick which accelerates the multiprecision computations.

Motivation, and a smooth integral formula

In [START_REF] Molin | Intégration numérique par la méthode double-exponentielle[END_REF], we showed that the trapezoidal scheme yields a very fast method for the computation of integrals of the form

R F(u)e -Q(u 2 ) du (1)
where

F(x) ∈ R(x) is a rational fraction with no pole on R and Q(x) ∈ R + [x] is a polynomial.
If we focus on the simplest case

f (α) = R e -u 2 u -α du, α ∈ R, (2)
we have the following formula, which enables us to compute the complementary error function to any accuracy.

Proposition 1.1

For all x such that Re(x) > 0,

erfc(x) = -ie -x 2 π f (ix). (3)
Proof : f is analytic outside R, derivating there under the integral and integrating by parts yields the differential equation

f ′ (α) + 2α f (α) + 2 √ π = 0. Writing f (α) = λ(α)e -α 2 , λ satisfies λ ′ (α) = -2 √ πe α 2 , so that for Im(α) > 0 we have λ(α) = K -2 √ π α 0 e u 2 du with K ∈ C. The limit lim x→∞ f (ix)e -x 2 = 0 determines K to be 2 √ π i∞ 0 e u 2 du. As a consequence, λ(ix) = 2 √ π i∞ ix e -u 2 du = iπ erfc x, hence the result.
More precisely, the general theory ensures that f (α) can be evaluated to a precision of p binary digits with a number of points n ∼ p log 2 π . The results of [START_REF] Molin | Intégration numérique par la méthode double-exponentielle[END_REF] also assert that this formula is optimal amongst integrationbased formulas, which means that any other integration formula will need a larger number of points (for theoretic reasons due to the uncertainty principle).

Computation

We adopt the following notation for the numerical approximation to p binary digits of absolute precision:

Definition 2.1 For a, b ∈ C, a = p b ⇔ |a -b| 2 -p . (4)
Moreover, the consideration of a residue term leads to the following:

Definition 2.2
Let x be a complex number with Re(x) 0, and h > 0, and define δ(x, h) to be

δ(x, h) = 1 if Re(x) < π h ; (5) = 0 otherwise. (6)
The numerical formula is then Theorem 2.3

Let x be a complex number such that Re(x) 1, and suppose p 2. Then for all h > 0 and n ∈ N such that h

π √ arcsinh(2 p √ π)+2
and nh p log 2 we have

erfc(x) = p e -x 2 π h x + 2hx n ∑ k=1 e -(kh) 2 x 2 + (kh) 2 - 2δ(x + 1, h) e 2πx h -1 . ( 7
)
We give a complete proof in the next paragraphs.

Trapezoidal formula

For α ∈ C such that Im(α) > 0 we define g(x) = e -x 2 x-α . The trapezoidal formula for f = R g reads

Lemma 2.4

Let α ∈ C such that Im(α) > 0. For all h > 0 and n 1, we have

f (α) = h n ∑ k=-n e -(kh) 2 kh -α + E t (n, h) -E q (h) (8)
where

E q (h) = ∑ k =0 ĝ( k h ) and E t (N, h) = ∑ |k|>N g(kh)
.

Proof : This is the Poisson formula written for g on the lattice hZ, isolating the first Fourier term ĝ(0) = f (α).

Quadrature error

Lemma 2.5

Let α ∈ C such that Im(α) > 0 and τ > 0 such that τ = Im(α), then we have for all

X > 0 ĝ(-X) + δ2iπe -α 2 -2iπXα √ π τ -Im(α) e τ 2 -2πXτ (9) ĝ(X) √ π τ + Im(α) e τ 2 -2πXτ (10)
where δ = 0 if τ < Im(α) and δ = 1 if τ > Im(α).

Proof : We have ĝ(±X) = R e -u 2 ∓2πuX u-α
du, and we shift the contour to R ± iτ thanks to the residue theorem. The first line corresponds to R → R + iτ, which goes though the pole at α if τ > Im(α), and the second line to the shift R → R + iτ. We then bound with the L 1 norm, with R e -(u±iτ) 2 ∓2iπ(u±iτ

)X du = e τ 2 -2πτX R e -u 2 du = √ πe τ 2 -2πτX . Lemma 2.6 Assume X 0 1 π arcsinh(2 p √ π) such that πX 0 ∈]Im(α) -1, Im(α) + 1[, then ∀X X 0 , ∑ k =0 ĝ(kX) + δ(x, 1/X 0 ) 2iπe -α 2 e -2iπαX -1 2 -p . ( 11
)
For all X

1 π arcsinh(2 p √ π) + 2 π , ∑ k =0 ĝ(kX) + δ(Im(α) + 1, 1/X) 2iπe -α 2 e -2iπαX -1 2 -p .
Proof : Fix τ = πX 0 in the preceding lemma; the hypothesis allows to remove the terms τ ± Im(α) . Summing on kX, k = 0 we obtain

∑ k =0 ĝ(kX) + δ 2iπe -α 2 e -2iπαX -1 2 √ π ∑ k>0 e -π 2 (kX 0 ) 2 √ π sinh(π 2 X 2 0 ) (12)
hence the first assertion of the lemma. Now take X as in the second part, and set

X 2 = X -2 π , X 1 = X -1 π . If Im(α) πX 1 , then we use the first part with X 0 = X and δ(Im(α), 1/X) = δ(Im(α), 1/X 1 ) = 1. If Im(α) > πX 1 we use the first part with X X 0 = X 2 and δ(Im(α), 1/X 0 ) = δ(Im(α), 1/X 1 ) = 0. Finally, we use δ(x, 1/X 1 ) = δ(x + 1, 1/X).

Truncation error Lemma 2.7

Assume Im(α) > 1, and let p 2. Then for all n such that nh p log 2, we have

h ∑ |k|>n e -(kh) 2 kh + α 2 -p (13)
Proof : Since Im(α) 1 we have (kh) 2 + α 1, and the left-hand side is bounded above by 2h

∞ n e -(th) 2 dt = 2 erfc(nh). The upper bound erfc(x)

2e -x 2 √ π(x+ x 2 + 4 π 2 )
gives the result (we have √ π(x + x 2 + 4 π 2 ) > 4 for x 2 log 2).

Proof of Theorem 2.3

We combine the preceeding lemmas for α = ix, writing

-ih n ∑ k=-n e -(kh) 2 kh + ix = h x + 2xh n ∑ k=1 e -(kh) 2 x 2 + (kh) 2 The value 1 h = X = √ arcsinh(2 p √ π) π + 2
π ensures the hypothesis of Lemma 2.6 is satisfied, so that

E q (h) + 2iπδ(x, h) e 2πx h -1 2 -p .
The lemma 2.7 bounds E t (n, h) , and this proves the computation of -i π f (ix) to absolute precision 2 π 2 -p < 2 -p under the hypothesis of Theorem 2.3. The fact that we obtain the right relative precision on erfc(x) follows from the bounds 1 2 |x| + 1 erfc(x) e -x 2 1 2 |x| -1 valid for any x ∈ C with |x| 1.

Practical algorithm

The formula of Theorem 2.3 is simple. We decribe here how to evaluate it efficiently for multiprecision values.

First, we write λ = x h to put formula (7) into the form

erfc(x) = p e -x 2 π 1 λ + 2λ n ∑ k=1 (e -h 2 ) k 2 λ 2 + k 2 -δ 2 e 2πλ -1 (14) (with h ∼ π √ p log 2 and n ∼ p log 2 π ).
We then write

U k = e -(kh) 2 , V k = e -(2k+1)h 2 and D k = λ 2 + k 2 , so that the main sum becomes n ∑ k=1 U k D k ,
subject to the recursions

U k+1 = U k V k (15) V k+1 = e -2h 2 V k (16) D k+1 = D k + 2k + 1. ( 17
)

Small integer trick

Finally, thanks to the loose condition we have on h, we can improve the computation if we constraint the factor e -2h 2 to be exactly a small precision rational u/2 v , as soon as h =log √ u/2 v satifies Theorem 2.3. This way, the computation of each term of the sum is reduced to the following significant operations

• a multiprecision multiplication U k V k ; • a small multiplication V k u/2 v ; • a multiprecision division U k /D k ;
and we have Theorem 3.1 (complexity)

The complex error function can be evaluated to p binary digits in complexity

p log 2(1+λ) π M(p) + o(pM(p))
, where M(p) denotes the complexity of a multiplication of size p and λ is the number of multiplications needed to perform a division.

Remark : Of course this complexity can be lowered for negative real part, since there the result approaches 2 and the precision required to compute the difference can be decreased.

Improvement for real argument

The trick above can be improved if we apply it to the denominator D k and simplify the division instead of the multiplication. We write

erfc(x) = p e -x 2 πλ 1 + 2 n ∑ k=1 (e -h 2 ) k 2 1 + k 2 /λ 2 - 2δ e 2πλ -1 (18)
and choose h to have exactly λ -2 = u2 -v , u, v ∈ N 2 with v much smaller than p, so that the computation becomes

erfc(x) = p e -x 2 πλ 1 + 2 v+1 n ∑ k=1 (e -h 2 ) k 2 2 v + k 2 u - 2δ e 2πλ -1 . ( 19
)
In order to compute with this formula, the only condition we have is that v > log 2 (x 2 /h 2 ). For the multiprecision range1 , this condition is very mild.

Remark : This trick assumes x to be real, otherwise there is no reason why both its real and imaginary parts should be exactly representable using small integers.

Extension for small real part

The theorem 2.3 can be extended to Re(x) < 1 if we shift the path of integration. Indeed, the function f defined for Im(α) > 0 in (2) extends analytically to Im(α) > -d for any d > 0 with

f (α) = R-id e -u 2 u -α du. ( 20
)
We then have Theorem 3.2

Let x be a complex number such that Re(x) 0. Then for all h > 0 and

n ∈ N such that h < π √ arcsinh(2 p √ π)+2
and nh (p + 1) log 2, we have

erfc(x) = p e -x 2 +1 π 1 λ + 2 n ∑ k=1 (λ cos(2kh) + k sin(2kh))e -(kh) 2 λ 2 + k 2 - 2δ(x + 1, h) e 2πλ -1 (21)
where λ = x+1 h .

Proof : Using (20), the main Riemann sum for f (α) takes the form (πX-d) 2 so that the value X 1 of Lemma 2.6 is enough; for X < 0 we apply the same bounds as in Lemma 2.6. This gives the result, fixing the value d = 1.

S(α) = h n ∑ k=-n e -(kh) 2 e 2idkh e d 2 kh - α , α = α + id = - he d 2 α + he d 2 n ∑ k=1 e -(kh) 2 2ikh sin(2dkh) + 2α cos(2dkh) (kh) 2 -α2 . If Re(x + d) 1,

Timings

We did a PARI/gp implantation which proves to be quite efficient. In Table 1 we compared its running time on a few inputs to MPFR [FHL + 07] and Maple 14. It turns out that the integration formula is not very interesting in the asymptotic range (near zero or infinity), where the asymptotic expansions of erf and erfc are very accurate. However, it sould be considered in the transition range, when the modulus of the argument comes around the square root of the precision (as in the line 10 4 digits, x = 200). When the precision increase, this range gets wider : at 10 5 digits our formula starts to be better from x = 10.

We also remark that the integration formula gives a running time which depends only on the required precision, and to a minor extent of the nature (real or complex) of the argument, contrary to Maple and MPFR.

  the truncation error is bounded by 2 ∞ nh e -t 2 +d 2 e d 2 E t (n, h), so that for d = 1 it is enough to take nh > (p + 3) log 2. (τ±d) 2 -2πXτ . For X > 0 the denominator is greater than d and this can be bounded by

	The quadrature error, once corrected by the residue
	is bounded by	2πδ(x + d, h)e x 2 e 2π(x+d) h -1 τ±Re(x+d) e √ √ π π d e -	,

We remark that current multiprecision library do not allow x to exceed 10 10 .

Algorithm 1: Computation of the complementary error function Input: x such that Re(x) 1 Input: binary precision p 1 Output: z such that (erfc xz)/z < 2 -p begin choose parameters