Xiang-Yu Li

W Q Chen

Hui-Ying Wang

Three-dimensional General Solutions for Thermoporoelastic Media and Its Application

Keywords: Thermoporoelasticity, Transversely isotropic, General solutions, Potential theory method, Penny-shaped crack

This paper presents a set of 3D general solutions for thermoporoelastic media for the steady-state problem. By introducing two displacement functions, the equations governing the elastic, pressure and temperature fields are simplified. The operator theory and superposition principle are then employed to express all the physical quantities in terms of two functions, one of which satisfies a quasi-Laplace equation and the other satisfies a differential equation of the eighth order. The generalized Almansi's theorem is used to derive the displacements, pressure and temperature in terms of five quasi-harmonic functions for various cases of material characteristic roots. To show its practical significance, an infinite medium containing a penny-shaped crack subjected to mechanical, pressure and temperature loads on the crack surface is given as an example. A potential theory method is employed to solve the problem. One integro-differential equation and two integral equations are derived, which bear the same structures to those reported in literature. For a penny shaped crack subjected to uniformly distributed loads, exact and complete solutions in terms of elementary functions are obtained, which can serve as a benchmark for various kinds of numerical codes and approximate solutions.

Introduction

The theory of poroelasticity was originally developed by [START_REF] Biot | General solutions of three-dimensional consolidation[END_REF] to consider the interaction between two phases composing a continuum. Now poroelasticity has many applications in various fields in science and engineering such as soil consolidation (Yue and Selvadurai, 1996), filtration [START_REF] Barry | Approximate and analytic solutions for deformation of finite porous filters[END_REF] [START_REF] Oomens | A mixture approach to the mechanics of skin[END_REF][START_REF] Cowin | Bone poroelasticity[END_REF].

It is well known that geomaterials are typically anisotropic (Cheng, 1997). In fact, the simple natural deposition of sedimentary rocks exhibit transverse isotropy and modeling anisotropic materials as an equivalent isotropic one can yield unexpected and erroneous results which has been demonstrated by analytical and numerical solutions [START_REF] Abousleiman | Mandel's problem revisited[END_REF][START_REF] Cui | Finite element analysis of anisotropic poroelasticity: A generalized Mandel's problem and an inclined borehole problem[END_REF]. Hence, it is natural to extend the isotropic theory to an anisotropic one, which was initially achieved by [START_REF] Biot | Theory of elasticity and consolidation for a porous anisotropic solid[END_REF], based on the formality of the generalized Hooke's law. To consider the effect of temperature, which arises in deep drilling and nuclear storage facilities, thermoporoelasticity model was founded [START_REF] Coussy | Poromechanics[END_REF][START_REF] Abousleiman | poro-thermo-elasticity in transversely isotropic porous materials[END_REF]. Accordingly, analytical solutions of fundamental problems for thermoporoelastic media with material anisotropy play an import role in obtaining a deep insight in the interaction between the phases involved.

Potential theory is the tool for studying problems governed by Laplace equations s in several areas of physics [START_REF] Kellogg | Foundations of potential theory[END_REF]. [START_REF] Fabrikant | Applications of potential theory in mechanics: A selection of new results[END_REF][START_REF] Fabrikant | Mixed boundary value problem of potential theory and their applications in engineering[END_REF]) developed a so called potential theory method by creatively representing the inverse of the distance of two points in Euclidean space in term of an integral (Chen and Ding, 2004a). Based on the general solutions proposed by [START_REF] Elliott | Three-dimensional stress distributions in aeolotropic hexagonal crystals[END_REF] for pure elasticity, [START_REF] Fabrikant | Applications of potential theory in mechanics: A selection of new results[END_REF][START_REF] Fabrikant | Mixed boundary value problem of potential theory and their applications in engineering[END_REF] then performed 3D exact analyses of punch and crack problems by employing the method. Chen and his co-authors extended Fabrikant's potential theory method to solve mixed boundary value problems mainly stemming from three-dimensional crack or contact problems of multi-field coupled media, for example, piezoelasticity [START_REF] Chen | Complete and exact solutions of a penny-shaped crack in a piezoelectric solid: antisymmetric shear loadings[END_REF][START_REF] Chen | On piezoelastic contact problem for a smooth punch[END_REF], thermoelasticity (Chen et al., 2004b), magneto-electro-thermo-elasticity (Chen, 2004c). It is noted that all the works by Chen and his coauthors were also based on the corresponding general solutions (Chen and Ding, 2004a). To the authors' best knowledge, no parallel work in the domain of thermoporoelasticity has been reported yet in literature, mainly due to the lack of the general solution which should be expressed in terms of quasi-harmonic functions as that for purely elastic materials.

The purpose of the present paper is to construct a 3D thermoporoelastic general 3 solution expressed by quasi-harmonic functions only. Two displacement functions are introduced in order for simplifying the basic 3D equations for transversely isotropic thermoporoelastic media in a steady-state. Then, a general solution expressed by two functions is derived by employing the operator theory and the superposition principle.

One function involved in the solution satisfies a quasi-Laplace equation and the other satisfies an eight-order partial differential equation. The general solution is further expressed in terms of five quasi-harmonic functions by virtue of the generalized Almansi's theorem. Four possible cases of material characteristic roots are considered.

To show its practical significance, an infinite medium containing a plane crack subjected to symmetric mechanical, pore pressure and temperature loads is considered for instance. The problem can be transformed to the corresponding mixed boundary-value problem of a half space and thus the potential theory method developed by Fabrikant can be utilized. One integro-differential equation and two integral equations are derived, which bear the same structures as those reported in the literature. For a penny shaped crack subjected to uniformly distributed loads, exact and complete solutions in terms of elementary functions are obtained. The present solution is compared with the results from the literature and good agreement is achieved. The solution presented in this paper can serve as a benchmark for various kinds of numerical codes and approximate analyses.

Governing Equations

In Cartesian coordinates ) , , ( z y x
, the constitutive equations for transversely isotropic thermoporoelastic media with the z-axis coincident with the axis of rotational material symmetry are listed as follows [START_REF] Coussy | Poromechanics[END_REF][START_REF] Abousleiman | poro-thermo-elasticity in transversely isotropic porous materials[END_REF])

