Previous Up Next

3  Examples

Data files are enclosed with the TOCHNOG distribution. These files have several purposes. First, they show the analysis capabilities of TOCHNOG. Second, these input files enable an easy start in performing calculations yourself. Third, the data items target_item and target_value are used in the data part, allowing you to test TOCHNOG on your specific computer system. This manual only contains results with some of the input files (examp1.dat, examp2.dat, etc..).

3.1  Example 1: Backward facing step (fluid flow)

The following specific case of the GAMM workshop [6] is studied here:

total length of channel L=22.

length inflow section l=3.

total height H=1.

height inflow section h=0.5.

viscosity ν=0.01.

density ρ=1.

maximum flow component of parabolic profile in inflow section Umax=1.

the flow volume is (2/3)hUmax=1/3.

Reynolds Re= Umax(Hh)/ν = 50.

The result for the minimum x-velocity vx at distance x=0.8 behind the step is studied. This result checks the correct modeling of the recirculation zone after the step. At the distance x=0.8 we find the minimum velocity is vx = −0.041. This happens at y=0.1 above the lower edge. The results of the GAMM workshop were in between −0.066 and −0.018.



3.2  Example 2: Confined compression of fluid filled porous material

In this 1D example, a porous solid saturated with fluid is compressed by a force F on the solid phase. The length of the test specimen is L, plane strain conditions hold, the Young’s modulus is E, Poisson’s ratio is ν, permeability kp, width of test specimen W, thickness T, area A=WT. Some non-dimensional variables are introduced

X = x/L           T = 
Hkp
L2
 t           P = 
Ap
F
H=
E(1−ν)
(1+ν)(1−2ν)

At X=1 free drainage holds (P=0) for T>0. At X=0 the boundary is not permeable (∂ P/∂ X = 0) for T>0. Using the non-dimensional variables, the analytical solution reads (see [1])

P = 
n=0
 
2
M
 sin(MX) exp(−M2T)
M = 
π
2
 (2n+1)

As a typical result, P at T=0.22 is checked; at X=0.4 the exact solution for the non-dimensional pressure P is 0.6.

Now we choose: L=1, W=1, T=1, E=1.e1, ν=0, kp=0.1 and F=1. Thus, for x=0.4 at t=0.22 we should find p=0.6. The numerical analysis (5 elements, time step size 0.01) gives p=0.607. The pressure distribution at T=0.22 is given below



3.3  Example 3: Plasticity in plate with circular hole

This example is taken from [13] (page 245). A plate (20   mm wide, 36   mm long) is perforated at the center with a hole with radius 5   mm. The plate is stretched in length direction by increasing the overall strain є in time. The other two boundaries are free to deform. Because of symmetry, only a quarter needs to be modeled. A mesh with 48 linear quadrilateral elements is used



The Young’s modulus E=7000   kg/mm2, Poisson’s ratio ν = 0.2, the Von Mises yield stress σy = 24.3   kg/mm2 and the strain hardening h=0.032E=224   kg/mm2. Plane stress conditions hold. The plasticity influence is characterized by σmean at the upper edge for v=0.0625. The elasto-plastic analysis gives σmean=12.4   kg/mm2 for v=0.0625. A purely elastic calculation would have given σmean=16.13   kg/mm2 .

3.4  Example 4: Two-dimensional convection and diffusion

This example demonstrates the effect of the spatial stabilization algorithm in 2D. A convection and diffusion of heat equation is analyzed on a 1 by 1 square. The two-dimensional mesh consists of distorted linear quadrilaterals



The convection velocity β=1 and the conductivity k=0.01. The boundary conditions for temperature are chosen such that the exact solution for a boundary layer in y-direction holds:

 T(y) = T(y=0) + (T(y=1)−T(y=0)) 
1− exp(β
y
k
)
1−exp(
β
k
)

where we choose T(y=0)=1 and T(y=1)=0. This is a severe test for the spatial stabilization algorithm. Many algorithms exist which solve this example exactly when using a one-dimensional domain, say with y-axis only, but few exist which do not show wiggles for irregular 2D grids. The node_dof records are initialized with temperature 1 as a first estimate for the solution field. The we check the results at x=0.7 and y=0.6. The exact solution is 1. The numerical solution with the 4-noded elements is 0.95. Splitting the elements in triangles (see control_mesh_split) would have given the solution 1.001. Triangles seem to behave better than distorted quads (in this example anyway). Both solutions are quite good however.

