TWO BOUNDARY VALUE PROBLEMS INVOLVING AN INHOMOGENEOUS VISCOELASTIC SOLID

. Recently, a thermodynamically consistent non-linear constitutive equation has been developed to describe the large deformation cyclic response of viscoelastic polyamides (see [17]). In this paper, two boundary value problems within the context of the above model, namely the stress relaxation of a right circular annular cylinder subject to twisting, and the inﬂation of a sphere are studied. In addition to solving the above problems numerically, investigation of the merits and pitfalls of studying the same boundary value problem for a special class of inhomogeneous body and its homogenized counterpart is undertaken. This study ﬁnds that for moderate strains the diﬀerences in relaxation time between the actual inhomogeneous body and its homogenized counterparts may be signiﬁcant.

1. Introduction. Two particles are said to be materially uniform if, for some placement of the two particles along with their neighborhoods, they exhibit indistinguishable mechanical response. If all the particles that constitute a given body are materially uniform with respect to a single placement then the body is said to be homogeneous (cf. [25]). A body that is not homogeneous is said to be inhomogeneous. Examples of inhomogeneous bodies abound. In fact, all bodies are inhomogeneous. Apart from naturally occurring bodies, modern manufacturing also often produces bodies with predetermined inhomogeneities, for example, composites.
For the ease in the analysis one assumes many inhomogeneous bodies to be homogeneous. It is expected that idealizing the body to be homogeneous would still result in an accurate estimate of the boundary traction required to realize a particular displacement or vice versa and the possibility of failure of the body for a given loading condition. This homogeneous approximation seems to be reasonable under some circumstances and not in some others. Hence, it is of interest to document the similarities and differences in the solution obtained for the inhomogeneous body and its equivalent homogeneous counterpart.
An analysis to check the reasonableness of the homogeneous approximation can proceed in a multitude of directions. In this study, we are interested in determining the appropriateness or otherwise of replacing the inhomogeneous composite by an equivalent homogeneous body with regard to how the body responds to external mechanical stimuli. Clearly, it is not expected that the stress and displacement time histories at various points in the homogeneous and inhomogeneous body will be the same. Hence, there is a need to define in what sense the response of the inhomogeneous body is to be compared to that of the corresponding homogeneous body. Appropriate use entails a clear understanding of the underlying assumptions.
In the micromechanics approach ( [9]), it is required that the volume averaged stress and strain over a given representative volume element (RVE) be the same for the inhomogeneous body and its homogeneous approximation. The micromechanics approach needs to be used with care. It is possible that the material properties of the homogeneous approximation depend on the shape of the representative volume element, implying an undesirable homogenization procedure. It is also possible that the maximum stress in the real body is larger than the volume averaged stress, thus making such a homogenization procedure unsuitable to study failure. Further, material properties obtained might depend on the boundary conditions applied on the RVE. These limitations have been understood within the context of linearized elasticity and linearized viscoelastic models. While for large deformation elasticity these issues are being studied [20,24,3], these issues within the context nonlinear viscoelastic models are largely unexplored. This study aims at addressing some of the questions that arise when homogenizing nonlinear viscoelastic bodies.
An alternative to micromechanics based homogenization is asymptotic homogenization. Asymptotic homogenization is used for studying heterogeneous bodies with the length of the inhomogeneity being L, where L denotes the length scale of the body and a non-dimensional constant. In asymptotic homogenization the solution to the relevant boundary value problem is sought as a function of . Then, the effective parameters are obtained from the limiting case where the parameter tends to zero. A detailed exposition on asymptotic homogenization procedure can be found in [18] and further advances can be found in [11,10]. The method has been applied for homogenization of elastic ( [12,2]) and viscoelastic bodies ( [26,5,1]). However, as pointed out in [4,3], while this is a reasonable methodology the existence and uniqueness of the homogeneous parameters is an issue.
There are a number of studies (see [13,7,6] and the references cited therein) that find bounds on the material parameters in the homogeneous approximation. For the bounds to be useful they should be tight. By solving a number of boundary value problems Saravanan and Rajagopal [22,19,20,21] have shown that such bounds cannot be tight at least for some popular constitutive models within the context of large elastic deformations. Whether such tight bounds exists for a particular nonlinear viscoelastic solid is investigated here.
The global response of a general inhomogeneous body depends on the direction of loading. Hence, even if the constituents of inhomogeneous body are isotropic, global response of the homogenized counterpart could be anisotropic in that it can have different responses along different directions.
Hill [8] and Chung et al [5] have reported anisotropic constitutive relations for inhomogeneous bodies made of isotropic constituents. Due to change in the material symmetry of the constitutive relation, the relationship between the principal direction of the stress and the left Cauchy-Green tensor is different. Quantification of the actual differences in the stress distributions for constitutive relations of varying material symmetry is required; this aspect is not investigated in this study.
Here, a systematic procedure to study twisting of an annular cylinder and inflation of a sphere for rate type constitutive relations, developed using the framework presented in [15], is outlined.
It is evident from the literature that the predictive capability of the homogeneous approximation of a given inhomogeneous body using a constitutive relation of the same form as that of the constituents of the inhomogeneous body is limited. In fact, Suquet [23] found that the inhomogeneous body made up of a material whose constitutive relation is given by the Maxwell viscoelastic model cannot be approximated as a homogeneous body made of a Maxwell material, especially when the relaxation times of the constituent materials are different. While Suquet [23] used a viscoelastic fluid like model, here a viscoelastic solid model that describes the response of polyamide 6 is used to study whether the constitutive relation for the homogeneous approximation can be of the same form as that of the constituents of the inhomogeneous body. A thermodynamically consistent nonlinear viscoelastic solid constitutive relation reported in [17] is used in this study. A parameter in the model developed in [17] which is related to the relaxation time is assumed to have a continuous spatial variation in case of the inhomogeneous body and to be a constant for its equivalent homogeneous approximation.
It is found that all the boundary loads required to be applied on various surfaces of the body to engender a given boundary displacement in the actual inhomogeneous body and its homogeneous approximation cannot be predicted equally well when the form of the constitutive relation for the homogeneous approximation is same as that of the constituents of the inhomogeneous body. Further, it is found that even if the mean of the spatial variation of the material parameter is the same, its value obtained from correlating the boundary force (or moment or traction) differs with the actual spatial variation of the material parameter and geometry of the body. Examining the stress distribution in the actual inhomogeneous body and its homogeneous counterpart having the same form for the constitutive relation as that of the constituents of the inhomogeneous body, it is found that qualitative features of the stress distribution such as the sense of the stresses or the gradient of the variation in the stresses are captured reasonably well for the magnitude of strains under consideration in the problems studied.
The organization of the paper is as follows. In section 2, the two boundary value problems are formulated, the governing equations are derived and solved for the thermodynamically consistent model developed in [17]. The procedure adopted to arrive at the constant material parameters in the homogeneous approximation is outlined in the section 3. The value of the material parameters obtained in the homogenized approximation for different cases are reported in section 4. The stress distributions in the actual inhomogeneous body and its homogenized counterpart are also presented in section 4. The implications of the simulations carried out and the results are detailed in section 5. The final section is devoted to concluding remarks.
2. Boundary value problem. Two boundary value problems are considered in this study. These boundary value problems are formulated in this section.
2.1. Twisting of an annular cylinder. First boundary value problem considered is that of twisting of an annular cylinder. A cylindrical polar coordinate system is used to study this problem. The annular cylinder in the reference configuration is assumed to occupy, where (R, Θ, Z) represent the coordinates of a typical point (X) in the reference configuration of the body, R i is the inner radius of the right circular annular cylinder, R o is the outer radius of the annular cylinder, H is the height of the annular cylinder. In the deformed configuration X moves to x.
This cylinder is assumed to be comprised of materials that are incompressible. While the rotation of the surface defined by Z = 0 is zero, the surface defined by Z = H rotates by an angle, ψH. In order to be able to realize this rotation at the surface defined by Z = H, a torsional moment, M is applied to this end. The cylinder is prevented from deforming axially by an axial load, L. By applying a suitable radial component of the normal stress at the inner surface of the cylinder (R = R i ), P i , radial displacement of this surface is prevented. The outer surface of the cylinder defined by R = R o is assumed to be traction free.
Appealing to a semi-inverse method, a motion of the following form is assumed, where (r, θ, z) represent the cylindrical polar coordinates of a material particle in the current configuration which occupied the point (R, Θ, Z) in the reference configuration. r(R, t) is an unknown function of R and time, t, represents the deformed radial location of a point, ψ(t) represents the angle of twist per unit length and is a prescribed function of time.
The cylindrical polar coordinate components of the deformation gradient, F κ R for the above motion (1) is obtained as, Since, the cylinder is incompressible, det(F κ R ) = 1. This incompressibility constraint yields a first order differential equation in r(R, t) which can be solved to obtain, where, r i (t) represents the deformed inner radius of the cylinder at time t. Since, there is no radial displacement at the inner surface, r i (t) = R i it follows that, r = R. Hence, the expression for deformation gradient (2) simplifies to, and the material time derivative of the deformation gradient,Ḟ κ R is computed as, whereψ = dψ dt . The cylindrical polar coordinate components of the spatial velocity gradient is found to be, The left Cauchy-Green deformation tensor is computed as, The torsional moment applied on the cylinder to effect a given time history of the angle of twist per unit length, ψ(t), is evaluated as, where the Cauchy shear stress T θz is to be computed from the constitutive relation for the assumed motion. The last equality holds since for the assumed motion (1), the incompressibility condition requires that r = R. The radial component of the normal stress that needs to be applied at the inner surface is computed as, by integrating the radial component of the equilibrium equation and using the fact that the outer surface of the cylinder is traction free and r = R from incompressibility constraint. Here T rr and T θθ are respectively the radial and circumferential component of the Cauchy stress to be determined from the constitutive relation for the assumed motion field. Similarly, the axial load that is applied to prevent the cylinder from expanding is given by, where T zz is the z component of the normal Cauchy stress to be determined from the constitutive relation for the assumed motion. The last equality is obtained using the result in [25] on observing that the components of the Cauchy stress is a function of only r and that r = R in order to satisfy the incompressibility condition.