M A N U S C R I P T A C C E P T E D
u v w w v c c c c P T c x y z y z u v w w u c c c c P T c x y z x z u v w u v c c P T c x y z y x σ α β σ σ α β σ σ α β σ ⎛ ⎞ ∂ ∂ ∂ ∂ ∂ = + - + - - = + ⎜ ⎟ ∂ ∂ ∂ ∂ ∂ ⎝ ⎠ ∂ ∂ ∂ ∂ ∂ ⎛ ⎞ = - + + - - = + ⎜ ⎟ ∂ ∂ ∂ ∂ ∂ ⎝ ⎠ ⎛ ⎞ ⎛ ⎞ ∂ ∂ ∂ ∂ ∂ = + + - - = + ⎜ ⎟ ⎜ ∂ ∂ ∂ ∂ ∂ ⎝ ⎠ ⎝ , ⎟ ⎠ (1) , 3 1 ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ + ∂ ∂ - ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ∂ ∂ + ∂ ∂ - = T z w y v x u M P m β α α ξ (2) in which ) , ( w v u
and ij σ are components of displacement and stress, respectively; In the following analysis, uncoupled thermoporoelastic theory is adopted by assuming that the thermoporoelastic loading varies slowly in time and the rate of both fluid mass content and entropy vanishes. Consequently, if the medium is in a steady-state, the pressure and temperature fields are constant in time and governed respectively by the following two quasi-Laplace equations , 0

P
2 2 33 11 = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ∂ ∂ + Δ P z κ κ (3) , 0 2 2 33 11 = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ∂ ∂ + Δ T z λ λ (4)
v w P T c c c u c c c c x y x z x x x y z u w P T c c c c c v c c x y y z y y x y z u v c c c c z x y z α β α β ⎛ ⎞ ∂ ∂ ∂ ∂ ∂ ∂ ∂ + + + + + + - - = ⎜ ⎟ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ⎝ ⎠ ⎛ ⎞ ∂ ∂ ∂ ∂ ∂ ∂ ∂ + + + + + + - - = ⎜ ⎟ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ⎝ ⎠ ⎛ ⎛ ⎞ ∂ ∂ ∂ ∂ + + + Δ + ⎜ ⎟ ∂ ∂ ∂ ∂ ⎝ ⎠ 3 3 0, P T w z z α β ⎞ ∂ ∂ - - = ⎜ ⎟ ∂ ∂ ⎝ ⎠ (5)
where 66 11 12 2c c c -= .

In the steady-state, Eqns. (1-4) indicate that the fluid content ξ is not an independent field variable and can be determined once the pressure, temperature and elastic fields have been obtained.

Static thermoporoelastic general solution

Introducing the following two displacement functions 0 Ψ and H as ,

0 x H y u ∂ ∂ - ∂ Ψ ∂ = , 0 y H x v ∂ ∂ - ∂ Ψ ∂ - = (6) Eqns. (3)-(5) can be transformed into , 0 0 2 2 44 66 = Ψ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ∂ ∂ + Δ z c c (7) 0 0 , 0 0 H w P T ⎧ ⎫ ⎧ ⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ = ⎨ ⎬ ⎨ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ⎭ ⎩ ⎭ M (8)
where M is the following operator matrix: .

( ) ( ) 2
0 0 0 0 0 0 c c c c z z c c c c z z z z z z α β α β κ κ λ λ ⎡ ⎤ ∂ ∂ Δ + - + ⎢ ⎥ ∂ ∂ ⎢ ⎥ ∂ ∂ ∂ ∂ ⎢ ⎥ - + Δ Δ+ - - ⎢ ⎥ ∂ ∂ ∂ ∂ = ⎢ ⎥ ∂ ⎢ ⎥ Δ + ⎢ ⎥ ∂ ⎢ ⎥ ∂ ⎢ ⎥ Δ + ⎢ ⎥ ∂ ⎣ ⎦ M (9)
The "determinant" of M reads The above derivation procedure is similar to that given by Chen et al. (2004b), where the equation of heat conduction and the other two equations governing the displacement functions are combined, instead of solving the temperature field independently and a prior as shown in the traditional general solutions with a particular part [START_REF] Ashida | General solution technique for transient thermoelasticity of transversely isotropic solids in cylindrical coordinates[END_REF][START_REF] Ding | Transversely Isotropic Elasticity[END_REF]. Here, the idea is employed to seek the 3D general solutions for thermoporoelastic media, although the pressure and temperature field can be solved independently.

M A N U S C R I P T A C C E P T E D ARTICLE IN PRESS X.-Y. LI, et al / European Journal of Mechanics A/Solids 6 4 2 2 2 2 0 0 0 1 1 3 3 1 1 3 3 4 2 2 2 , a b c z z z z κ κ λ λ ⎛ ⎞ ⎛ ⎞ ⎛ ⎞ ∂ ∂ ∂ ∂ = + Δ+ Δ Δ+ Δ+ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ∂ ∂ ∂ ∂ ⎝ ⎠ ⎝ ⎠ ⎝ ⎠ M (10)
According to the operator theory, the following general solutions can be obtained

1 2 3 4 , , , , ( 1, 2, 3, 4)
,

i i i i H C F w C F P C F T C F i = = = = = (12)
where ij C are the algebraic cofactors of M , and F satisfies , which are actually the same to those for pure elasticity [START_REF] Elliott | Three-dimensional stress distributions in aeolotropic hexagonal crystals[END_REF][START_REF] Ding | Transversely Isotropic Elasticity[END_REF];

2 2 2 2 2 2 2 2 1 2 3 4 0, F z z z z ⎛ ⎞ ⎛ ⎞⎛ ⎞ ⎛ ⎞ ∂ ∂ ∂ ∂ Δ + Δ + Δ + Δ + = ⎜ ⎟ ⎜ ⎟⎜ ⎟ ⎜ ⎟ ∂ ∂ ∂ ∂ ⎝ ⎠⎝ ⎠ ⎝ ⎠ ⎝ ⎠ (
3 = i corresponds to the general solution 1 X (say) with 0 = T
, which is identical to that for poroelasticity; Taking 4 = i corresponds to the general solution 2 X (say) with 0 = P , which is equivalent to the one given by Chen et al. (2004b) for thermoelasticity.

Due to the linear nature of the thermoporoelastic theory adopted in this paper, superposing 1 X and 2 X leads to

M A N U S C R I P T A C C E P T E D ARTICLE IN PRESS X.-Y. LI, et al / European Journal of Mechanics A/Solids 7 4 2 2 0 1 1 1 4 2 4 2 2 0 1 1 1 4 2 4 2 2 2 2 2 4 2 , , , F u a b c y z z x F v a b c x z z y F w a b c z z z ⎛ ⎞ ∂Ψ ∂ ∂ ∂ = - + Δ+ Δ ⎜ ⎟ ∂ ∂ ∂ ∂ ⎝ ⎠ ⎛ ⎞ ∂Ψ ∂ ∂ ∂ = - - + Δ + Δ ⎜ ⎟ ∂ ∂ ∂ ∂ ⎝ ⎠ ⎛ ⎞ ∂ ∂ ∂ = + Δ+ Δ ⎜ ⎟ ∂ ∂ ∂ ⎝ ⎠ 6 4 2 2 3 3 3 3 3 6 4 2 , P a b c d F z z z ⎛ ⎞ ∂ ∂ ∂ = + Δ+ Δ + Δ ⎜ ⎟ ∂ ∂ ∂ ⎝ ⎠ 6 4 2 2 3 4 4 4 4 6 4 2 , T a b c d F z z z ⎛ ⎞ ∂ ∂ ∂ = + Δ+ Δ + Δ ⎜ ⎟ ∂ ∂ ∂ ⎝ ⎠ (15)
where

( )( ) ( ) ( )( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )( ) (
) )(

1 1 3 4 4 3 3 3 3 3 3 3 3 3 3 1 3 3 1 1 13 44 11 3 11 3 44 33 1 33 1 33 11 1 11 1 1 4 4 1 1 1 1 1 1 2 4 4 3 3 3 3 3 3 2 11 33 3 33 3 44 11 3 11 3 13 44 33 1 33 1 2 1 1 1 1 3 1 1 3 1 , , , , , ( a
c c c b c c c c c c a c b c c c c c c c κ β λ α κ β λ α κ β λ α κ β λ α κ β λ α κ β λ α κ β λ α κ β λ α κ β λ α κ β λ α κ β λ α = + + - + = + + - + - + = - + = + = + + + - + + = + -3 4 4 1 1 1 1 1 1 3 0 33
a b a b c b c d c a a b a b c b c d c κ β λ α λ λ λ λ λ λ κ κ κ κ κ κ + + = = + = + = = = + = + = ), , , , , , , , . c a
Taking advantage of the generalized Almansi's theorem [START_REF] Ding | General solutions for coupled equations for piezoelectric media[END_REF], F can be expressed by four quasi-harmonic functions

1 2 3 4 1 2 3 4 1 2 3 4 1 2 1 2 3 4 1 2 3 4 2 3 1 2 3 4 1 2 3 4 , for distinct , (Case 1) , s s , (Case 2) , s s , ( Case 3)
, s s , (C i F F F F s F F F zF s s s F F F zF z F s s F zF z F z F s s + + + + + + ≠ ≠ = ≠ = + + + ≠ = = + + + = = = ase 4) ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ (17)
where i F respectively satisfy 2 2 0, ( 1, 2, 3, 4).

i i F i z ⎛ ⎞ ∂ Δ + = = ⎜ ⎟ ∂ ⎝ ⎠ (18)
For Case 1, the general solution in Eqn. (15) takes on the simplest form as

5 5 4 4 0 0 1 1 4 4 1 1 5 6 6 4 3 4 2 3 3 4 4 5 6 6 1 3 4 , , , , , i i i i i i i i i i i i i F F u v y x z x y z F F F w s P T z z z η η η η η = = = ∂Ψ ∂ ∂Ψ ∂ = - =- - ∂ ∂∂ ∂ ∂∂ ∂ ∂ ∂ = = = ∂ ∂ ∂ ∑ ∑ ∑ (19) with M A N U S C R I P T A C C E P T E D ARTICLE IN PRESS X.-Y. LI, et al / European Journal of Mechanics A/Solids 8 ( ) ( ) 4 2 4 2 2 33 0 3 0 3 0 33 3 11 4 2 2 44 0 3 0 3 0 33 4 11 ( 1,2), ( ) , ( ) . i i ij j j j a s b s c j a s b s c s a s b s c s η η λ λ η κ κ = - + = = - + - = - + - (20)
Letting ( )

4 1 4 , 1, 2,3, 4 i i i i F i z η ∂ Ψ = = ∂ (21)
Eqn. ( 19) can be rewritten as

4 4 4 0 0 1 1 1 1 2 2 4 4 2 3 2 2 1 1 , , ,
, ,

i i i i i i i i i i i i i i i i u v w y x x y z P T z z μ μ μ = = = = = ∂Ψ ∂Ψ ∂Ψ ∂Ψ ∂Ψ = - =- - = ∂ ∂ ∂ ∂ ∂ ∂ Ψ ∂ Ψ = = ∂ ∂ ∑ ∑ ∑ ∑ ∑ (22)
where

1 2 1 ( 1 ) 2 3 3 4 1 , ( )( ) , i i i i i j i j j j i i i s μ λ λ μ λ δ δ δ δ λ + = = + + ( 2,3, 1 4 j i = = -) and ij δ
is the Kronecker delta, and 2 2 0, ( 0, 1, 2, 3, 4),

i i i z ⎛ ⎞ ∂ Δ + Ψ = = ⎜ ⎟ ∂ ⎝ ⎠ (23)
i i i i u u ρ φ ρ φ ρ ρ ρ φ = = ∂Ψ ∂Ψ ∂Ψ ∂Ψ = - =- - ∂ ∂ ∂ ∂ ∑ ∑ (24)
The general solutions for the other three cases are given in Appendix B for the theoretical completeness.

Once the general solutions are available, it is very convenient to conduct threedimensional analyses of the mixed boundary-value problems associated with cracks and punches by virtue of the potential theory method proposed by [START_REF] Fabrikant | Applications of potential theory in mechanics: A selection of new results[END_REF][START_REF] Fabrikant | Mixed boundary value problem of potential theory and their applications in engineering[END_REF] as shown in the parallel works by Chen and his coauthors (2000Chen and his coauthors ( , 2004bChen and his coauthors ( , 2004c).

In the following section, a flat crack problem is considered as an illustrative example to show the practical significance of the general solution for Case 1.

Generalized potential theory method for thermoporoelastic crack problem

To generalize the potential theory method suggested by [START_REF] Fabrikant | Applications of potential theory in mechanics: A selection of new results[END_REF]

U u iv i i σ σ σ σ σ σ σ τ σ σ = + = + = -+ = + (25)
where 1 -= i

. Then from Eqns (1) and ( 22), the expressions of the thermoporoelastic field take on a more compact form

2 2 4 4 4 0 0 0 1 1 2 2 2 1 1 1 4 4 2 0 2 6 6 0 3 4 40 1 1 0 , , , 2 , , i z j j j j j i i j j z j j j j U i z z c i i c s z z σ γ σ γ σ τ γ = = = = = ⎛ ⎞ ∂ Ψ ∂ Ψ = -Λ Ψ + Ψ = = ⎜ ⎟ ∂ ∂ ⎝ ⎠ ⎛ ⎞ ∂Ψ ⎛ ⎞ ∂Ψ = - Λ Ψ + Ψ = Λ - ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ∂ ∂ ⎝ ⎠ ⎝ ⎠ ∑ ∑ ∑ ∑ ∑ (26)
where

y i x ∂ ∂ + ∂ ∂ = Λ
, and

1 1 3 3 3 1 3 3 2 1 1 1 6 1 3 1 1 2 1 3 3 4 4 1 , 2 , ( ), 1 4. j jj j j jj j j j j j c c s c c c s c s j γ μ βμ γ μ αμ βμ γ μ ⎡ ⎤ = + - = -+ - - ⎣ ⎦ = - = - (27)
It is note that the identity 3 1 ,( 1 4)

π τ ∈ = Θ =Θ =Θ ∂ ∂ ∈ = = = ∂ ∂ ∈ = (28)
where the factor π 2 is introduced for the algebraic convenience. In order to extend the potential theory method [START_REF] Fabrikant | Applications of potential theory in mechanics: A selection of new results[END_REF] to thermoporoelasticity, assume that ( )

3 0 1 0, , , , ( 1 4), i i j j i j e
x y z i

= Ψ = Ψ = Ω = - ∑ (29)
where ij e are constants to be determined, and

[ ] { } 1 0 1 0 0 0 ( ) ( , , ) , ( , ) ( , , ) ( ) ln ( , ) ( , ) , ( 2,3), S j j S N x y z dS R M N x y z N z R M N z R M N dS j ω ω Ω = Ω = + - = ∫∫ ∫∫ (30)
where the kernels of the potentials 1 Using the well known property of a simple layer potential, one can arrive at

3 1 0 0 3 1 0 1 3 2 0 2 0 3 3 3 0 3 0 3 ( , ) : 0, 0, ( 2, 3) ( , ) : 2 2 ( , ,0), ( , , ) ( , ) : 2 2 , ( , , ) ( , ) : 2 2 . j z z z z z z z x y S j z z x y S w x y z P x y z x y S z z T x y z x y S z z πω π πω π πω π = = = = = = = ∂ Ω ∂Ω ∈ = = = ∂ ∂ ∂Ω ∈ =- =- ∂ ∂ Ω ∂ ∈ =- =- ∂ ∂ ∂ Ω ∂ ∈ =- =- ∂ ∂ (32)
Eqns. ( 22), (32) and the second condition in Eqn. , , ( 2,3; 1, 2,3). 2 2

j j k i ij ik ij i i i e es k j δ δ μ μ π π = = = - = - = = ∑ ∑ (33)
Thus combining Eqns. ( 31) and ( 33) yields

γ γ γ γ δ μ μ μ μ δ μ μ μ μ π δ μ μ μ μ - ⎧ ⎫ ⎧ ⎫ ⎡ ⎤ ⎪ ⎪ ⎪ ⎪ ⎢ ⎥ ⎪ ⎪ ⎪ ⎪ ⎢ ⎥ = - = ⎨ ⎬ ⎨ ⎬ ⎢ ⎥ ⎪ ⎪ ⎪ ⎪ ⎢ ⎥ ⎪ ⎪ ⎪ ⎪ ⎣ ⎦ ⎩ ⎭ ⎩ ⎭ (34)
To satisfy the first condition in Eqn. ( 28), we have gives arise to

3 0 1 0 11 1 1 2 0 0 ( ) ( ) 2 ( ), ( , ) ( , ) j j j S S N N m d S m d S N R N N R N N ω ω π = Δ - = -Θ ∑ ∫∫ ∫∫ (35) 0 1 0 ( ) 2 ( ), ( 2,3), ( , ) j j j S N dS s N j R N N ω π + = - Θ = ∫∫ (
1 0 1 0 ( ) 2 ( ), ( , ) S N dS N R N N ω π Δ = -Θ ∫∫ (37)
where [ ]

1 1 12 3 2 13 4 3 11 ( ) ( ) ( ) ( ) , N N m s N m s N m Θ = Θ + Θ + Θ (38)
which can be viewed as the generalized mechanical load. Once ) , ( y x i ω is available, the thermoporoelastic field can be obtained by differentiating the potential function It is seen that the pressure and thermal loadings can effect the elastic field (see Eqns. 37 and 38), while the mechanical loadings can change neither the thermal nor pressure field (see Eqn. 36); furthermore, no mutual influences exist between the thermal and pressure fields. This is due to the uncoupled theory employed herein. ), which provides a way to verify the solutions in this paper by comparing with those in literature.

It is interesting to note that the structure of integral equations in Eqn. ( 36) is the same as the one for punch problems in elasticity while the structure of the integro-differential equation in Eqn. ( 37) is identical to that for crack problems in elasticity [START_REF] Fabrikant | Applications of potential theory in mechanics: A selection of new results[END_REF]Chen, 2004b). Due to this similarity, the splendid results

given by [START_REF] Fabrikant | Applications of potential theory in mechanics: A selection of new results[END_REF][START_REF] Fabrikant | Mixed boundary value problem of potential theory and their applications in engineering[END_REF] can be directly employed to seek the solutions to the governing equations in some particular cases. This will be shown in the next section, where the cylindrical coordinate system ) , , ( z φ ρ with the origin coincident to the center of the crack is alternatively used for simplicity.

A penny-shaped crack subjected to uniformly distributed loads

Consider a penny-shaped crack subjected to uniformly distributed mechanical, pressure and temperature loads

0 0 0 0 1 1 2 2 3 3 1 1 ( , ) , ( , ) , ( , ) , ( , )
, x y x y x y x y Θ =Θ Θ =Θ Θ =Θ Θ = Θ (39)
where the quantities with the superscript 0 are prescribed constants. Under this condition, the exact solutions to Eqns. ( 36) and ( 37) can be obtained by using the results of [START_REF] Fabrikant | Applications of potential theory in mechanics: A selection of new results[END_REF] as

0 0 1 2 2 1 1 2 2 2 2 ( , ) , ( , ) , ( 2,3)
,

j j j s a j a ω ρ φ ρ ω ρ φ π π ρ + Θ Θ = - =- = - (40)
where a is the radius of the crack. Inserting Eqn. (39) into Eqn. (30) leads to ( ) ( )

2 2 0 2 2 2 1 2 2 1 1 1 2 2 2 0 2 2 1 1 2 2 2 2 2 2 2 1 2 2 2 2 3 ( , , ) 2 2 sin , ( , , ) 2 sin 2 2 3 2 ln , ( 2,3) j j j a l a z a z l a l a a z s a z l a l l a az l l a j a ρ φ ρ ρ ρ φ - - + ⎡ ⎤ ⎛ ⎞ - Ω =Θ + - - - ⎢ ⎥ ⎜ ⎟ ⎝ ⎠ ⎣ ⎦ ⎡ ⎛ ⎞ ⎛ ⎞ Ω =- Θ + - ⎢ ⎜ ⎟ ⎜ ⎟ ⎝ ⎠ ⎝ ⎠ ⎣ ⎤ - - -+ + - = ⎥ ⎦ (41) where M A N U S C R I P T A C C E P T E D ARTICLE IN PRESS X.-Y. LI, et al / European Journal of Mechanics A/Solids ( ) ( ) 2 2 2 2 1 1 2 2 2 2 2 2 1 ( , , ) ( ) ( ) , 2 1 ( , , ) ( ) ( ) . 2 l l a z a z a z l l a z a z a z ρ ρ ρ ρ ρ ρ = = + + - - + = = + + + - + (42)
If there are no mechanical and pressure loads, 3 ω and 3 Ω vanish identically, and expressions for j ω and ) 2 , 1 ( = Ω j j are exactly the same as these given by Chen (2004b). From Eqns. ( 41)-( 42), the thermoporoelastic field is independent of the angular coordinate φ because of the axisymmetry pertinent to the problem.

At this stage, one can derive the whole thermoporoelastic field through

differentiation 2 2 4 2 0 1 1 1 2 1 2 2 2 2 3 4 1 0 2 2 1 1 2 2 1 2 2 s i n 4 s i n , 2 2 j i j j j j i j k k j k j j k j j a l a a U e e l l l e a s e l a a az a l φ φ ρ ρ ρ - = - + = = ⎡ ⎤ - ⎛ ⎞ ⎢ ⎥ = Θ - ⎜ ⎟ ⎜ ⎟ ⎢ ⎥ ⎝ ⎠ ⎣ ⎦ ⎡ ⎤ ⎛ ⎞ ⎛ ⎞ - Θ - - - + ⎢ ⎥ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎢ ⎥ ⎝ ⎠ ⎝ ⎠ ⎣ ⎦ ∑ ∑ ∑ ( ) 4 0 1 2 2 1 1 1 1 1 2 3 4 0 1 2 2 2 2 1 1 1 2 2 2 1 2 4 s i n 4 s i n l n , j j j j j j k k j j k j j j j k j j a w e z a l l a s e z a l a l l l μ μ ρ - = - + = = ⎡ ⎤ ⎛ ⎞ = Θ - - ⎢ ⎥ ⎜ ⎟ ⎜ ⎟ ⎢ ⎥ ⎝ ⎠ ⎣ ⎦ ⎡ ⎤ ⎛ ⎞ - Θ - -+ + - ⎢ ⎥ ⎜ ⎟ ⎜ ⎟ ⎢ ⎥ ⎝ ⎠ ⎣ ⎦ ∑ ∑ ∑ 2 2 4 3 4 2 0 1 0 1 1 1 1 1 1 2 2 1 2 1 2 2 1 2 2 2 4 3 2 0 1 0 1 1 1 2 1 1 2 2 2 1 2 2 2 1 2 4 sin 4 sin , 4 s i n 4 s i n j z j j k k j j k j k j j j j j j j j k k j jk j k j j j j a l a a a e s e l l l l a l a a a e s e l l l l σ γ γ σ γ γ - - + = = = - - + = = ⎡ ⎤ - ⎛ ⎞ ⎛ ⎞ ⎢ ⎥ = Θ - - Θ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ - ⎢ ⎥ ⎝ ⎠ ⎝ ⎠ ⎣ ⎦ ⎡ ⎤ - ⎛ ⎞ ⎛ ⎞ ⎢ ⎥ = Θ - - Θ ⎜ ⎟ ⎜ ⎜ ⎟ ⎜ - ⎢ ⎥ ⎝ ⎠ ⎝ ⎠ ⎣ ⎦ ∑ ∑ ∑ ∑ ∑ 4 1 2 2 2 4 1 1 2 0 2 2 6 6 1 2 2 2 1 2 2 1 2 3 4 0 2 2 2 2 66 1 2 1 2 2 1 2 2 4 3 4 3 1 1 0 2 0 2 2 1 1 3 1 2 2 2 1 2 1 2 2 1 , 8 ( ) 4 ( 2 ) , 4 4 ( ) ( ) j j j j i j j j j i k k j k j j j j k j i j j j i z k k j j k j j k j j
+ = = = ⎟ ⎟ - = - Θ - + Θ -+ -- - = Θ + Θ -- - ∑ ∑ ∑ ∑ ∑ ∑ , ∑ 0 1 0 1 2 3 23 24 2 2 sin , sin , a a P T l l π π - - ⎛ ⎞ ⎛ ⎞ = Θ = Θ ⎜ ⎟ ⎜ ⎟ ⎝ ⎠ ⎝ ⎠ (43)
where ) , , (

j ij ij z a l l ρ = , ( 4 1 , 2 , 1 - = = j i
) defined in Eqn. ( 42). As expected, the ( ,) .

z m a a S a S ρ φ ρ σ ρ φ ρ ρ φ - ⎧ ⎛ ⎞ Θ Θ - ∈ ⎪ ⎜ ⎟ = - ⎝ ⎠ ⎨ ⎪ Θ ∈ ⎩ 0 1 2 0 2 0 1 3 0 3 2 sin , (
a S P S a S T S ρ φ ρ φ π ρ ρ φ ρ φ ρ φ π ρ ρ φ - - ⎧ ⎛ ⎞ Θ ∈ ⎪ ⎜ ⎟ = ⎨ ⎝ ⎠ ⎪ Θ ∈ ⎩ ⎧ ⎛ ⎞ Θ ∈ ⎪ ⎜ ⎟ = ⎨ ⎝ ⎠ ⎪ Θ ∈ ⎩ (45)
It is evident that the boundary conditions have been satisfied. It is also noted that the normal stress z σ on S is the sum of two terms which are singular and non-singular parts on the crack tip, respectively. This is similar to the crack problem in elasticity [START_REF] Fabrikant | Applications of potential theory in mechanics: A selection of new results[END_REF].

If the stress intensity factor is defined as 0 lim 2 ( ) ,

I z z a K a ρ π ρ σ = → = - (46)
then from Eqn. ( 45), one can obtain ( )

0 0 0 0 11 1 1 3 12 2 4 13 3 4 4 . I K m a a sm s m π π = - Θ = - Θ + Θ + Θ (47)
It is interesting to note that if the pressure and temperature loads are absent, the stress intensity factor is independent of the material property and is consistent with the result given by [START_REF] Fabrikant | Applications of potential theory in mechanics: A selection of new results[END_REF].

If the mechanical and pressure loads are missing and the temperature is solely presented, the forms of the displacements and stresses shown in Eqn ( 43) are the same as that in Chen (2004b) with an exception of the coefficients. However, numerical calculation reveals that the results presented in this section agrees well with those shown in the previous work. Furthermore, the pressure field in Eqn. ( 43) can be verified by solving the mixed boundary-value problems of a half space. These two aspects will be presented in the next section. ,,) ,

a P z d d R π ω ρ φ ρ φ ρ ρφ ′ = ′ ∫ ∫ (48)
where

3 2 2 2 2 0 0 0 0 3 3 1 , 2 cos( ) . 2 z P R z z ω ρ ρ ρ ρ φφ π = ∂ ′ ′ = = + - - + ∂ (49)
It is noted that the boundary condition ( , ) 0, ( , )

P S z ρ φ ρ φ ∂ = ∈ ∂ is automatically
meted. The boundary condition on S is also satisfied by letting

2 2 0 0 0 0 0 2 2 2 0 0 0 0 0 ( , ) ( , ). 2 cos( )
a d d π ω ρ φ ρ ρ φ ρ φ ρ ρ ρρ φ φ ′ =Θ + - - ∫ ∫ (50)
By virtue of the results given by [START_REF] Fabrikant | Applications of potential theory in mechanics: A selection of new results[END_REF], the solution to Eqn. ( 49) can be written out as

( ) ( ) ( ) 2 2 2 2 2 2 2 2 0 1 1 ( , )
, ,

a x d xdx d ydy L L Ly y d x d x x x y ρ ω ρ φ ρ φ π ρ ρ ρ ⎛ ⎞ ′ = - Θ ⎜ ⎟ ⎝ ⎠ - - ∫ ∫ (51)
where ) (• L is an operator defined as:

2 2 0 0 2 0 0 1 1 ( ) ( ) ( ) , (0 1). 2 1 2 cos( ) k L k f f d k k k π φ φφ π φ φ - = ≤ < + - - ∫ (52)
Substituting Eqn. (52) into Eqn. ( 48) yields ( )

2 1 3 2 0 0 0 0 0 3 2 0 0 1 ( , , ) tan ( , ) , a z R h P z d d h R R π ρ φ ρ φ ρ ρ φ π - ′ ⎡ ⎤ ⎛ ⎞ = + Θ ⎜ ⎟ ⎢ ⎥ ′ ′ ⎝ ⎠ ⎣ ⎦ ∫ ∫ (53)
π -⎛ ⎞ = Θ ⎜ ⎟ ⎝ ⎠ (55)
which is identical to the expression in Eqn. ( 43). The temperature field can be verified in a similar manner.

To validate the solutions presented in the last section, numerical calculation is performed for a penny-shaped crack subjected to a constant thermal load

C 0 0 3 100 = Θ
only. The material constants are listed in Table 1 for a deep-water Gulf shale of Mexico [START_REF] Kanj | Porothermoelastic analyses of anisotropic hollow cylinders with applications[END_REF], where the physical meanings are given in Appendix A. Under this situation, the solution in the preceding section should be identical to those given by Chen et al. (2004b). Comparison is made in Table 2 and good agreement is achieved. This further validates the thermoporoelastic fields in Eqn.

(43).

Tables 12

Concluding remarks

The basic equations in static thermoporoelasticity are simplified by introducing two displacement functions. Then, a set of general solutions in terms of five quasi-harmonic functions are derived by using the operator theory and the superposition principle along with the generalized Almansi's theorem. Four possible cases of material characteristic roots are considered for the theoretical completeness.

No particular solutions associated with the pressure and temperature fields appear in the general solutions, making it possible to extend the potential theory method to thermoporoelasticity.

Based on the general solution, the potential theory method developed by [START_REF] Fabrikant | Applications of potential theory in mechanics: A selection of new results[END_REF] for the mixed boundary-value problems was extended to the transversely isotropic thermoporoelastic plane crack problems. Two potentials with the same mathematical structure are introduced to account for the thermal and pressure effects.

The governing equations have the same structures as these reported in literatures. For the problem of a penny-shaped crack under the actions of uniformly distributed mechanical, pressure and temperature loadings, exact and complete thermoporoelastic fields are expressed explicitly in terms of elementary functions. The present solution is verified by comparing to those given by Chen et al. (2004b). Furthermore, the solutions pertinent to the pore pressure are completely new to literature.

It should be emphasized that the thermoporoelastic field is only valid for the materials with distinct characteristic roots. In the case of multiple roots, the general solutions will takes on different forms as shown in Appendix B. The potential theory method can also be employed if appropriate potentials are chosen. However, as pointed out by [START_REF] Fabrikant | Applications of potential theory in mechanics: A selection of new results[END_REF], the corresponding solutions in the case of multiple roots can be derived directly from the solutions, as presented in this paper, with the aid of L'Hostpital rule.

Although the loadings on the penny-shaped crack in this paper are independent of the angular coordinate, it is not the limitation inherent to the potential theory method.

In fact, the problems of a penny-shaped crack subjected to the point loading can be successfully solved and the Green's functions which are significant to the boundary element method can be derived. This will be reported in another paper.

) ( ) , , (1 )
( 2 ) (1 )( 2 ) (1 ) , , , 2(1 ) 2 2 E E E E E E c c E E E E E E EE E E c c c E E E E E E ν ν ν ν ν ν ν ν ν ν ν ν ν ν ν ν ′ ′ ′ ′ - + = = ′ ′ ′ ′ ′ ′ + - - + - - ′ ′ ′ ′ - = = = ′ ′ ′ ′ ′ ′ ′ + - - - - (A.1)
where E is the drained Young's modulus in the plane of isotropy ( y xplane), E′ the drained Young's modulus perpendicular to that plane ( z -direction), respectively. If the assumption of micro-homogeneity and micro-isotropy of the solid skeleton is adopted, 1 where φ′ is the porosity.

Appendix B

The general solutions for the cases of multiple roots of , ( 1, 2), ,

i i i i F F F F i s s z z z z η η η η ∂ ∂ ∂ ∂ ′ Ψ = = Ψ = + Ψ = ∂ ∂ ∂ ∂ (B.1) where ) 4 1 , ( , 2 4 2 - = - = ′ j i b s a i i i ij η , the general solutions are ( ) ( ) 0 0 1 2 3 3 4 1 2 3 3 4
u z v z y x x y w z P T z z z z z z μ μ μ μ μ μ μ ∂Ψ ∂Ψ ∂ ∂ = - Ψ + Ψ + Ψ + Ψ = - - Ψ + Ψ + Ψ + Ψ ∂ ∂ ∂ ∂ ∂Ψ ∂Ψ ∂Ψ ∂Ψ ∂Ψ ∂Ψ = + + + + Ψ = = ∂ ∂ ∂ ∂ ∂ ∂ (B.
η μ μ η η μ μ η η ′ ′ ⎛ ⎞ = - + ⎜ ⎟ ⎝ ⎠ - + - + = = (B.3) It is noted that ( 0 4) i i Ψ = - satisfies
F F F F s z z z z F F F s z z s z η η η η η η η ∂ ∂ ∂ ∂ ′ ′ ′ Ψ = Ψ = + + ∂ ∂ ∂ ∂ ∂ ∂ ∂ ′ Ψ = + Ψ = ∂ ∂ ∂ (B.4)
u z z y x z v z z x y z w z z z z z z z z P z T z z z μ μ μ μ μ μ μ μ μ ⎛ ⎞ ∂Ψ ∂Ψ ∂ = - Ψ +Ψ + Ψ + ⎜ ⎟ ∂ ∂ ∂ ⎝ ⎠ ⎛ ⎞ ∂Ψ ∂Ψ ∂ = - - Ψ + Ψ + Ψ + ⎜ ⎟ ∂ ∂ ∂ ⎝ ⎠ ⎛ ⎞ ∂Ψ ∂Ψ ∂Ψ ∂ Ψ ∂Ψ = + + + + Ψ + + Ψ ⎜ ⎟ ∂ ∂ ∂ ∂ ∂ ⎝ ⎠ ⎛ ⎞ ∂Ψ ∂ Ψ ∂Ψ = + + = ⎜ ⎟ ∂ ∂ ∂ ⎝ ⎠ 2 3 4 4 35 2 36 2 2 2 2 2 , z z z z μ ⎛ ⎞ ∂Ψ ∂ Ψ ∂Ψ + + ⎜ ⎟ ∂ ∂ ∂ ⎝ ⎠ (B.
η η μ η μ η μ η η η η η η η μ μ μ μ η η η η η η ′′ ′ ′′ ′ ⎡ ⎤ + - - ⎣ ⎦ = ′′ ′ ′ ′ ′ ′ ′ ⎛ ⎞ ⎛ ⎞ = - + + = - + + ⎜ ⎟ ⎜ ⎟ ⎝ ⎠ ⎝ ⎠ (B.6)
In view of Eqns. (B.4) and ( 18), it is also seen that ( 0 4)

i i Ψ = - satisfies Eqn. (23) Case 4 ( 4 3 2 1 s s s s = = = ): Introducing 2 4 3 3 1 2 4 1 1 1 1 1 1 1 1 1 1 1 4 3 2 1 1 1 1 , 3 3 2 , , ,
F F F F s a z z z z F F F F F F s s s z z z s z z s z η η η η η η η η η η ∂ ∂ ∂ ∂ ′ ′ ′ ′ ′ Ψ = + + + ∂ ∂ ∂ ∂ ′′ ′ ∂ ∂ ∂ ∂ ∂ ∂ ′ Ψ = + + Ψ = + Ψ = ∂ ∂ ∂ ∂ ∂ ∂ (B.7) into Eqn. (15) leads to 2 2 3 0 4 4 1 1 2 1 1 2 1 1 2 2 3 0 4 4 1 1 2 1 1 2 1 1 2 3 2 3 1 2 4 4 11 1 1 1 2 3 1 2 1 1 2 3 3 4 4 15 2 1 1 16 3 1 2 1 1 1 , , 2 3 3 u z z z y x z z v z z z x y z z w z z z z z z z z z z z z z μ μ μ ⎛ ⎞ ∂Ψ ∂Ψ ∂ Ψ ∂ = - Ψ + Ψ + + ⎜ ⎟ ∂ ∂ ∂ ∂ ⎝ ⎠ ⎛ ⎞ ∂Ψ ∂Ψ ∂ Ψ ∂ = - - Ψ + Ψ + + ⎜ ⎟ ∂ ∂ ∂ ∂ ⎝ ⎠ ⎛ ⎞ ∂Ψ ∂Ψ ∂ Ψ ∂ Ψ = + + + ⎜ ⎟ ∂ ∂ ∂ ∂ ⎝ ⎠ ⎛ ⎞ ⎛ ∂Ψ ∂ Ψ ∂Ψ + Ψ + + + Ψ + ⎜ ⎟ ⎜ ∂ ∂ ∂ ⎝ ⎝ ⎠ 17 4 2 3 2 3 3 3 3 2 4 25 1 1 26 1 27 4 2 3 2 1 1 2 1 1 2 3 2 3 3 3 3 2 4 35 1 1 36 1 37 4 2 3 2 1 1 2 1 1 , 2 3 3 , 2 3 3 , P z z z z z z z z T z z z z z z z z μ μ μ μ μ μ μ ⎞ + Ψ ⎟ ⎠ ⎛ ⎞ ⎛ ⎞ ∂ Ψ ∂ Ψ ∂Ψ ∂Ψ ∂ Ψ = + + + + + Ψ ⎜ ⎟ ⎜ ⎟ ∂ ∂ ∂ ∂ ∂ ⎝ ⎠ ⎝ ⎠ ⎛ ⎞ ⎛ ⎞ ∂ Ψ ∂ Ψ ∂Ψ ∂Ψ ∂ Ψ = + + + + + Ψ ⎜ ⎟ ⎜ ⎟ ∂ ∂ ∂ ∂ ∂ ⎝ ⎠ ⎝ ⎠ (B.8) where ( ) ( ) ( ) 3

  Y. LI, et al / European Journal of Mechanics A/Solids 2 bio-science

  Y. LI, et al / European Journal of Mechanics A/Solids

  and T are changes in the pore pressure and temperature, respectively; 0 = P and 0 = T correspond to the stress-free state; ξ is variation of fluid content; elastic moduli, Biot's effective stress coefficients and thermal constants, respectively. It is noted that ij c , i α and i β can be expressed in terms of engineering constants such as Young's moduli, Poisson's ratio, etc., which are shown in Appendix A (Kanj and Abousleiman, 2005).

  and thermal conductivity. The equations governing the displacement field can be easily derived by inserting Eqn. (1) into the equilibrium equations in terms of stress components and are shown to be

  of w , P and T keep unchanged as given in Eqn. (22), while the expressions of the radial and circumferential displacements turn out to be

  true, which can be verified by direct substitution. Now consider an infinite transversely isotropic thermoporoelastic body weakened by a flat crack whose surface is parallel to the plane of isotropy. It is further assumed that the radius of the borehole 0 r is much smaller than the crack length 0 L , so that the borehole effect can be ignored. The coordinates are chosen such that the y xplane coincides with the crack surface and the origin is located on the surface of the crack. Two symbols S and S are introduced to denote the regions on the plane 0 = z , occupied respectively by the crack and its exterior. It may be noted that and lower crack faces. Due to the symmetric condition occupied by the problem, the problem can be transformed to a mixed boundary-value problem of a half space 0 ≥ z subjected to the following boundary conditions on the plane 0 = z (Chen, 2004b; Wang, 2000): M A N U S C R I P T A C C E P T E D ARTICLE IN PRESS

  which lie on S . It is noted that 11 m depends on the elastic constants ij c only; 12 m is independent of the coefficients ( pertinent to the pressure field. This property can be verified by any standard symbol computation software. Substituting Eqn. (36) into Eqn. (35)

  Y. LI, et al / European Journal of Mechanics A/SolidsConsequently, the solutions thus obtained degenerate to those of the crack problems in thermoelasticity

  Y. LI, et al / European Journal of Mechanics A/Solids6. Verification of the solutionAs mentioned before, the pressure and temperature fields can be solved independently. Here, the potential theory method was presented to solve the pressure field governed by Eqn. (3) of a half-space subjected to an arbitrarily

  Y. LI, et al / European Journal of Mechanics A/Solids 17

(

  α and 3 α can be formulated in terms of ij c and the bulk modulus of the solid grain s

  Eqn. (23) as well as that in Case 1.

  Y. LI, et al / European Journal of Mechanics A/Solids 20 In view of Eqns. (B1.7) and (18), it is again seen that (

  Y. LI, etal / European Journal of Mechanics A/Solids 3

Fig. 1

 1 Fig. 1 Vertical (a) and horizontal (b) cross sections of a flat crack embedded in an infinite medium.

Fig. 2 A

 2 Fig. 2 A penny-shaped crack and the coordinate systems.

Fig

  Fig.3The dimensionless temperature

Fig

  Fig.4The dimensionless temperature

Table 1 :

 1 material parameters for thermoporoelastic analysis

Table 2 :

 2 The dimensionless stress E

	σ	z	and displacement	w	a	on I

Table 1 :

 1 Material parameters for thermoporoelastic analysis

	Parameter	Value	Unit
	E ( E′ )	1854 (927)	[ ] MPa
	K	s	12180

Table 2 :

 2 The dimensionless stress E

				σ	z	and displacement	w	a	on I
	ρ	σ	z	E		ρ		w	a
	a	Present		Chen et al.		a	Present	Chen et al.
	1.001	-0.025909095 -0.025909123 0.000	-0.012331027786	-0.012331041131
	1.010	-0.008174811 -0.008174820 0.100	-0.012269217734	-0.012269231012
	1.020	-0.005766139 -0.005766145 0.200	-0.012081890432	-0.012081903508
	1.030	-0.004696422 -0.004696427 0.300	-0.011763050799	-0.011763063529
	1.050	-0.003620044 -0.003620048 0.400	-0.011301573645	-0.011301585876
	1.100	-0.002529100 -0.002529103 0.500	-0.010678983317	-0.010678994875
	1.500	-0.001036622 -0.001036623 0.600	-0.009864822229	-0.009864832905
	2.000	-0.000669137 -0.000669138 0.700	-0.008806115238	-0.008806124769
	5.000	-0.000236576 -0.000236576 0.800	-0.007398616672	-0.007398624679
	10.00	-0.000116482 -0.000116482 0.900	-0.005374970399	-0.005374976216
	20.00	-0.000058022 -0.000058022 1.000	0	0
	50.00	-0.000023184 -0.

Acknowledgment

The work was supported by the National Natural Science Foundation of China (Nos. 10725210 and 10832009), and the National Basic Research Program of China (No. 2009CB623204). One author (Xiang-Yu LI) would like to express his sincere thanks to Mrs.

Yang Zhong-Xuan and Zhou Zhi-Qiang for valuable help.

Appendix A

The material coefficients in Eqns. (1) can be denoted by the engineering constants as follows


H Szrensen

A Danfaer

A B Strathe
email: strathe@life.ku.dk

H Sørensen

A new mathematical model for combining growth and energy intake in animals: The case of the growing pig

Keywords: Growth, Energy intake, Pigs, Stochastic differential equations

A simultaneous model for analysis of net energy intake and growth curves is presented, viewing the animal's responses as a two-dimensional outcome. The model is derived from four assumptions: 1) the intake is a quadratic function of metabolic weight; 2) the rate of body energy accretion represents the difference between intake and maintenance; 3) the relationship between body weight and body energy is allometric and 4) animal intrinsic variability affects the outcomes so the intake and growth trajectories are realizations of a stochastic process. Data on cumulated net energy intake and body weight measurements registered from weaning to maturity were available for 13 pigs. The model was fitted separately to 13 datasets. Furthermore, slaughter data obtained from 170 littermates was available for validation of the model. The parameters of the model were estimated by maximum likelihood within a stochastic state space model framework where a transform-bothsides approach was adopted to obtain constant variance. A suitable autocorrelation structure was generated by the stochastic process formulation. The pigs' capacity for intake and growth were quantified by eight parameters: body weight at maximum rate of intake (149 -281 kg); maximum rate of intake (25.7 -35.7 MJ/day); metabolic body size exponent (fixed: 0.75); the daily maintenance requirement per kg metabolic body size (0.232 -0.303 MJ/(day× kg 0.75 )); reciprocal scaled energy density (0.192 -0.641 kg/MJ 6 θ ); a dimensional exponent, ș 6 (0.730 -0.867); coefficient for animal intrinsic variability in intake (0.120 -0.248 MJ 0.5 ) and coefficient for animal intrinsic variability in growth (0.029 -0.065 kg 0.5 ). Model parameter values for maintenance requirements and body energy gains were in good agreement with those obtained from slaughter data. In conclusion, the model provides biologically relevant parameter values, which cannot be derived by traditional analysis of growth and energy intake data.

A c c e p t e d m a n u s c r i p t 1. Introduction

Growth in animals has been studied for centuries both from a biological and a mathematical view point and many growth functions are available [START_REF] France | Mathematical models in agriculture[END_REF].

These functions are especially valuable if growth and energy intake are studied for long periods of time where nonlinearities are visually detectable in data, i.e. sigmoid shapes and diminishing return behaviours. Analysis of intake and growth data is usual done separately by fitting nonlinear functions of time [START_REF] Kanis | Daily gain, food intake and food efficiency in pigs during the growth period[END_REF][START_REF] Lopez | A generalized Michaelis-Menton equation for the analysis of growth[END_REF]. Another approach for modelling growth and intake curves in animals is to describe the body weight (BW) as function of the cumulated intake [START_REF] Andersen | Growth and feed intake curves for group-housed gilts and castrated male pigs[END_REF]. This has an advantage because the efficiency of

A c c e p t e d m a n u s c r i p t

utilizing nutrients for growth follows directly as the derivative with respect to intake [START_REF] Andersen | Growth and feed intake curves for group-housed gilts and castrated male pigs[END_REF]. Growth and intake processes are interrelated because it can be assumed that the animal grows as a consequence of the nutrients and energy consumed (push effect) or it eats to meet the demands of its growth potential (pull effect) [START_REF] Emmans | A Method to Predict the Food Intake of Domestic Animals from Birth to Maturity as a Function of Time[END_REF]. Simultaneous modelling of growth and energy intake offers the possibility to formulate a mathematical model where the equations are linked so that more detailed information is extracted, e.g. maximum intake rate, maintenance requirements, body energy (BE) accretion. If complete records of growth and energy intake trajectories from individual animals are available, estimation of model parameters could be possible. In the following, we define dietary energy as net energy because net energy systems are currently the recommended energy evaluation system for formulating diets for growing mammals [START_REF] Just | The net energy value of balanced diets for growing pigs[END_REF]Noblet et al., 1994a;NRC, 2000).

Models in animal biology have usually been defined in terms of ordinary differential equations (ODEs) [START_REF] France | Mathematical models in agriculture[END_REF]. Another option is to use stochastic differential equations (SDEs). The main attraction of the SDEs and the primary difference to the corresponding ODEs is the inclusion of a diffusion term which accounts for model uncertainty. Model uncertainty is biological meaningful because there are many factors affecting intrinsic growth processes which cannot be explicitly modelled. Modelling growth and intake processes using SDE models provides the possibility of quantifying both the trend and variation in the processes. Furthermore, using SDE models instead of ODE models is a powerful approach to deal with serially correlated residuals that are likely to occur when analysing growth data [START_REF] Sandland | Stochastic growth curve analysis[END_REF][START_REF] Garcia | A stochastic differential equation model for height growth of forest stands[END_REF].

The objective of the present study is to develop a mathematical model for simultaneous analysis of serial net energy intake and BW measurements. From these easily obtainable data, the model should quantify the dynamics of growth and yield information on energy intake capacity, maintenance requirements and body energy gain.

Material and methods

A c c e p t e d m a n u s c r i p t 2.1 Animals and diets

The data used in this analysis is part of a serial slaughter experiment set up to determine growth capacity, energy and nutrient utilization in growing pigs (Danfaer and Fernandez, in preparation). A total of 183 pigs participated in the experiment. Data on intake and growth from 13 of these pigs were used in the present paper to estimate the parameters of the proposed model. The experiment was set up as a fully randomized block design with three blocks and six different treatment groups.

Each experimental group referred to a combination of two factors, i.e. gender (male, female and castrate) and a genetic sub-index for feed efficiency (high and medium level). The experimental pigs used in this study were crosses of YL sows (Landrace and Yorksire) and Duroc boars. The pigs were fed seven highly nutritious diets in the corresponding intervals: 4-7 weeks, 7 weeks-25 kg, 25-45 kg, 45-65 kg, 65-100 kg, 100-150 kg and 150 kg to maturity. The chemical composition of the seven diets is presented in Table 1. The pigs were housed individually under thermoneutral conditions and given ad libitum access to feed during the entire growth period to maximize growth.

All pigs were weighed at birth, weaning, weekly until approximately 150 kg BW and then every second week until the time of slaughter. The cumulated feed intake was record twice weekly from feed dispensers and is summarized to the corresponding dates of BW measurement. The total net energy intake is calculated as the net energy content of the diets (Table 1) multiplied by the cumulated feed intake in the corresponding periods of use. For the present analysis, we assume that the methods for prediction of net energy contents of the seven diets [START_REF] Just | The net energy value of balanced diets for growing pigs[END_REF]Boisen and Fernandez, 1997) are valid in the entire growth period from 7 -450 kg BW. This is somewhat simplistic, but is nevertheless supported by the findings of Noblet et al. (1994b). After slaughter, the chemical composition and energy content of the body was determined (Danfaer and Fernandez, in preparation).

Insert Table 1.

Mathematical derivation of the growth model

A c c e p t e d m a n u s c r i p t

In the following, we assume that animal growth and feed intake are driven by the state of the animal, which is described by the state variable BW [START_REF] France | Mathematical models in agriculture[END_REF]. Furthermore, it is assumed that allometry ( a X Y ∝

) can describe the differential growth of BW and BE. Note that the

θ θ × × - × = BW BW dt dNE (1)
where ș 1 , ș 2 , ș 3 are parameters. This functional relationship (1) has previously been used to describe protein accretion curves in boars of different breed [START_REF] Tauson | Pattern of protein retention in growing boars of different breeds, and estimation of maximum protein retention[END_REF]. The assumption that ingestion rate is a function of metabolic body size has also been used in dynamic energy budget models in animal ecology [START_REF] Nisbet | From molecules to ecosystems through dynamic energy budget models[END_REF]. Using the above function (1) implies that the rate of intake is symmetric in BW, but it does not imply that the rate of intake is symmetric in time because the growth rate is not constant over time. Using a quadratic function of metabolic size for modelling the rate of intake allows to parameterize the model in terms of the maximum rate of intake (ș 2 * ) and the BW at which the maximum rate (ș 1 * ) occurs. There are several advantages to reparameterize the model because the new parameters can be assigned biological meaning and reduce the intracorrelation between the two parameter estimates. The new parameters can be derived by differentiation of (1) with respect to

BW BW BW BW BW dBW dt dNE θ × × - × - = × - = BW ș BW ș ș BW dt dNE dt dBE (5)
where ș 4 is the maintenance energy requirement scaled to metabolic body size. The growth function derived here is without any assumption of an upper limit for the size of the animal. Since the mature size of an animal is an important biological growth parameter, the proposed growth model should accommodate calculation of such a trait. The basic idea built into the model is that the rate of intake approaches maintenance levels as the animal matures and thus the retention of BE should approach zero. Setting equation ( 5) equal to zero and solving with respect to BW yields the following equation, which can be used to calculate the mature BW of an animal based on the estimated parameters of the model

1 • © § 1 • © § = 3 2 4 1 max / log exp θ θ θ θ - BW (6)
A link function is needed in order to convert the state variable BE into BW and by assuming that the allometric relationship between BW and BE is valid:

6 6 1 5 5 /ș ș ș BW BE BE ș BW 1 • © § = ⇔ × = (7)
From this expression body energy density can be quantified throughout the entire growth period as M t = BE t /BW t , which is an indirect measure of the fatness of the animal. Different animals at the

× - × - × × × = × = BW BW BE dt dBE dBE dBW dt dBW (9)
However, BE is a latent state variable because the energy content of the animal was measured only once, i.e. after slaughter. Therefore further manipulations of ( 9) are needed for parameter estimation, which can be done by substituting (7) into (9). The manipulations produce two differential equations:

( ) ( )

0 2 2 4 1 6 1 5 5 0 2 2 1 ) 0 ( , ) 0 ( , 3 3 6 6 3 3 BW BW BW BW BW dt dBW NE NE BW BW dt dNE = × - × - × × 1 • © § × = = × - × = (10)
The parameter vector ș = (ș * 1 , ș * 2 , ș 3 , ș 4 , ș 5 , ș 6 ) T is to be estimated from data and the initial values y 0 = (NE 0 , BW 0 ) are obtained by using the measurements at t = 0 as starting values for the filtering and estimation procedures. We investigate if the parameter (ș 3 ) can be fixed at 0.75 because this is the traditional scaling factor applied in energy metabolism studies. This in turn makes it easier to compare the estimated maintenance requirement with published values. A complete list of all the parameters, their physical interpretation as well as dimensions is presented in Table 2.

Insert Table 2.

The stochastic growth model

Two sources of variation are considered in the model: 1) a dynamic noise term, which is part of the system (the animal) such that the value of the process at time t depends on this noise up to the time t, and 2) a measurement noise term, which does not affect the process itself, but only its observations. Separation of intra animal variation into two noise components requires additional A c c e p t e d m a n u s c r i p t explanation. The system noise can be thought of as animal intrinsic variability due to the following; 1) animal growth is always embedded in a randomly varying environment (no matter how well the experimental conditions are controlled), 2) growth processes are subjected to external and internal influences that change over time (e.g. shifting diets, sub-clinical diseases, ambient temperatures, hormonal influences, emotional stress etc.) which may randomly affect the growth. Modelling all of these aspects that disturb the system would produce a very large and complicated model that renders model identifiability and parameter estimation given data. In the foregoing section, we derived an ODE model based on biological assumptions regarding intake and growth. It is now translated into a stochastic state space model. Thus, the equations governing the growth and intake processes can be written as the following continuous-discrete time state-space model where the two system equations are:

( ) ( ) t BW NE dw dt BW BW BW BW BW dBW dNE 1 • © § + 1 • © § × - × - × × 1 • © § × × - × = 1 • © § × - θ θ θ 0 0 3 3 6 6 3 3 2 2 4 1 6 1 5 5 2 2 1 (11)
The first term on the right hand side of (11) is commonly called the drift term and the second is commonly called the diffusion term in which ı NE and ı BW are diffusion coefficients and w t is a two dimensional Wiener process with independent increments. The two measurement equations can be written as:

1 • © § + 1 • © § = 1 • © § k k t t k k e e BW NE y y k k , 2 , 1 , 2 , 1 (12)
The observation equations describe what is actually measured at discrete time points ( ) instead. For linear models it is possible to use the Box-Cox transformation approach, i.e. find the optimal (in a likelihood sense) transformation among the power and logarithmic transformations. A similar theory is not available for SDE models. Instead, we used a "transform both-sides" (TBS) approach, i.e. both responses and expected values in the model are transformed using the same transformation. The appealing feature about the TBS approach is that the parameter interpretations are unchanged. More specifically, the ODE system is transformed and then system noise is added to the transformed ODE system. Hence, system noise is assumed to have a constant intensity on the transformed scale. This approach can be generalized to accommodate any type of transformation.

N N k t t t t t < < < < = -1 1 0 ...
We may introduce z t = h(x t ) where x t denotes the model states (NE and BW). If the original ODE ( 10) is written as

dt t x f dx t t ) , , ( θ =
then the ODE for the transformed response

t t z x h = ) (
is obtained by the chain rule:

( ) ( ) ( )dt t z h f z h h dz t t t θ , ), ( 1 1 - - × ′ = ( 13
)
where h(x t ) is the transforming vector function, h´(x t ) is the derivative with respect to the states x t , h - 1 (z t ) is the inverse of h(•), and f(•) is the functional expression for the original ODE system (10). For

instance, if ) , ( ) ( BW NE x h t = then ) , ( ) ( 2 2 1 BW NE z h t = - and ) 2 / 1 , 2 / 1 ( ) ( BW NE x h t = ′ .
Using (13) and adding system noise yields the following SDE system ( )

( ) t BW NE t dw BW BW BW BW BW BW NE dz 1 • © § + 1 • © § × - × - × × 1 • © § × × × - × × = × θ θ θ 0 0 2 1 ) ( 2 1 3 3 6 6 3 3 4 2 2 4 1 6 ) 1 ( 2 5 5 4 2 2 1 (14)

A c c e p t e d m a n u s c r i p t

The transformed measurements are obtained as the square root of the original and then the additive measurement noise is added on the transformed scale. Visual inspection of the innovations plotted as function of time is used to choose between the two transformations (logarithmic or square root).

Parameter estimation

Given the model structure ( 11) and ( 12) approximate maximum likelihood (ML) estimates of the unknown parameters can be determined by finding the parameters ș that maximize the approximate likelihood function of a given sequence of measurements ( )

N N N y y y y Y , ,..., , 1 1 0 - = .
The likelihood function can be expressed as the following by conditioning on the observations at time t 0 i.e., ( )

0 , 2 0 , 1 0 ; y y y = , ( ) ( ) ( ) 1 • © § = = ∏ = - N k k k N N Y y p Y p Y L 1 1 , | θ θ θ (15)
where p denotes the probability density function conditional on the parameters and the previous observations. Since the SDE in ( 12) are driven by a Wiener process and since increments of a Wiener process are Gaussian, it is assumed that the conditional densities can be approximated by Gaussian densities, which means that a method based on the so-called extended Kalman filter (EKF) can be applied [START_REF] Kristensen | Continuous Time Stochastic Modelling -CTSM 2.3[END_REF]. The EKF filter is a recursive algorithm. This means that only the estimated state from the previous time step and the current measurement are used to compute the estimate for the current state because the previous state contains all the information up to that time point. The EKF has two distinct phases: predict and update. The predict phase uses the state estimate from the previous time step to produce an estimate of the state at the current time step. In the update phase, measurement information at the current time step is used to refine this prediction to arrive at a new, more accurate state estimate, again for the current time step (see appendix for full details regarding the EKF algorithm). Assuming a Gaussian density, which is completely characterized by its first and second order moments, the EKF is used because EKF describes the evolution of the first and second

A c c e p t e d m a n u s c r i p t

order moments of the conditional probability densities in terms of ODEs and algebraic equations. In order to formally introduce the assumption of Gaussianity, the following notation is defined

{ } θ , | ˆ1 1 | - -= k k k k Y y E y (16) { } θ , | 1 | 1 | - -= k k k k k Y y V R (17) 1 | ˆ- - = k k k k y y e (18)
and the likelihood function is rewritten in the following way:

( ) ( ) 1 • © § × 1 • © § - = ∏ = - - - N k k k k k k T k N R e R e Y L 1 1 | 1 1 | det 2 2 1 exp | π θ (19)
The parameter estimates can now be determined by solving the following nonlinear optimisation problem:

( ) ( ) { } N Y L | log min arg ˆθ θ θ - = Θ ∈ (20)
where, for each value of ș in the optimization, k e ˆ and

1 | - k k R
, are computed recursively by the EKF.

The estimation scheme has been implemented in the parameter estimation software Continues Time Stochastic Modelling (CTSM). More details of the parameter estimation, filtering and smoothing are giving in the CTSM User and Mathematics guides [START_REF] Kristensen | Continuous Time Stochastic Modelling -CTSM 2.3[END_REF].

The growth and energy intake trajectories obtained from the 13 pigs are analysed separately. Population parameter estimates and their standard deviations are then calculated as the average and standard deviation of the 13 individual estimates. This procedure is called a two stage approach [START_REF] Davidian | Nonlinear models for repeated measures data[END_REF].

Results and discussion

At a preliminary stage initial parameter estimates and their effect on convergence and final parameters estimates were monitored. In general, for all 13 animals in the data set the

A c c e p t e d m a n u s c r i p t

convergence properties of the estimation procedure were good because all estimations converged in less than 150 object function evaluations. The two transforming functions were applied and diagnostic plots were inspected for all 13 animals, which showed that the best variance stabilizing function was the square root. Therefore, the parameter values and likelihood tests presented in the following is based on the square root transformation. The parameters were all estimated by maximizing the likelihood, and the improvement of fit when ș 3 was allowed to vary freely compared to being fixed at 0.75 was measured by the likelihood value at convergence. This is presented in Table 3. In ten of the thirteen analysed data sets, the value of ș 3 did not deviate significantly from 0.75 (p > 0.05). As a consequence, ș 3 was fixed at 0.75 for all animals in the subsequent analyses. Parameter values for each animal obtained by fitting the model to data is presented in Table 4 together with population averages and coefficient of variation (CV) for the drift part of the model.

Insert Table 3.

Population based methods (mixed models) have been more widely applied in analysis of pig growth (Criag and [START_REF] Craig | Nonlinear mixed effects model for swine growth[END_REF][START_REF] Kebreab | Comparative evaluation of mathematical functions to describe growth and efficiency of phosphorus utilization in growing pigs[END_REF]. It is assumed that all animals follow the same functional form with parameters varying according to multivariate normal distribution. Thirteen animals as available in this data set are too few for application of a population based approach and thus the data from the 13 pigs were analysed separately. Population based methods are especially useful if the population values are of main interest. In the present study, the animal specific values are of primary interest (e.g. comparison with slaughter data), so the data was analysed pig by pig, but a population based analysis would be of interest, too. Such an approach would yield a complete description of the dynamics with animal-to-animal variation described by key structural parameters and the within-animal variation described by a stochastic process. An algorithm for non-linear mixed effects in combination with SDE-models has recently been developed [START_REF] Mortensen | A matlab framework for estimation of NLME models using stochastic differential equations -Applications for estimation of insulin secretion rates[END_REF][START_REF] Mortensen | PSM: Non-Linear Mixed-Effects modelling using Stochastic Differential Equations[END_REF], but needs further investigation.

Moreover, note that the current model has a two-dimensional response. None of the standard

A c c e p t e d m a n u s c r i p t

software packages, which are usually applied for nonlinear mixed modelling, will fit multivariate response models because these are developed for the univariate case [START_REF] Pinheiro | Mixed effects models in S and S-PLUS[END_REF][START_REF] Littell | SAS System for Mixed Models, Second Edition[END_REF].

Measures made on the individual animal are likely to be more correlated than measures made on different individuals, and serial autocorrelation is present since measurements made closely in time tend to be highly correlated compared to measures made further distant in time [START_REF] Sandland | Stochastic growth curve analysis[END_REF]). An autocorrelation structure is generated "automatically" due to the stochastic process formulation. Moreover, unequally spaced observations do not create a problem because the model is formulated in continuous time. This is an additional advantage compared to application of ARIMA models for modelling correlated residuals.

Insert Table 4.

Modelling growth and energy intake patterns

Growth curves and feed intake curves are often represented as a function of age because age can be recorded more easily and more precisely than for instance weight. The parameters ș 1 * (BW at maximum rate of intake) and ș 2 * (maximum rate of intake) are estimated in the range of 149 -281 kg BW and 25.7 -35.7 MJ NE/day. Our results indicate that the maximum intake was reached at a much later stage in life than usually reported in studies quantifying feed intake [START_REF] Andersen | Growth and feed intake curves for group-housed gilts and castrated male pigs[END_REF][START_REF] Lorenzo | Comparison of linear and nonlinear functions and covariance structures to estimate feed intake pattern in growing pigs[END_REF]. This is probably caused be the fact that measurements were conducted outside the normal BW range (20-120 kg). According to the knowledge of the authors there is no information available on patterns of feed intake of pigs beyond 150 kg BW, which excludes the possibility for comparison. However, our results agree very well with the theory proposed by [START_REF] Emmans | A Method to Predict the Food Intake of Domestic Animals from Birth to Maturity as a Function of Time[END_REF], which suggests that intake is maximised at about half of the mature size. A prerequisite for the validity of model is of course that the framework is able to describe the observed patterns of intake and growth. Those are presented in Figure 1, for the drift part of the model and the one-step-ahead predictions generated by the EKF together with observed data from two representative pigs in the dataset. In Figure 1 we have plotted NE, BW, intake, and

A c c e p t e d m a n u s c r i p t

growth rate against time to show that the model described data well. [START_REF] St | Use of genetic size-scaling in evaluation of animals growth[END_REF] suggested that intake of animals with a large appetite might decline after having reached a maximum. [START_REF] Fowler | Growth in mammals for meat production[END_REF] found that the daily intake reached a peak at about 120 kg BW and then declined to very fluctuating values. The feed intake data presented by [START_REF] Thompson | The effect of genotype and sex on the patterns of protein accretion in pigs[END_REF] suggested that the feed intake did not stabilize at 200 days of age. Moreover, large fluctuations in the feed intake during the course of time were observed as well, which corresponds well with the observations reported in this study. A feature of the model is that the rate of intake increases quickly in the early part of the growth period and then plateaus as the pig matures, which is a consequence of the energy being diverted towards maintenance. Moreover, it has the flexibility to identify that rate of intake has a peak and then decreases towards a plateau level, which corresponds to maintenance of body functions at maturity.

Insert Figure 1.

Modelling BW as a function of age requires a four parameter function like the Bridges, generalized Michalis-Menten or Richards functions for adequately describing the growth trajectory of pigs from birth to maturity (Kreab et al., 2007;Strathe et al., submitted). [START_REF] Lorenzo | Comparison of linear and nonlinear functions and covariance structures to estimate feed intake pattern in growing pigs[END_REF] presented an analysis of linear and nonlinear models and reported that a three parameter logistic function was appropriate for characterizing the pattern of feed intake in pigs. Thus, a total of seven structural parameters are needed when the responses are modelled separately. The proposed bivariate model uses six parameters in its complete form. However, the dimensionless exponent ș 3 was fixed 0.75 during parameter estimation and thus restricting the number of free parameters to only five plus optional two initial state parameters. In terms of the number of structural parameters that has to be estimated the bi-variate approach is not much different than modelling the responses separately. However, the proposed framework offers other quantities (maintenance requirement, energy gain) to be derived compared to the traditional analysis of growth and intake, which are usually confined to estimate initial value, maximum value and inflexion points. The modelled

A c c e p t e d m a n u s c r i p t

growth curves presented in Figure 1 have clear sigmoid shapes and their derivatives with respect to time increase in the early part of life towards a maximum and then decrease towards zero.

Figure 2 presents goodness of fit plots. The observed rates of energy intake and BW gain are plotted against the corresponding one-step-ahead predictions generated by the EKF for the 13 animals and the line Y = X is also included. The model fits well to the observed data because the model dynamically tracks the fluctuations in the observations. A comprehensive review of [START_REF] Whittemore | Technical review of the energy and protein requirements of growing pigs: Food intake[END_REF] concluded that prediction of short time fluctuations in feed intake was difficult, although the information is essential for the practical purposes of pig's day to day nutrition. Thus, application of SDE models may be a way to obtain precise predictions of short time fluctuations in intake and growth rates.

Insert Figure 2.

The biological implications of stochasticity

In Table 5 parameter estimates related to the stochastic part of the model, e.g. diffusion and measurement error terms are presented. First, it should be noted that the estimated system noise was significant at the 5% level (t-score around 10) for all modelled animals. The magnitude of difference between the system and measurement noise indicates that the measurement error is negligible compared to the system error for most of the pigs. Note however that it is very difficult to distinguish between the two sources of variation, in particular for equidistant observations. For the bodyweight, the total variance of the difference between two consecutive measurements 1 , 1 y and 2 , 1

y is approximately BW BW S t 2 2 + Δ σ .
If ǻt is the same for all pairs of consecutive measurements, then only the sum, not the two distinct terms is identifiable. The irregular sampling scheme with weekly and bi-weekly sampling improves identification somewhat, but simulations have indicated that there are still problems for this sampling scheme. Still we can conclude that the animal intrinsic variability is more pronounced than the measurement noise, and fitting a model without measurement errors gives very similar results for the fixed effects. This

A c c e p t e d m a n u s c r i p t

suggests that the incorporation of model uncertainty through a diffusion process is successful from both a biological and a statistical point of view.

Insert Table 5.

Previous attempts to model growth and intake processes are based on models, which are solutions to ODEs [START_REF] France | Mathematical models in agriculture[END_REF][START_REF] Lopez | A generalized Michaelis-Menton equation for the analysis of growth[END_REF]. SDE models provide a much more realistic representation of reality than their deterministic counterparts because system noise reflects that we do not have a full understanding of the complex processes, which regulate intake and growth. For this reason, the aim of the current model formulation has been focussed on modelling both the trend of the growth processes and the random fluctuations by means of SDEs. It further implies that an individual animal does not have a fixed growth curve from birth as growth is subjected to disturbances during life time. This can be demonstrated by a simulation study where the mean parameter values for the modelled growth and intake trajectories are used. The mean parameter vector is given as ș = (ș 1 * =221, ș 2 * =31.5, ș 3 =0.75, ș 4 =0.27, ș 5 =0.32, ș 6 =0.81, ı NE =0.17, ı BW =0.07) T . An Euler-Maruyama scheme was used with an integration step of 0.5 day [START_REF] Kloeden | Numerical solution to stochastic differential equations[END_REF]. The results of the simulation analysis are shown in Figure 3. The figure is organized as follows: plots A and B show the simulated paths of the stochastic growth model along with the mean curves for the stochastic process; plots C and D show the derivatives with respect to time in order to display the hypothesis build into the model. In contrast to previously published deterministic approaches (e.g. [START_REF] Lopez | A generalized Michaelis-Menton equation for the analysis of growth[END_REF] to model growth and intake, the trajectories are not smooth. In fact, a deterministic model assumes that: 1) the mathematical processes generating the observed intake and growth trajectories are smooth (continuous and continuously differentiable) within the considered time frame; 2) the variability of the actual measurements is due to observation error. The present approach results from the hypothesis that the underlying mathematical process is not smooth, but is subjected to random intrinsic perturbations.