3.5  Example 5: Inverse modeling of indenter test

In an experiment, a plate (size 1 by 1) is clamped at the bottom and loaded at the middle top point by a concentrated load (−0.000 1 in y-direction)



The y-displacement at the point of loading is measured (−0.000 25571) and used to estimate the Young’s modulus of the material. Four quadratic elements are used to model the indenter test. Initially, the Young’s modulus is estimated to be 0.2. The top force and the measured displacement are such that the Young’s modulus should become 1 in the inverse iteration process (with this element mesh). The table below shows the estimated Young’s modulus versus the inverse iteration number

iterationYoung’s modulus
00.200
10.276
20.371
30.487
40.619
50.760
60.887
70.970
80.998
91.000

3.6  Example 6: Impact of metal plate, local mesh refinement

A plate has an initial velocity and impacts on a rigid disk. The initial mesh is plotted below



With this initial mesh, several time steps are taken. Then a local mesh refinement is applied. This local mesh refinement uses the residue in the equations to decide where the mesh should be refined. This strategy is applied several times. The final mesh is plotted below:



3.7  Example 7: Global mesh refinement over polynomial domain

This example shows that mesh refinement follows the boundaries of a domain; more precisely, if the nodes of an element edge are all placed on a specific geometrical entity (geometry_line, geometry_circle, etc.) then new generated nodes along that edge will also be placed on that geometrical entity; see refine_globally_geometry for this. Here we mesh a domain bounded at the left and right side by vertical lines (x=0 and x=10), bounded at the bottom by the polynomial y = −10 − 10x + x2 and bounded at the top by the polynomial y = +10 + x2 − 0.1x3. Initially only one 4-noded quadrilateral is used



The four nodes of the element are placed on the intersection points of the polynomials and the lines. After 4 global refinements the mesh looks like



Each refinement is accompanied by some remeshing in order to obtain a more regular mesh.

3.8  Example 8: Residual remeshing near heat spot

The boundaries of a rectangle (size 1 by 1) have prescribed temperature 0.



The middle point is heated so that it gets temperature 1. The residual remeshing is used to capture the temperature distribution near the heated spot in the middle. After some time, the temperature in the middle is imposed, and the residual remeshing did locate nodes more closely to the middle. The number of nodes is not changed; the nodes are replaced in the direction of high equation residual.



3.9  Example 9: 2D Impact of metal plate, rebuilding of mesh after large deformations.

A metal plate impacts on a rigid sphere. We study here the large deformation of the plate. The initial mesh is show below



In each time step, a new mesh is built to prevent the elements from becoming too distorted. State variables (velocities, stresses, etc.) are mapped from the old mesh onto the new mesh. The mesh after these time steps is shown below:



Finally we show the equivalent plastic strain distribution:



3.10  Example 10: Shear of low-tension material (cracking)

The top edge of a 1 by 1 plate is sheared over 0.01 to the right, and the bottom edge of the plate is sheared over a distance 0.01 to the left. The material cannot sustain positive principal stress exceeding 0.002 — it cracks. This cracking is modeled with a plasticity model which does not allow for principal stresses exceeding 0.002; the group_materi_plasti_tension model is used for this. The plate is modeled with linear triangles (-tria3). The first plot below shows the (deformed) mesh.



The second plot below shows the largest positive principal stress; notice that this principal stress did become equal to 0.002 over the entire tension diagonal.



The third plot below shows the plastic strain; plastic strain does show up especially in the corners of the tension diagonal.



3.11  Example 11: Propagation of a disturbance in the wave equation.

The wave equation is used over a one-dimensional domain (between x=0 and x=1). The domain is divided into 128 linear elements. At time t=0, s=0 over the entire domain. After time t=0, at x=0 the scalar s is prescribed to hold the value 0 and at x=1 the rate of the scalar s (s) is prescribed to have the value 1.   10 −4. This disturbance at the right edge propagates into the domain with the speed of sound (c=1). At time 0.5 the rate of s over the entire domain is monitored; at this time point s should have become 1.   10 −4 in the right half of the domain, whereas nothing should have happened yet in the left half of the domain. The first plot shows s if we use purely explicit time stepping (control_timestep_iterations is set to 1).



It is clear that the disturbance did propagate into half of the domain, but quite some oscillations do show up. The oscillations are greatly reduced if two iterations are used (control_timestep_iterations is set to 2); see the second plot.