2.2.
Inflation of a spherical shell. The second boundary value problem studied is that of the inflation of a spherical shell. A spherical coordinate system is used to study this problem. A sphere in the reference configuration is assumed to occupy, represent the coordinates of a typical point in the reference configuration of the body. R i represents the inner radius of the sphere. R o represents the outer radius of the sphere. This sphere is assumed to be made of an incompressible material. The sphere is subjected to a radial component of the normal stress at the inner surface defined by R = R i and the outer surface (R = R o ) is traction free. The motion is assumed to be, where (r, θ, φ) represent the spherical coordinates of a material particle in the current configuration which occupied the point (R, Θ, Φ) in the reference configuration r(R, t) represents the inflation of the sphere and is a yet to be determined function of radial coordinate in the reference configuration (R) and time (t). The spherical coordinate components of the deformation gradient F κ R for this motion field is, Since, the material is assumed to be incompressible, det(F κ R ) = 1. Using this condition, the resulting first order differential equation in r is solved along with the boundary condition that r(R i , t) = r i (t) to obtain, Substituting Eq. (13) in Eq. (12), The material time derivative of the deformation gradient,Ḟ κ R is evaluated as, The spherical coordinate components of the spatial velocity gradient is computed as, The left Cauchy deformation tensor is obtained as, Integrating the radial component of the equilibrium equation and using the boundary condition that the outer surface of sphere is traction free, the radial component of the normal stress acting on the inner surface, P sph is obtained as, where, T rr , T θθ , T φφ are the respective spherical components of the Cauchy stress to be determined from the constitutive relation for the assumed motion.