This system noise represents the complexed effect of many factors, each with a small individual effect, which are not explicitly represented in the deterministic part of the model (the drift term of

A c c e p t e d m a n u s c r i p t

the model). Figure 3 also clearly illustrates the consequence of assuming that the system noise is additive on a square root scale (NE and BW) because once the simulated trajectories are transformed back to the original scale, it can be seen that the fluctuations in the intake process increase as the animal matures. This phenomenon has previously been described for intake [START_REF] Fowler | Growth in mammals for meat production[END_REF][START_REF] Thompson | The effect of genotype and sex on the patterns of protein accretion in pigs[END_REF].

Insert Figure 3.

Estimation of maintenance requirement

Although the functional relationship might look complicated, the model was derived from basic biological ideas and thus biologically relevant results can be extracted from the analysis.

A maintenance component is estimated and expressed on a net energy basis. The maintenance component ș 4 in the model was estimated to be in the range of 232 -303 kJ/(day× kg 0.75 ) with a mean of 268 kJ/(day× kg 0.75 ) and a CV of 8.6%. In net energy systems the maintenance component is usually estimated by regressing the daily energy retention or heat production on daily ME intake, then the intercept represents the maintenance requirement [START_REF] Just | The net energy value of balanced diets for growing pigs[END_REF]Noblet et al., 1994). 4). There are several explanations to this fact: 1) the maintenance requirement was represented by one parameter in the model and it is likely that this representation is to simplistic because the chemical composition of the body, organ size etc. changes during the course of growth [START_REF] Tess | The Effects of Body Composition on Fasting Heat Production in Pigs[END_REF]; 2) fasting heat production measurements are usually conducted in the BW range of 20-110 kg and thus little information is available beyond this point; 3) the net energy contents of the seven diets used in the present study were calculated from dietary characteristics using the method proposed by [START_REF] Just | The net energy value of balanced diets for growing pigs[END_REF] and Boisen and Fernandez (1997). We will direct focus towards estimating a maintenance Where q represents the maintenance requirement and ǻBE is the measured increment BE from weaning (t 0 ) to the time of slaughter (t 1 ). The integral in ( 21) can be approximated be the trapezium rule, i.e.