3.12  Example 12: Contact frictional heat generation.

A rectangular piece of material (size 1 by 1; Young’s modulus 1) is rubbed against a solid wall. In this example, we will analyze the stresses and temperature profile caused by frictional heat generation. The free edges of the material are prescribed to have temperature 0.



The deformed mesh (with 16 quadratic elements) at time 1 looks like:



Initially the right edge of the material did penetrate the solid wall over a distance 0.01, but the contact algorithm did eliminate this penetration. This resulted in a normal stress -sigxx of about size −0.01. The normal stress causes a frictional stress -sigxy of about size −0.0001 at the interface between material and solid wall, due to a friction coefficient of 0.01:



This frictional stress, in turn, causes frictional heat generation. The stationary temperature profile is plotted below:



3.13  Example 13: Automatically embedded tendon with local mesh refinement

A tendon is embedded over a diagonal of a circular solid.



This tendon has an initial stress. Due to this initial stress, the tendon will contract and deform the solid. The solid is modeled with linear triangular elements. TOCHNOG will automatically distribute the tendon over the solid elements (automatic embedment).

Time steps are taken, after which the mesh is automatically locally refined; this is done several times. The first plot shows the final deformed mesh. Notice that most refinement takes place where the tendon enters into the solid, but there is also some refinement of elements along the tendon.



The second plot shows the residue in the equations. Notice that the residue is high at the entrance of the tendon in the solid, but high residues also appear on both sides of the tendon.



3.14  Example 14: Continuous metal forming.

A metal sheet is extruded in a continuous process. The initial thickness of the metal is 1, and the thickness after extrusion is 0.9. The Young’s modulus is 7000 and the Poisson ratio is 0.2. The Von Mises plasticity model is used to model permanent deformations (yield stress 24.3). The first flow shows the finite elements mesh:



The extrusion is modeled by prescribing a velocity 1.0 at the right edge of the computational domain. On the left edge, the stresses and plastic strains are set to zero, in order to specify that virgin material enters the domain. Friction at the upper wall is not taking into account in this example. The result for κ, measuring the build up of plastic strains, is plotted below



3.15  Example 15: Generation of two holes in 2D mesh.

This example demonstrates some mesh manipulation. We start with one 4-noded quadrilateral element (of size 1 by 1). After four global refinements of the mesh, and splitting the 4-noded elements into 3-noded triangles the following mesh is obtained:



Finally a two circular holes are taken out of the mesh, and 10 remeshing passes are applied to get more nicely shaped elements



3.16  Example 16: Bearing capacity of foundation.

In this plane-strain example the Mohr-Coulomb plasticity law is used to calculate the bearing capacity of a foundation. The foundation is symmetric, so that only half of the problem is analyzed. The width of (half of) the foundation is 1.52 m. The size of the soil domain is taken to be 6 m by 6 m. The Young’s modulus of the soil is 207 ×   106 Nm−2 and the Poisson ratio is 0.3. For the Mohr-Coulomb law, both the yield rule angle and the flow rule angle are 20   degrees, and the cohesion is 0.069   106 Nm−2. At the bottom of the domain, all displacements are assumed to be fixed. At the left edge (the symmetry axis) and at the right edge, the horizontal displacement is fixed while the vertical displacement is free.

The classical solutions from Prandtl, Coulomb and Terzaghi give for the maximal average pressure σyy at the foundation values in the range 0.98   106 Nm−2 up to 1.20   106 Nm−2.

To get the same order of accuracy as in the classical solutions, a coarse mesh with only sixteen quadratic elements is used. For optimal numerical stability, the elements are integrated fully (integration points in the nodes). The foundation is prescribed a downward velocity of 0.05 m/s in the calculation. With time steps of 10−3 s the solution is advanced up to time 0.6 s. First the deformed mesh at time 1 s is plotted below.



Secondly, the development for the average pressure at footing over time is plotted



The maximum value of average pressure is near 1.1   106 Nm−2.

3.17  Example 17: Hertz contact problem.

A circular solid is compressed to a rigid surface. Due to symmetry only half of the circle needs to be analyzed. The compression is imposed, by prescribing the displacement of the middle line of the circle, so that in fact only one quarter of the circle needs to be analyzed. Plane strain conditions are assumed. The circle radius is 254   mm, the prescribed displacement of the middle line is 10.16   mm, the Young’s modulus is 206000   N mm−2 and the Poisson ratio is 0.3. The mesh contains 64 quadratic -quad9 elements.