Constitutive relation.
To establish the thesis of this study, a nonlinear viscoelastic solid like model proposed in [17] to model polymeric materials is used. The constitutive equation for this model is developed using the framework of materials with evolving natural configurations, the evolution of the natural configuration being determined by the maximization of the rate of entropy production (see [15]): and the evolution of B κ p(t) (the left Cauchy-Green tensor associated with the evolving natural configuration 1 ) is specified as, (20) where the upper convected Oldroyd derivative is defined as In this study, the material parameters µ 2 , ν and β are assumed to be constants; specifically: µ 2 = 0.609 MPa, ν = 3.705 MPa(s (2β−1) ) and β = 0.568. These parameters are the same as that reported in [17] which are obtained so that the stress relaxation behavior of polyamide 6 could be captured. The material parameter µ 1 is assumed to vary spatially in the inhomogeneous body and to be a constant in its homogeneous approximation. The intent here is to demonstrate that effect of variation in one parameter in the model. This could have been the parameter that influences the stiffness or the relaxation time of the material. It was felt that varying stiffness would cause a variation similar to that reported for elastic response and hence the parameter that influences the relaxation time is chosen for the analysis. A more comprehensive study could undertake the study of influence of variation in each of the other parameters and their interactions. It can be seen that the relaxation time, T relax for this constitutive relation is 2.3.1. Spatial variation of µ 1 . For the purpose of our study here, the material parameter µ 1 is assumed to vary only along the radial direction. This dependance on the radial spatial coordinate R, is prescribed using, where R i is the inner radius of the cylinder or sphere, R o is the outer radius of the cylinder or sphere. Note that, the normalized parameterR varies between 0 and 1 over the domain of the body. µ mean 1 denotes the mean of the spatial variation, that is, The following forms of inhomogeneity are considered: 2.3.2. Linear variation.
where δ and k represents the amplitude and frequency of material property variation.
where δ and k represents the amplitude and frequency of material property variation.
In this study, unless otherwise specified µ mean 1 = 1.978 MPa, a value that corresponds to the stress relaxation experiments on polyamide 6. Small variations from this mean value (about 5%, to 15%) are considered in this study. However, even this small 5 to 15 percent variation in the parameter µ 1 results in the relaxation time of the material (see equation (22)) changing by 50 to 100 percent.