( ) ( )

¦ = - - + × × - ≈ N k k k k k BW BW q t t NEM 2 75 . 0 75 . 0 1 1 2 / ( 22
)
The relationship between NE and ǻBE + NEM is shown in Figure 4. When a q value of 326 kJ/(day× kg 0.75 ) [START_REF] Just | The net energy value of balanced diets for growing pigs[END_REF] was used in the left plot, NE was overestimated in relation to measured NE intake. In the right plot, a q value was estimated from data by minimizing the weighted least sum of squares so that estimated NE would be in accordance with the measured NE intake. This procedure converged at 261 kJ/(day× kg 0.75 ) and the results including the line Y=X is shown in the right plot in Figure 4. The derived maintenance requirement can be regarded as a population requirement because each point in the graphs corresponds to one animal. This can be compared to the mean maintenance value of 268 kJ/(day× kg 0.75 ) obtained with the model (Table 4), which indicates that the model predicted a realistic value. Rather than interpreting this as an absolute maintenance value (i.e. fasting heat production), the maintenance parameter should be considered as an apparent maintenance value because it depends upon the net energy system that has been used to calculate the dietary values.

Insert Figure 4.

Modelling relation between body energy retention and body weight

An allometric relationship was chosen to represent the relation between BE and BW due to its biological interpretation. The reciprocal scaled energy density ș 5 differed the most