The first plot shows the deformed mesh. Notice that the nodes at the bottom of the circle did not penetrate the line y=0.



The second plot shows the Mises stresses (the size of the deviatoric stress matrix) in N mm−2..



3.18  Example 18: Thermally induced stresses in plate.

A plate (of size 1 by 1) is initially stress free and in thermal equilibrium with its environment (all initial temperatures are 0). Then the left edge of the plate is prescribed a temperature 1 which will heat the plate. The plate starts convecting energy at its other edges; the coefficient for heat convection is 103, the convection environmental temperature is 0 and the conductivity of the plate is 50. The first plot shows the stationary temperature distribution. A mesh with 16 by 16 linear finite elements is used for the numerical analysis.



Due to the temperatures, thermal stresses will be induced in the plate. The plate cannot deform at its left edge, but is otherwise free to deform. The Young’s modulus is 206   109, the Poisson ratio is 0.3, and a plane stress situation is assumed (σzz is 0). The thermal expansion coefficient is 8.4   10−4. The second plot shows the size of the deviatoric stresses (the von Mises stress) which indicates if plasticity is to be expected.



3.19  Example 19: Nonlocal plasticity in softening bar.

A bar is clamped at the left edge and is subjected at the right edge to a prescribed deformation. A plasticity model is used to limit the possible tensile stress in the bar. This limit of the tensile stress is softened to zero with the effective plastic strain κ. Due to softening, unlimited localization would be possible. This is prevented, however, by applying a viscoplastic power law in combination with a nonlocal yield rule, see options_nonlocal. A nonlocal radius of 0.2 is used in the determination of the nonlocal yield rule.

At the right side of the bar we introduce a relatively weak spot, so that plasticity occurs first there. In the plots below the results for the velocity field vx, the effective plastic strain κ and the nonlocal yield rule fn respectively are shown. Please notice that the softening zone is of a limited, non mesh dependent, size.







3.20  Example 20: Three dimensional local mesh refinement near heat source.

A sphere has zero temperature at the outer surface, and contains a heat source in the middle such that the temperature in the middle is 1. The temperature distribution is calculated. A wire frame plot of the initial mesh (with -TET4 elements) is given below



Using the local refinement algorithm, 10 percent of the elements are refined after the time stepping shown in the figure. The residue in the temperature equation is used to select which elements need to be refined. Since the residue in the temperature equation is most large near the heat source in the middle, the local refinement algorithm places most elements near the middle. This can be seen in the wire frame plot below



3.21  Example 22: Constant volume cylindrical compression with Modified Cam Clay.

In soil mechanics, a constant volume cylindrical compression represents the so-called undrained triaxial test. The calculation was performed using a single -quad4 element under axisymmetric conditions. Only small deformations were considered. The parameters of the Modified Cam Clay model correspond to the Weald clay according to [2]: M=0.882, λ=0.088, κ=0.031, G=3000 kN/m2. Initial (effective) stress was isotropic with the mean stress p′=207 kN/m2. The void ratio e=0.627 remained constant during the compression due to the imposed constant volume condition. The initial preconsolidation mean pressure was p0=207 kN/m2.

The calculated stress-strain curve (q is the stress deviator and εa is the compressive axial deformation) in the left figure corresponds exactly to the numerically integrated equations of the model. The stress path in the right figure coincides with the analytical solution, see [12]:

  
pi
p
 = 


M22
M2i2
 


Λ



 

where η=q/p′ and Λ=(λ−κ)/λ. The subscript i denotes initial values.





3.22  Example 23: Automatic local mesh refinement near shear band.

A test specimen is subjected to a vertical compression. The initial mesh is given in the first plot



A hole is present in the down left corner; this hole triggers a shear band. Symmetry boundary conditions are applied, so that in fact the quarter part of a total test specimen is modeled. After the test specimen is compressed, automatic local mesh refinement is applied; the elements with the highest plastic strains are refined. Then the calculation restarts with the new mesh automatically. This refinement + automatic restart is done three times. The final plot shows the deformed mesh after the three refinements.



3.23  Example 24: 3D Impact of metal plate, rebuilding of mesh after large deformations.

This is essentially the same calculation as example 9, but now in 3D however. The metal plate is modeled with -tet4 elements. At the end time, the sphere starts completely penetrating the plate. We only show the equivalent plastic strain distribution in the final mesh.




Previous Up Next