Governing equations and their solution.
In this section the governing equations that needs to be solved for each of the boundary value problems is documented. The solution procedure adopted to solve the governing differential equations is outlined.

2.4.1.
Twisting of an annular cylinder. Motivated by the non zero cylindrical polar components of the corresponding B κ R given in Eq. (7), the non zero cylindrical polar components of B κ p(t) for twisting of an annular cylinder is assumed to be, where A(r, t), B(r, t), C(r, t) and D(r, t) are yet to be determined functions of the radial spatial coordinate and time.
Using the Eq. (28) in Eq. (20) and using the incompressibility condition, the following equations are obtained where and use is made of the relation r = R obtained from the incompressibility constraint. Thus, for a given R, the first order differential equations in time, (30) through (32) is numerically solved using ODE15s, a built in solver in MATLAB . Then, the torsional shear stress is computed as, is used in Eq. (8) to obtain the torsional moment applied to realize the angle of twist per unit length, ψ(t). The integrand in Eq.(9) is used to find the radial component of the normal stress at the inner surface, Similarly, the integrand in Eq.(10) is used to find the magnitude of the axial load applied, Inflation of a spherical shell. Motivated by the non-zero spherical components of B κ R in Eq. (17), the non-zero spherical components of B κ p(t) for inflation of a sphere are assumed to be, where, E(R, t) is an yet to be determined function of the radial spatial coordinate and time. Note that the assumed B κ p(t) satisfies the incompressibility requirement. Using the Eq.(37) in Eq.(20) results in following differential equation, where Eq.(13) is used to convert the differential equation with r and t as variables to R and t as independent variables. For a given R, the first order differential equations in time, (38) is numerically solved using ODE15s, a built in solver in MATLAB .
The integrand in Eq. (18) is then evaluated as from which the radial component of the normal stress applied at the inner surface determined.
3. Finding material parameters in the homogenized approximation. As discussed in introduction, there are various measures by which one can seek a homogeneous approximation for a given inhomogeneous body. Each of these measures have their own merits and demerits. Here, the equivalence is sought such that the boundary force (or moment or traction) time history required to realize a given boundary displacement time history is same between the actual inhomogeneous body and its homogeneous approximation. Also, as discussed in the introduction, the form of the constitutive relation in the homogeneous approximation is assumed to be the same as that of its constituents, except that the material parameter which varied spatially is now assumed to be constant. The boundary displacement time history considered corresponds to that of a stress relaxation experiment. Here too as in the case of a real experiment, the boundary displacement time history is not given by a heaviside function. The boundary displacement is varied linearly from zero to the targeted boundary displacement. This is done because as discussed in [14], it is possible that smooth solutions such as those sought here might not exist for the nonlinear model when the boundary displacement is not a continuous function of time. Again, since the viscoelastic model is nonlinear, it is possible that a different boundary displacement time history could yield an entirely different homogeneous approximation for a given inhomogeneous body. This aspect of dependance of the material parameters in the homogeneous approximation on the boundary displacement time history is investigated only to a limited extent. Particularly, only the influence of rate of loading to reach the targeted boundary displacement is studied.
The procedure detailed below applies for both the twisting of an annular cylinder and the inflation of a sphere: 1. The boundary force and moment (or traction) required to engender the given boundary displacement corresponding to that of a stress relaxation experiment (see Figure 1) for the inhomogeneous body corresponding to a given spatial variation of µ 1 is computed (see Figure 2). 2. The inhomogeneous body is replaced by an homogenized body of identical geometry. The constitutive relation for the homogenized approximation is assumed to be the same as that given in Eq. (19) and Eq. (20) except that the spatially varying material parameter µ 1 is replaced with a constant µ mean−hom 1 . The value of other material parameters (µ 2 , ν, β) are assumed to be same as that in case of the inhomogeneous body. 3. The value of µ mean−hom 1 is obtained through optimization such that the overall error, tot is less than 2 percent and the maximum error max is less than 1 percent for the boundary force (or moment or traction) required to engender the given boundary displacement corresponding to that of a stress relaxation  Figure 2. Applied boundary force on the inhomogeneous body to realize the desired boundary displacement experiment. The overall error, tot is defined as, where N is the number of time instances when the boundary force (or moment or traction) on the inhomogeneous body and its homogenized counterpart are compared. F inhom denotes the boundary force (or moment or traction) in the inhomogeneous body and F hom the same boundary force (or moment or traction). The maximum error max is defined as, It is noted that other limiting values may also be chosen for these error measures depending on the requirement of closeness of agreement between response of homogenized body to that of the inhomogeneous body. But the values used here appear to give a good approximation with the value of the inferred µ mean−hom and, where T relax−inhom is computed using µ mean 1 in Eq. (22) and T relax−hom is computed using µ mean−hom 1 in Eq.(22), a study is conducted on the variation of µ mean−hom with respect to the two boundary value problems studied here. Specifically for the twisting of an annular cylinder, investigations are carried to quantify the variation of µ mean−hom (also the percentage variation in the relaxation time (Eq.(43))) with (a) the thickness of the cylinder (b) actual spatial variation of µ 1 (R) (c) correlation of the torsional moment versus axial load versus radial component of the normal stress at the inner surface of the cylinder (d) magnitude of the maximum twisting angle (e) rate of twisting. Similarly, for inflation of a spherical shell, quantification of the variation of µ mean−hom with respect to (a) the thickness of the sphere (b) actual spatial variation of µ 1 (R) is reported. With this information conclusions are drawn on the feasibility of homogenizing the viscoelastic solid like body comprised of a material whose response is captured by the constitutive relation studied here.