A c c e p t e d m a n u s c r i p t

between animals (0.192 -0.641 kg/MJ 6 θ ). All parameter estimates of ș 6 (0.730 -0.864) in the model were less than unity, which means that the derivative with respect to BE is continuously decreasing. This corresponds to increasing energy density of the pig, which is an obvious consequence of the BE gain changing towards lipid accretion as the animal matures. The BE state variable is a latent variable in the model and is therefore eliminated from the differential equations due to the functional relationship between BE and BW. However, growth curves were modelled accurately, which means that the model should be able to predict energy retention at various stages of growth. We evaluated the proposed relationship ( 7) between BE and BW by comparing it with data from the 170 pigs in the serial slaughter experiment. The results are presented in Figure 5 on log-log scales where line ( ) log( ) log( ) log( 65

BE BW θ θ + =
) is added to the plot. The average of the 13 estimates of ș 5 and ș 6 is used to construct the population relation. In general, the line agrees well with the observed data, which indicates that the pattern of energy gain is reflected in the NE and BW measurements. However, close inspection of Figure 5 also reveals that BW is underestimated for low BE values and overestimated for high BE values. Furthermore, considering that the estimated maintenance level corresponds closely to the data obtained from serial slaughter experiments suggests that the rules for partitioning of net energy for maintenance and energy gain are adequate.

