NUMERICAL ESTIMATE OF ERROR CONNECTED WITH THE APPLICATION OF IMPERFECT TRANSMISSION
CONDITIONS IN COMPOSITE STRUCTURES
Gennady Mishuris1, Andreas Öchsner2,3
1Department of Mathematics
Rzeszow University of Technology, Rzeszow, Poland
2Centre for Mechanical Technology and Automation
3Department of Mechanical Engineering
University of Aveiro, Aveiro, Portugal
ABSTRACT: The evaluation of imperfect transmission conditions for modeling thin elastic intermediate layers in composite materials has been discussed in the framework of a FEM analysis. The type of the conditions under consideration depends essentially on the mechanical properties of the interphase and has been analytically derived earlier by asymptotic analysis in the case of soft, stiff and comparable value interfaces. However, such analysis enables researchers to extract the real error connected with introducing imperfect interface models in numerical calculations especially in zones of high gradients. In this paper, the accuracy of the transmission conditions has been verified by FEM analysis for different types of interfaces, loadings and relations between mechanical parameters. It could be verified that the errors are in the analytically predicted range or even better and that the transmission conditions fail only in an extremely small region near the free boundary.
1. INTRUDUCTION
Composite materials are usually considered as nonhomogeneous solids with perfect bonding between different phases of the composites (e.g. (1,3)). On the other hand, such structures, in fact, contain thin intermediate layers matching materials of the phases together. Moreover, features of the layer may play an important role and influence the composite properties. When a structure consists of components of essentially different sizes and properties, FEM analysis of the structure becomes very difficult. This fact follows from the necessity to construct a complicated mesh structure which, in turn, may provide unstable numerical calculations (e.g. (5,9)). The aim of the investigation is to verify transmission conditions numerically at least in cases when it is possible to receive stable calculations in order to estimate in value the appearing error and compare it with that asymptotically predicted.
Figure 1: Schematic representation of the problem.
We consider geometrically symmetrical elastic structures consisting of two thick layers matched together by a thin interphase layer exhibiting different material properties in comparison with the bonded materials (Fig. 1). Various loading conditions of Dirichlet type on the top of the dissimilar structure are applied.
Transmission conditions along imperfect interfaces of zero thickness have been evaluated analytically by asymptotic analysis in (2,7,8) and are collected within Table 1. Under the additional assumption that the material properties of the interface do not change in direction perpendicular to the interface, functions aj = aj(x) in Table 1 are calculated according to formulas from Table 2, where Young's modulus, E = E(x), Poisson's ratio, n = n(x), and interphase thickness, h = h(x), can vary along the imperfect interface. Let us note that, in case of the soft weakly compressible interface (n ~ 0.5) in plane strain states, the second condition in Table 1 is not yet valid and has to be replaced by the other one: [uy] = 0.
Table 1: Possible sets of transmission conditions along the line y = 0 depending on the relative properties of the thin intermediate layer. 2D case.
interface / transmission conditionssoft / [ux]-a2sxy = 0 / [uy]-a1sy = 0 / [sxy] = 0 / [sy] = 0
comparable / [ux] = 0 / [uy] = 0 / [sxy] = 0 / [sy] = 0
stiff / [ux] = 0 / [uy] = 0 / [sxy] + ¶/¶x(a3¶ux/¶x) / [sy] = 0
Table 2: Parameters aj(x) for the plane strain and plane stress case.
case / a1 / a2 / a3plane strain / / /
plane stress / / / 2hE
We will show in the paper by FEM analysis of the structure that the numerical error connected with these transmission conditions is comparable with that predicted by theory or even better depending on the applied loading. Finally, such important values, for practical numerical calculations dealing with biomaterial structures with thin interfaces, as the range of the edge effect zone, the validity zone of the mentioned transmission conditions and the singularity dominated zone will be discussed.
2. FEM SIMULATION
A commercial finite element software (MSC.Marc) is used for simulating the mechanical behavior of the thin intermediate layer which has the dimension 2h=H/100=0.01 and L = 10. The two-dimensional FE-mesh is built up of four-node, isoparametric elements with bilinear interpolation functions. In order to investigate the edge effect on both sides of the rectangle (cf. Fig. 2), a strong mesh refinement in this region is performed. Furthermore, the mesh is generated in such a way so that it is possible to evaluate the displacements and stresses along the axes of symmetry: lines A and B (cf. Fig. 1) and along the transition zone of the materials: lines C and D. (The lines Ce and De belong to the bonded material and the lines Ci and Di to the interphase). We have checked the constructed FEM mesh by comparing results with theoretical ones for some simple loadings applied to the completely homogeneous rectangle (see, for example, (4)). It turned out that the constructed mesh structure is appropriate for the tensile and shear loading, but not for the bending one.
Figure 2: 2D FE-mesh: strong mesh refinement in the investigated area.
First numerical simulations have been made under the simple uniform tensile and the simple shear load case (corresponding components of the given displacements v are equal to 0 and 1) for the symmetrical sample when the elastic constants, Young's modulus E± = E* = 72700 MPa and Poisson's ratio n± = n* = 0.34, of the aluminum alloy AlCuMg1 (2017) are assigned to the bonded material (cf. Fig. 1). Differing elastic constants are assigned to the intermediate layer and the calculations are carried out for plane stress and plane strain cases. We will discuss the obtained numerical results in details for only one case as an example: plane stress under uniform external tensile loading with the following interface constants: E = 813 MPa, n = 0.3. In Fig. 3 distribution of displacements and stresses along the line A (cf. Fig. 1) are presented. Note, that the stress component sx(0,y) is discontinuous at the interface boundaries, as it should be expected, while all other components are continuous. Although the ratio E/E* can be estimated only as 0.1 while 2h/H = 0.01, one can see that the distribution of the displacements within the interface exhibits a linear character which coincides completely with theoretical results (7).
Figure 3: Displacement and stress distribution along line A (cf. Fig. 1).
We checked numerically the validity of the transmission conditions at point x = 0. Therefore, the value of Duy/sy is compared with the constant a2 = 2h/(2m+l) in the plane stress tensile case along line A. Although both values have order 10-7 for E = 8138 MPa and n = 0.4999, the relative error takes the same magnitude of 10-7 which is essentially better then one can expect from the asymptotic estimate O(e) where e ~ 10-2 in our case. Note also that the value of Duy(x)/sy(x) practically does not change along the entire interface, and only near the external boundary the edge effect becomes essential. To show this fact, distributions of the displacements and stresses along five lines B, Ci, Ce, Di, De (cf. Fig. 1) are presented for the same example in Figs. 4 and 5. Let us note that the first, second and third transmission condition (cf. Table 1) is practically satisfied identically for the entire interface as it follows from Figs. 4 and 5 due to tensile loading. In the case of the simple shear loading the same results have been obtained for the last three conditions. These facts are the simple consequences of the symmetry in geometrical and mechanical properties of the sample under consideration.
Figure 4: Displacement distribution along lines B, C and D (cf. Fig. 1).
Figure 5: Normal and shear stress distribution along line B (cf. Fig. 1).
One can think that the displacements and stresses should behave in opposite way at the free edge in comparison with that presented in Figs 4b and 5. At the first glance, the interphase stiffness seems to increase near the free edge, because the displacement decreases. In Fig. 6 a final shape of the free edge boundary after the deformation is presented for a sample with simple tensile loading. It is clear that such a behavior of the displacements near the free edge is reasonable because of the contraction. In order to compare two limiting cases of Poisson's ratio, the shape is drawn for the same tensile sample in Fig. 9b with another Poisson's ratio of n = 0.0001. Now, there is practically no contraction of the interphase material. However, the distance from the free edge, measured where the stress component sy is no longer a straight line, is practically the same for both cases with an accuracy of 1 % in the reference coordinate system.
Figure 6: Shape of the interphase near the free edge for different Poisson's ratios.
Similar numerical results are obtained for numerous combinations of elastic constants, loadings and plane states. For the shear loading in the plane strain case all results are practically the same as they are in the case of the shear loading in the plane stress. This is evident from the problem formulation with taking into account the sample symmetry. For this reason, we do not present here corresponding results.
To estimate an influence of the edge effect, the relative length of this zone is presented in Fig. 7 for different cases under consideration. Except the case of the weakly compressible interface (n = 0.4999) in plane strain case for the tensile loading, all curves in Fig. 6 exhibit a similar behavior d = d(E/E-). Namely, for small values of the ratio E/E- ~ 10-3, d takes a value near 0.05 and for smaller values it becomes comparable with accuracy of d due to the definition (1 % accuracy in displacement compared from the center of the sample). Moreover, magnitude of Poisson's ratio influences the value of d slightly, except the mentioned case of the weakly compressible interface in the plane strain case under tensile loading. However, d = 0.05 means that l = 50h, so that the depth of the edge zone is 25 times longer than the thickness of the interphase (or even more).
Figure 7: Relative length of edge effect zone.
Let us consider now for comparison a more complicated case: a dissimilar structure consisting of two different materials: steel E+ = 210000 MPa, n+ = 0.3, and aluminum E- = 72700 MPa, n-= 0.34, matched together by the thin interphase with Young's modulus E = 813 MPa as above, but with the less favorite value of Poisson's ratio n = 0.4999 (weakly compressible interface). Thus, the considered interface is the soft one. Instead of the uniform external loading, we apply now a complex tensile loading defined as follows: uy(x,h/2) = vy(x) = x2/25, vx(x) = 0 in the plane stress case. In Figs. 8-10 the same corresponding results as in Figs 3-5 are presented. One can see that displacements still change linear within the interphase which corresponds to the theory (7).
Figure 8: Displacement and stress distribution for asymmetric case.
However, distributions of displacement and stress components along the interface are no longer practically constants. Nevertheless, the vector of stresses is continuous through the interface as it follows from Fig. 10b and as it has to be due to the transmission conditions (Table 1). In this case it is essentially more difficult to find the edge effect zone. To do this, we have introduced additional perturbations in stresses of zero main vector at the right and left hand sides of the rectangle and observed the difference in the solutions. Obtained results have been similar to those reported in Fig. 7.
Figure 9: Displacement distribution for asymmetric case along lines B, C and D.
Figure 10: Normal and shear stress distribution for asymmetric sample and asymmetric load along horizontal lines B, C and D.
To define a region where the transmission conditions are still valid, the jump of the displacement components and the respective normalized stresses from the conditions (cf. Table 3) are shown in Fig. 11 and 12.
Figure 11: Verification of the first transmission conditions along the imperfect interface.
Figure 12: Verification of the second transmission conditions along the imperfect interface.
From the first glance the same excellent agreement between the compared values exists. However, it is difficult to directly extract the relative error at the symmetry point x = 0 from the numerical calculation at all. To clarify this fact, Table 3 is presented where the relative error for the second and the fourth transmission conditions is presented in two different points (in point x = 0 the error has been calculated by extrapolation). Fortunately, it is still possible to determine the zone of validity of the transmission conditions with the 1 % accuracy criterion starting from the point separated from zero (cf. Figs 11b, 12b). Corresponding values are presented in Table 4.
Table 3: Point by point error in the second and fourth transmission condition for the asymmetric case in plane stress, E = 813 MPa, n = 0.4999.
x / / / Rel. error /0.0 / 8.621439 10-6 / 9.226322 10-6 / -7.016 10-2 / -3.088 10-2
3.0 / 9.221909 10-6 / 9.226322 10-6 / -4.785 10-4 / -2.728 10-4
Table 4: Estimates of different zones near the free edge of the dissimilar rectangle in plane stress. E = 813 MPa, n = 0.4999.