4.1.1.
Twisting of an annular cylinder. In all the simulations carried out R o = 1 and H = 1. Since, the outer radius of the cylinder is taken as the characteristic length for this problem, R o = 1. H = 1 allows no distinction between angle of twist per unit length and angle of twist. Figure 3 shows the variation of µ mean−hom with respect to the thickness of cylinder for different spatial variations -linear, sine and cosine variation -of the parameter µ 1 when the cylinder is twisted at a rate of 0.04 rad/s to a maximum twist angle of 0.01 rad, so that the torsional moment required to engender a given angle of twist is the same both in the actual inhomogeneous body and its homogeneous counterpart. It can be seen from the Figure 3 that the µ mean−hom tends to µ mean 1 as the thickness of the cylinder decreases, irrespective of the spatial variation of µ 1 . A 15 percent spatial variation in µ 1 causes 9 percent variation in µ mean−hom (100 percent variation in relaxation time, from Eq. (22)) just due to the change in thickness of the cylinder even though µ mean 1 for these various thicknesses is same (see  ). For a given thickness of the cylinder depending on whether µ 1 increases or decreases with radial location changes the value of µ mean−hom from 0.91 to 1.09 for the thickest annular cylinder studied (see Figure 3a). It can be inferred from the Figure 3b that as the frequency of the spatial variation increases µ mean−hom tends to µ mean 1 , irrespective of the amplitude of the variation. For the same frequency and amplitude of spatial variation of µ 1 , the value of µ mean−hom changes depending on whether the spatial variation is sine or cosine. Thus, it can be inferred from Figure 3 that the value of µ mean−hom depends on the spatial variation of µ 1 and the thickness of the cylinder. The changes in relaxation time with thickness are far greater than the changes in µ mean−hom (see Figure 4).
The variation of µ mean−hom with the thickness of cylinder, for the linear variation of the parameter µ 1 along the radial direction, when the cylinder is twisted at a rate of 0.04 rad/s to a maximum twist angle of 0.01 rad, so that the radial component of the normal stress applied at the inner surface to prevent the displacement of the inner surface for the given angle of twist is equal (for both inhomogeneous and corresponding homogenized body), is plotted in Figure 5a. When the cylinder is twisted at a rate of 0.04 rad/s to a maximum twist angle of 0.01 rad and the parameter µ 1 varies linearly along the radial direction, equating the axial load required to maintain the cylinder at its native length while twisting, for the inhomogeneous body and its homogeneous counterpart, requires µ mean−hom to vary with the thickness of the cylinder as shown in Figure 5b. Comparing Figures 3a, 5a and 5b it is apparent that different correlations for the same boundary value problem and spatial variation of µ 1 requires the value of µ mean−hom to be different. The maximum observed difference is about 5 percent in the value of µ mean−hom between different correlations when the spatial variation of µ 1 is about 15 percent. An interesting as the thickness of the cylinder reduces for cases where µ mean−hom is obtained by equating the moment or the radial component of the normal stress at the inner surface, in the case of µ mean−hom obtained from equating the axial load required to maintain the length of the cylinder while twisting, the deviation from µ mean 1 increases as the thickness of the cylinder decreases. Thus, comparing Figures 3a and 5 it can be seen that the value of µ mean−hom obtained for a given thickness of the cylinder and spatial variation of µ 1 , depends on what is correlated -moment, stress or force -between the actual inhomogeneous body and its homogenized counterpart to engender identical displacement fields. Figure 6a depicts the variation of percentage change in relaxation time with rate of twisting for different thicknesses of the cylinder and linear variation of µ 1 by equating the torsional moment required to engender a given maximum angle of twist (0.01 rad). Figure 6b portrays the variation of the percentage change in relaxation time with rate of twisting for different thicknesses of the cylinder and linear variation of µ 1 by equating the torsional moment required to engender a given angle of twist when the time of twisting is held constant (0.25 s). Thus, while Figure  6a portrays the variation of the percentage change in relaxation time with rate of loading (0.02 to 0.049 rad/s), Figure 6b shows the variation of the percentage change in relaxation time with rate of loading as well as magnitude of angle of twist (0.005 to 0.01225 rad). For the time histories of boundary displacement studied, there is no significant variation in the estimated value of the percentage change in relaxation time, obtained by equating the torsional moment between the inhomogeneous body and its homogenized counterpart with rate of twist. While the results are not sensitive to strain rate for a fixed variation in µ 1 and for a fixed geometry, it is clear that changing the cylinder thickness and the magnitude of variation in µ 1 can significantly alter the relaxation times.

4.1.2.
Inflation of a spherical shell. In all the simulations carried out as part of this study, R o is taken as the characteristic length for the problem and hence, R o = 1. Figure 7 shows the variation of µ mean−hom with respect to the thickness of sphere for different spatial variations -linear, sine and cosine variation -of the parameter µ 1 (R) when the sphere is inflated at a rate of 0.004R o /s to a maximum deformed inner radius of 1.001R i , so that the radial component of the normal stress applied at the inner surface of the sphere to realize a given deformed inner radius of the sphere is same between the actual inhomogeneous body and its homogenized counterpart. It can be seen from Figure 7 that µ mean−hom tends to µ mean 1 as the thickness of the sphere decreases, irrespective of the spatial variation of µ 1 . A 15 percent spatial variation in µ 1 causes 8 percent variation in µ mean−hom just due to change in thickness of the sphere even though µ mean 1 for these various thickness is same (see Figure 7b). For a given thickness of the sphere depending on whether µ 1 increases or decreases with radial location changes the value of µ mean−hom from 0.93 to 1.07 (see Figure 7a). It can be inferred from the Figure 7b that as the frequency of the spatial variation increases µ mean−hom tends to µ mean 1 , irrespective of the amplitude of the variation. For the same frequency and amplitude of spatial variation of µ 1 , the value of µ mean−hom changes depending on whether the spatial variation is sine or cosine. Thus, it can be inferred from Figure 7 that the value of µ mean−hom depends on the spatial variation of µ 1 and the thickness of the sphere, as in the case of cylinder.
In this case too, µ mean−hom did not vary with the rate of inflation in the range 0.002R o /s to 0.0049R o /s or with the maximum deformed inner radius in the range 1.0005R i to 1.0012 R i . However, for brevity these results are not presented here.      of the Cauchy stress predicted by the homogenized approximation (obtained by equating the moment) is in excellent agreement with that obtained from the actual inhomogeneous body when both are subjected to identical deformation field time histories. Figure 9 displays the radial variation of the spherical components of the Cauchy stress at the instant when the maximum deformed inner radius of the sphere is attained (t = 0.25s) for a thick annular sphere with R i = 0.5R o subjected to inflation for different radial variations of the parameter µ 1 in the constitutive relation. It can be seen from Figure 9 that the magnitude and sense of the spherical components of the Cauchy stress predicted by the homogenized approximation (obtained by equating the normal stress on the inner surface) is in good agreement with that obtained from the actual inhomogeneous body when both are subjected to identical deformation field time histories.
For all the radial variations of the parameter µ 1 studied, deformation field time histories considered and different time instances at which the Cauchy stress determined, the magnitude and the sense of the Cauchy stress predicted by the homogenized approximation is in agreement with that obtained from the actual inhomogeneous body. These results are not shown here.

5.
Discussion. The material parameter in the homogenized approximation is arrived at by either requiring the torque employed to engender a given twist or the axial force required to maintain a given axial length for a specified twisting angle or the radial component of the normal stress to be applied for a given twist to prevent the displacement of the inner radius, to be the same for both the inhomogeneous and "equivalent" homogenized body. It is found that the material parameter obtained by different correlations are not the same. This variation in the material parameter depending on whether the forces or moments or traction are equated for the same boundary value problem suggests the inappropriateness of the constitutive relation used for the homogenized approximation, that is the homogeneous model cannot be described by the same class of constitutive relations as that of the inhomogeneous body. Also, in this study, contrary to the usual practice of formulating the constitutive relation for the homogenized approximation from the response of a representative volume element subjected to a uniform state of stress or a homogeneous deformation, inhomogeneous deformation fields are used to find the material parameters in the constitutive relation for the homogenized approximation.
The response of a given geometry of the body to three different spatial variations of the relaxation time having the same mean is compared. It is found that even though the mean of the spatial variations is the same, the responses of the respective inhomogeneous bodies are different. Similarly, for a given spatial variation of the relaxation time if the geometry of the body is varied such that the mean relaxation time is the same, the material parameters in the constitutive relation for the homogenized approximation changes. This dependance on the material parameters in the constitutive relation of the homogenized approximation on the size and shape of the body is consistent with the observation in [22,19,20,21]. Thus, erroneous idealization of an inhomogeneous body as a homogeneous body results in the material parameters in the homogeneous model sought varying with the geometry of the body. Further studies are needed to check whether this variation of the homogeneous material parameters with geometry explain the size and shape effect reported in the literature for heterogeneous bodies like concrete when modeled as a homogeneous body.
The material parameters in the homogenized approximation are insensitive to the temporal variation of the boundary displacement for the cases studied here. This is because even though the constitutive relation is nonlinear, the numerical simulations are carried out for small deformations. In case of small deformations, the nonlinear viscoelastic model studied can be approximated by a linear viscoelastic model. Also, the temporal variations considered are limited to changes in the rate of twisting to achieve a constant angle of twist. It is known that [16] for linear viscoelastic models the relaxation of the stress with time does not depend on the rate of twisting. Hence, it is not surprising that the material parameters found for the homogeneous approximation is independent of the rate of deformation. More studies are required to understand the role of temporal variation on the inferred homogeneous constitutive relation.
The agreement of the magnitude and sense of the Cauchy stress components determined using the homogenized approximation with that of the actual inhomogeneous body is in contradiction with the results documented in [22,19]. This may be because the simulations in this study, are limited to small deformations or because of the form of the inhomogeneities considered in this study, or both. It is worthwhile to point out that Saravanan and Rajagopal [22,19] observed qualitative differences in the stress field for piecewise constant variation of the material parameters, a form of inhomogeneity not studied here. Hence, this qualitative agreement in the stress distribution has to be cautiously interpreted and more studies are warranted before arriving at general conclusions. 6. Summary. Two boundary value problems corresponding to that of pure twisting of an annular right circular cylinder and inflation of a sphere were formulated for an inhomogeneous body comprised of incompressible materials. The nonlinear viscoelastic constitutive relation proposed in [17] is used to describe the material that the inhomogeneous body is comprised of. One of the parameters in the model which is related to the relaxation time is assumed to vary spatially and simulations were performed for three special spatial variations. Using the chosen nonlinear viscoelastic constitutive model, geometry of the body and boundary value problems, the limitations in arriving at homogenized models in describing inhomogeneous bodies is investigated.