Insert Figure 5.

Other functional relationships, like the augmented allometric, have been proposed for relating different body components, and for some datasets they may yield a better fit than the allometric [START_REF] Schinckel | Characterization of growth parameters needed as inputs for pig growth models[END_REF]. The estimated standard errors for ș 5 were large relative to the estimates for many of the animals, and the correlation between the estimators of ș 5 and ș 6 was generally large, around 0.8. This indicates that the estimation procedure has difficulties in identifying and distinguishing these two parameters. An even more flexible function would enhance this problem. The statistical procedure presented here is suited for multi-response modelling [START_REF] Kristensen | Continuous Time Stochastic Modelling -CTSM 2.3[END_REF] and thus serial measurements of body composition, e.g. ultra sound or computed tomography scans can be directly integrated as well. Adopting multivariate response

A c c e p t e d m a n u s c r i p t

modelling techniques provides the possibility to formulate structural equations that are fitted simultaneously to data. Analysis of energy balance data has also been approached by multivariate techniques and more information can be extracted from such data [START_REF] Van Milgen | Energy partitioning in growing pigs: The use of a multivariate model as an alternative for the factorial analysis[END_REF].

Conclusion

The focus in the present investigation was to develop a mathematical model for energy intake, body energy accretion and body weight gain. The model is derived from simple biological assumptions and it can be used to study the interplay between net energy intake and body weight gain. The model is based on two coupled differential equations that describe the evolution in energy intake and body weight gain during the life time of animals. Animal intrinsic variability affects intake and growth and is build into the model through a stochastic differential equation approach. Experimental data obtained on growing pigs is used for parameter estimation and verification of the biological significance of the parameter values. We showed that eight parameters quantifies the time trend and variation in the two dimensional outcome, i.e. body weight at maximum rate of intake (149 -281 kg); maximum rate of intake (25.7 -35.7 MJ/day); metabolic body size exponent (fixed: 0.75); the daily maintenance requirement per kg metabolic body size (0.232 -0.303 MJ/(day × kg 0.75 )); reciprocal scaled energy density (0.192 -0.641 kg/MJ 6 θ ); a dimensional exponent, ș 6 (0.730 -0.867); coefficient for animal intrinsic variability in intake (0.120 -0.248 MJ 0.5 ) and coefficient for animal intrinsic variability in growth (0.029 -0.065 kg 0.5 ).

In conclusion, the new mathematical model provides an important alternative to traditional analysis of energy intake and growth curves by treating the animal's response as two dimensional.

1 | k k k k k t x h y - -= (A.3) S C CP R T k k k k + = - - 1 | 1 | (A.4)
which predict the mean and covariance of the output at time t k given all past information available at time t k-1 . Here S is the variance of the measurement error and C is the linearization of the measurement equation, i.e.

1 | ˆ- = ∂ ∂ = k k x x x h C (A.5)
The second step of the recursions involves the innovation equation: .6) which determines the one-step-ahead prediction residual at time t k . The third step of the recursions involves the Kalman gain equation:

1 | ˆ- - = k k k k y y e (A
1 1 | 1 | - - - = k k T k k k R C P K (A.7)

A c c e p t e d m a n u s c r i p t

This is basically a matrix of weights that determines the degree of updating in the fourth step because updating is based on a compromise between the observation and current model state. In a situation where the model is good, but the observations are dominated by measurement error, the state estimate should rely more on the model as opposed to fitting the observations. On other the hand, if the model is incomplete the states should rely more on the observations than the model. The updating equations are given by:

k k k k k k e K x x + = -1 | | ˆ (A.8) T k k k k k k k k k K R K P P 1 | 1 | 1 | | - - -- = (A.9)
where the predicted mean and covariance of the states at time t k are updated based on the observed value via the innovation and the Kalman gain. The fifth and final step of the EKF algorithm involves the state prediction equations:

( ) θ , , ˆ1 | 1 | t x f dt x d k t k t - -= (A.10) T T k t k t k t A P AP dt dP σσ + + = - - - 1 | 1 | 1 | (A.

ˆ-

= ∂ ∂ = k t x x x f A (A.12)
The new predicted mean and covariance of the states are then used to predict new values for the next observation and thus steps from A.3 to A.12 are repeated. [START_REF] Just | The net energy value of balanced diets for growing pigs[END_REF] and Boisen and Fernandez (1997). b Likelihood ratio for comparison of the two models M 0 and M 1 .

A c c e p t e d m a n u s c r i p t

c

The corresponding p-value in the chi-squared distribution with one degree of freedom.

  Model diagnostics (residual and autocorrelation plots) are employed and visual inspection of the innovations plotted as a function of time is used to pinpoint if variance stabilization is required. It turned out that on the original scale variance was increasing with expected value and it was decided to work with transformed responses (square root or logarithmic)

  [START_REF] Schiemann | Energetische Futterbevertung und Energienormen[END_REF] estimated the maintenance requirement as 280 kJ/(day× kg0.75 ) for barrows at approximately 95-185 kg BW, whereas[START_REF] Just | The net energy value of balanced diets for growing pigs[END_REF] estimated 326 kJ/(day× kg 0.75 ) for growing pigs weighing from 20 to 90 kg BW. The maintenance estimates produced by the present model are on average slightly underestimated compared to those values (Table

  slaughter data because these can be used for validation. According to the definition of NE we may write the following equation:

Fig. 1 .Fig. 2 .

 12 Fig.1. The relationship between total net energy intake, body weight, net energy intake, growth rate

Fig. 5 .

 5 Fig. 5. Validation of the proposed relationship between body weight and body energy

  Figure 1.

  Figure 2.

  Figure 3.

  Figure 5.

Table 1 .

 1 Chemical composition of experimental diets.

	7	150 kg-	884	7.80		147	26.1	30.6	529	159	6.33	2.25	2.83	4.52		44.3	7.3	6.0
	6	100-150 kg	888	8.11		162	26.8	31.4	492	157	7.88	2.68	3.13	5.75		45.5	7.7	5.7
	5	65-100 kg	886	8.26		179	44.4	32.4	458	178	9.51	3.02	3.26	6.74		50.3	7.9	6.1
	4	45-65 kg	883	8.41		204	54.5	35.5	406	171	10.8	3.34	3.46	7.31		53.0	8.2	5.9
	3	25-45 kg	889	8.96		213	86.3	36.9	358	176	12.4	3.80	3.58	8.48		56.2	8.7	6.0
	2	7 weeks-25 kg	896	9.26		240	91.3	31.8	337	172	14.0	4.69	3.49	9.48		59.4	9.4	6.5
	1	4-7 weeks	899	9.57		267	93.7	27.4	317	157	15.6	5.49	3.58	10.4		61.4	10.2	7.1
	Diet	Period of use	Dry matter, g/kg diet	Net energy, MJ/kg diet a	Nutrients, g/kg DM:	Crude protein	Crude fat	Sugars b	Starch	Dietary fibre c	Lysine	Methionine	Cystine	Threonine	Minerals, g/kg DM	Ash	Calcium	Phosporus

a Estimated according to

Table 2 .

 2 Definitions of model parameters.

	Parameters Unit	Parameter interpretation
	ș 1	*	kg	Body weight at maximum rate of intake
	ș 2	*	MJ/day	Maximum rate of intake
	ș 3			Metabolic body size exponent
	ș 4		MJ/(day × kg 3 θ ) The daily maintenance requirement per kg metabolic body size
	ș 5		kg/MJ 6 θ	Reciprocal scaled energy density
	ș 6			Dimensionless exponent
	ı NE	MJ	Diffusion coefficient related to animal intrinsic variability in intake
	ı BW	kg	Diffusion coefficient related to animal intrinsic variability in
				growth
	S NE	MJ 2	Measurement noise
	S BW	kg 2	Measurement noise

Table 3 .

 3 Fit statistics for models M 0 (varying ș 3 ) and M 1 (fixed ș 3 ).

	Pig	Logl 0	a	Logl 1	a	LR b	p-value c
	5365	16.35	16.39	0.080	0.777
	5374	28.77	29.16	0.780	0.377
	5711	58.80	62.30	7.000	0.008
	5713	53.79	54.28	0.980	0.322
	5720	38.73	38.97	0.480	0.488
	6053	45.80	51.80	12.00	0.001
	6056	25.89	32.14	12.50	0.000
	6205	10.69	11.97	2.560	0.110
	6206	-7.178	-7.155	1.250	0.264
	6211	20.60	20.83	0.460	0.498
	6710	-0.270	0.815	2.170	0.141
	6713	22.57	23.27	1.400	0.237
	6716	18.86	19.76	1.800	0.180

a Log-likelihood value at final convergence for the models.

A c c e p t e d m a n u s c r i p t

The approximate maximum likelihood (ML) estimation of parameters in SDEs involves the use of state filtering. Consider the general SDE model with the following system and measurement equation:

The extended kalman filter (EKF) is a set of equations that provides an efficient recursive approach to approximate the conditional distributions with Gaussian distributions of a sequence of measurements ( )

The conditional densities are thus completely characterized by the one-step prediction error e k and the associated one step prediction covariance matrix R k|k-1 which are then used to construct the likelihood function.

Given the parameters ș, initial states

, and initial state covariance

, the first step of the EKF involves the output predictions equation: