License

Claus Johansen, Sønderborg, Denmark, 2015-05-01 -> 2020-06-15


Torsion as a thermal analogy in FEA




Overview



0) Contents:

1) Preface

2) The analogy

3) A simple "text-book" example, given twist

4) Finding the (reaction) torque

5) Finding stress by given torque

6) Multiple symmetry

7) Hollow sections, ex. 1, no symmetry used

8) Hollow sections, ex. 2, some symmetry used




Top Up Down Bot.

1) Preface

Treating torsion as a thermal analogy can give tremendous reductions in problem size, in up to 3 steps:

The method is rarely described in FEA and if then often not in depth. E.g. twist pr. unit length used as input, but mostly a torque is known, not the twist - and in the cases where the twist is known, the corresponding torque will still be asked. How to take advantage of (multiple) symmetry? How to treat holes/hollow sections? These problems will be treated below.

Basis for § 3: "A simple "text-book" example, given twist" is an ANSYS® manual. All FEA are done in the ANSYS® FEA program - the "Classic/APDL" version.

The analytical values are all calculated by the formulas in: ROARK "Formulas for Stress and Strain".




Top Up Down Bot.

2) The analogy

I first discovered the analogy in an ANSYS® manual, but it can also be found in various places on the Internet.

The shear stress under torsion can be computed as heat conduction, by simply swapping the letter designations. The reason is that it's basically the same differential equation that describes the phenomena:

Heat conduction: Eq.1   (1)
Torsion: Eq.2   (2)

The analogy is hereafter:

Heat conduction Torsion
3a ~ Temperature     3b ~ Stress function
4a ~ (½) Heat generation pr. unit mass     4b ~ Twist pr. unit length
5a ~ Reciprocal heat conductivity     5b ~ Shear modulus
6a ~ Temperature gradient x-dir.     6b ~ Shear stress yz
7a ~ (negative) Temperature gradient y-dir.     7b ~ Shear stress xz

The analogy covers also phenomena as: Electrostatics, membranes under pressure (soap film), seepage (permeability) through a porous medium, ideal flow (no rotation), etc.



Top Up Down Bot.

3) A simple "text-book" example, given twist

The solution to this example is as I found it in an ANSYS® manual, but similar descriptions can be found various places on the Internet.

There are no difficulties in this example, only moderate hurdles can be to define the material property and find the twist pr. unit length when the total twist is known for the physical structure.

Below are the steps in a typical analysis: Geometry, material, element type, mesh, support, load and post-processing.


3.1) Geometry

A cross-section in the bar is a solid quarter circle, radius ~ length of the straight sides is 10 mm. Here the geometry of the cross-section is represented by a plot of the area:

fig.3.1.1
Fig. 3.1.1: Geometry.


3.2) Material

The only material property needed is the heat conductivity. According to the analogy, as described in § 2, the heat conductivity "K" is simply the reciprocal of the shear modulus "G".

K = 1/G

If the specific value of the shear modulus isn't known the modulus of elasticity (Young's modulus) together with Poisson's ratio can be used as a basis. This will allow for consistent results if the results are compared with 3-D models.

G = E/(2(1+ν))

In a case with steel and units in N and mm, the values can be:

E = 205.e+3 [N/mm2].

ν = 0.3 [-].

G = E/(2(1+ν)) = 205.e+3/(2(1+0.3)) = 78846. [N/mm2].

K = 1/G = 1/78846 = 0.1268292683E-04


3.3) Element type

The element types shall preferably be a thermal 2-D element type. In ANSYS® these element types are no 55 or 77. No 55 is a triangular element and no 77 is quadrilateral.


3.4) Mesh

In this example the area is meshed with the quadrilateral element, no 77 - ANSYS® does a good job in meshing with quad's in 2-D, so why not use this element where ever possible, as quad's in general are better/more efficient than triangles.

fig.3.4.1
Fig. 3.4.1: Mesh.

The shown mesh, with element type 77, has 3229 elements, 9910 nodes and thereby also 9910 degrees of freedom. Though a relatively fine mesh still a small model.


3.5) Support

The support consists simply of defining a temperature of zero along the outer boundary of the cross-section.

fig.3.5.1
Fig. 3.5.1: Support.

The number of active degrees of freedom will be reduced moderately (to 9460) but it's probably without importance.


3.6) Load

Assume a physical bar, with the cross-section described previously, and with a length of one meter and with a forced deflection of 10°.

Then the natural steps in finding the twist per unit length will be:

The forced defection in radians - convert the angular deflections in degrees to radians:
θr = θ°·π/180 = 10·π/180 = 0.1745329252 [-].

The forced defection pr. unit length "twist" - The cross-section is modelled in millimetres, so the length must be converted into mm and the twist will be pr. mm:
θ = θr/1000 = 0.1745329252/1000 = 0.1745329252E-03 [1/mm].

Heat generation pr. unit mass - according to the analogy, as described in § 2, the heat generation pr. unit mass "q" is simply two times the twist:
q = 2·θ = 2·0.1745329252E-03 = 0.3490658504E-03.

The heat generation is evenly distributed (constant) over cross-section:

fig.3.6.1
Fig. 3.6.1: Load.


3.7) Post-processing

A plot of the temperature will give the "shape" of the potential:

fig.3.7.1
Fig. 3.7.1: The potential "temperature".

A plot of the thermal gradient will give the sear stress:

fig.3.7.2
Fig. 3.7.2: The shear stress "thermal gradient".

The max shear stress can finally be converted to an equivalent stress by the traditionally used relationship, in this case:
σ = τ·√3 = 82.26887512·√3 = 142.4938716 ~ 142. [N/mm2].



Top Up Down Bot.

4) Finding the (reaction) torque

During my long service as an engineer I never met a case, where stresses as a function of angular displacement "twist" only were asked. The problem was at least always followed by the question of the corresponding torque?

In a case as the text-book example in the previous section, where the geometry isn't naturally modelled around the geometric centre of the cross-section "CoG", there must be created a local coordinate system in this position.

In the Pre-processor, find the coordinates of the mass centre "CoG" by the commands:

 asum                                 ! sum area properties.
 *get,name,area,,cent,x $ xc= name    ! x-pos. of mass centre.
 *get,name,area,,cent,y $ yc= name    ! y-pos. of mass centre. 

Afterwards create a local coordinate system in this position - it's a large advantage if this local coordinate system is cylindrical, as shown here:

 csys,0 $ dsys,0
 local,11,1,xc,yc,0                   ! local coord.sys. at CoG.
 csys,0 $ dsys,0 

A plot with the local coordinate system indicated by the number "11":

fig.4.1
Fig. 4.1: The local coordinate system in CoG.

The torque is calculated with the facilities in the etable procedure, the principles sketched below:

fig.4.2
Fig. 4.2: Finding the torque.

The procedure that takes place in the commands below can be described as:

The commands are:

 csys,11 $ dsys,11 $ rsys,11          ! cyl. coord. sys. at centre of mass.
 etable,rcen,cent,x                   ! element centre radius.
 etable,tauy,tg,x                     ! therm. gradient, radial component.
 etable,area,volu                     ! volume (depth == 1) eq. area.
 smult,df,area,tauy                   ! tangential force pr. element.
 smult,dt,df,rcen                     ! torque pr. element.
 ssum                                 ! all etable values are summed up.
 *get,name,ssum,dt                    ! torque from sum to parametre "name".
 torq= -name                          ! torque in the shaft. 

The radial (X-dir.) thermal gradient corresponding to the tangential (Y-dir.) shear stress, that is used in the operations above, look like this:

fig.4.3
Fig. 4.3: Radial thermal gradient ~ tangential shear stress.

In the text-book example the torque that corresponds to the angular deflection of 10° for the 1000 [mm] long steel bar of a ¼ circular sector of radius 10 [mm] is calculated to [Nmm]:

 PARAMETER TORQ =     11329.10468 

The stiffness of the bar, or more precisely, the stiffness of a unit length of the bar, can now be calculated. The general formula is:

θr = (TL)/(KG)

Where:
T ~ Torque.
L ~ Length.
K ~ Stiffness.
G ~ Shear modulus.

By re-arranging the general formula K can be found to:

K = (TL)/(θrG) = T/(θG)

In this text-book example the stiffness "the effective torsional moment of inertia" of the bar of a ¼ circular sector of radius 10 [mm] is calculated to [mm4]:

 PARAMETER KEFF =     823.2613159 

This section "Finding the (reaction) torque" is essential for practical use of the analogy, and not only being of academic interest. It will also allow for a comparison with analytical formulas - stiffness and stresses at a given torque. I will below give an example of how the problem can be solved when the torque is known, and the stresses and deflection are the desired results. This represents to me a typical procedure in a torsional problem.



Top Up Down Bot.

5) Finding stress by given torque

The geometry is as in the first example, § 3, and the procedures of finding the corresponding torque are as described in § 4, but as here the torque, not the twist, is known, it calls for a 2 step analysis. I'll describe it from the bottom again:


5.1) Geometry

A cross-section in the bar is a solid quarter circle, radius ~ length of the straight sides is 10 mm. Here the geometry of the cross-section is represented by a plot of the area:

fig.5.1.1
Fig. 5.1.1: Geometry.


5.2) Material

The only material property needed is the heat conductivity. According to the analogy, as described in § 2, the heat conductivity "K" is simply the reciprocal of the shear modulus "G".

K = 1/G

If the specific value of the shear modulus isn't known the modulus of elasticity (Young's modulus) together with Poisson's ratio can be used as a basis. This will allow for consistent results if the results are compared with 3-D models.

G = E/(2(1+ν))

In a case with steel and units in N and mm, the values can be:

E = 205.e+3 [N/mm2].

ν = 0.3 [-].

G = E/(2(1+ν)) = (205.e+3)/(2(1+0.3)) = 78846. [N/mm2].

K = 1/G = 1/78846 = 0.1268292683E-04

Now if K is defined to this value in ANSYS®, the output shows:

 PARAMETER K =    0.1268292683E-04 

Then the command for the heat conductivity is:

 mp,kxx,1,k 

And the output shows:

 MATERIAL         1     KXX  =  0.1268293E-04 

5.3) Mesh

In this example the area is meshed with the quadrilateral element, no 77 - ANSYS® does a good job in meshing with quad's in 2-D, so why not use this element where ever possible, as quad's in general are better/more efficient than triangles.

fig.5.3.1
Fig. 5.3.1: Mesh.

The shown mesh, with element type 77, has 3229 elements, 9910 nodes and thereby also 9910 degrees of freedom. Though a relatively fine mesh still a small model.


5.4) Preparing load and torque calculation

Before leaving the pre-processor we can as well prepare load and torque calculation by issuing the "asum" command and extract polar moment of inertia (about CoG) and the position of CoG and place the coordinate system for (reaction) torque calculation here:

 asum                                 ! sum area properties.
 *get,name,area,,imc,z  $ ip= name    ! polar moment of inertia.
 *get,name,area,,cent,x $ xc= name    ! x-pos. of mass centre.
 *get,name,area,,cent,y $ yc= name    ! x-pos. of mass centre. 

And the output shows:

 PRINT GEOMETRY ITEMS ASSOCIATED WITH THE CURRENTLY SELECTED AREAS

 *** NOTE ***                            CP =       0.687   TIME= 22:12:13
 Density not associated with all selected areas.  Geometry items are     
 based on a unit density.                                                

  TOTAL NUMBER OF AREAS SELECTED =      1  (OUT OF      1 DEFINED)

  TOTAL SURFACE AREA OF ALL SELECTED AREAS =   78.540    

  TOTAL VOLUME OF ALL SELECTED AREAS =   78.540    

  CENTER OF MASS:  XC=  4.2441      YC=  4.2441      ZC=  0.0000    

                    *** MOMENTS OF INERTIA ***
          (BASED ON A UNIT DENSITY AND A UNIT THICKNESS)

        ABOUT ORIGIN     ABOUT CENTER OF MASS    PRINCIPAL
    IXX =    1963.5            548.78            713.49    
    IYY =    1963.5            548.78            384.07    
    IZZ =    3927.0            1097.6            1097.6    
    IXY =   -1250.0            164.71    
    IYZ =    0.0000            0.0000    
    IZX =    0.0000            0.0000    

    PRINCIPAL ORIENTATION VECTORS (X,Y,Z):
      0.707  0.707  0.000   -0.707  0.707  0.000    0.000  0.000  1.000
      (THXY=   45.000  THYZ=    0.000  THZX=    0.000)

 *GET  NAME      FROM  AREA  ITEM=IMC  Z     VALUE=  1097.56639    

 PARAMETER IP =     1097.566392    

 *GET  NAME      FROM  AREA  ITEM=CENT X     VALUE=  4.24413177    

 PARAMETER XC =     4.244131767    

 *GET  NAME      FROM  AREA  ITEM=CENT Y     VALUE=  4.24413177    

 PARAMETER YC =     4.244131767    

Afterwards create a local coordinate system in this position - it's a large advantage if this local coordinate system is cylindrical, as shown here, the commands are:

 csys,0 $ dsys,0
 local,11,1,xc,yc,0                   ! local coord.sys. at CoG.
 csys,0 $ dsys,0 

And the output shows:

 ACTIVE COORDINATE SYSTEM SET TO      0  (CARTESIAN)  

 DISPLAY COORDINATE SYSTEM SET TO      0  (CARTESIAN)  

 COORDINATE SYSTEM     11 DEFINITION.  TYPE= 1  (CYLINDRICAL)
   XC,YC,ZC=  4.2441      4.2441      0.0000    
   ANGLES=    0.00    0.00    0.00   PARAMETERS=   1.000   1.000
   ORIENTATION=  1.00  0.00  0.00   0.00  1.00  0.00   0.00  0.00  1.00

 ACTIVE COORDINATE SYSTEM SET TO     11  (CYLINDRICAL)

 ACTIVE COORDINATE SYSTEM SET TO      0  (CARTESIAN)  

 DISPLAY COORDINATE SYSTEM SET TO      0  (CARTESIAN)   

A plot with the local coordinate system indicated by the number "11":

fig.5.4.1
Fig. 5.4.1: The local coordinate system in CoG.


5.5) Support

The support consists simply of defining a temperature of zero along the outer boundary of the cross-section.

fig.5.5.1
Fig. 5.5.1: Support.

The number of active degrees of freedom will be reduced moderately (to 9460) but it's probably without importance.


5.6) Initial load

Let the torque be 1.e+4 = 10 000 [Nmm].

As the formula for angular deflection of a circular bar is:

θr = (TL)/(IpG)

Where:
T ~ Torque.
L ~ Length.
Ip ~ Polar moment of inertia.
G ~ Shear modulus.

Then the expression for angular deflection of a circular bar pr. unit length "twist" is:

θ = T/(IpG)

With:
T = 1.e+4 = 10000. [Nmm]. Given torque.
Ip = 1098. [mm4]. Polar moment of inertia as described in § 5.4.
G = 78846. [N/mm2]. Shear modulus of steel as described in § 5.2.

The first guess on twist in the bar is the same as if the bar was circular - this is not the case, but it gives a first guess with some physical sense and a value that can have a reasonable size. Alternatively, a unit value could be used. In the end it should make no difference as the result of the initial load set will be scaled to obtain the correct final result.

The values below are only first guesses and will be designated by index 1.

The angular deflection of a circular steel bar pr. is unit length "twist" is:

θ1 = 10000/(1098·78846) = 0.1155549853E-03 [1/mm].

Heat generation pr. unit mass - according to the analogy, as described in § 2, the Heat generation pr. unit mass "q" is simply two times the twist:

q1 = 2·θ1 = 2·0.1155549853E-03 = 0.2311099705E-03

The heat generation is evenly distributed over cross-section:

fig.5.6.1
Fig. 5.6.1: Initial load.


5.7) Finding the initial (reaction) torque

The procedure is described in depth in § 4.

Solve the problem and go to the post-processor. The commands are:

 csys,11 $ dsys,11 $ rsys,11          ! cyl. coord. sys. at centre of mass.
 etable,rcen,cent,x                   ! element centre radius.
 etable,tauy,tg,x                     ! therm. gradient, radial component.
 etable,area,volu                     ! volume (depth == 1) eq. area.
 smult,df,area,tauy                   ! tangential force pr. element.
 smult,dt,df,rcen                     ! torque pr. element.
 ssum                                 ! all etable values are summed up.
 *get,name,ssum,dt                    ! torque from sum to parametre "name".
 torq1= -name                         ! torque in the shaft. 

The radial (X-dir.) thermal gradient corresponding to the tangential (Y-dir.) shear stress, that is used in the operations above, look like this:

fig.5.7.1
Fig. 5.7.1: Radial thermal gradient ~ tangential shear stress.

And the output shows:

 ACTIVE COORDINATE SYSTEM SET TO     11  (CYLINDRICAL)

 DISPLAY COORDINATE SYSTEM SET TO     11  (CYLINDRICAL)

 RSYS KEY SET TO  11

 USE COORDINATE SYSTEM  11 FOR SOLUTION RESULTS

 STORE RCEN     FROM ITEM=CENT COMP=X     FOR ALL SELECTED ELEMENTS

 STORE TAUY     FROM ITEM=TG   COMP=X     FOR ALL SELECTED ELEMENTS

 STORE AREA     FROM ITEM=VOLU  FOR ALL SELECTED ELEMENTS

 OPERATION SMUL  RESULT= DF        OPERAND1= AREA      OPERAND2= TAUY    
  FACTOR1=  1.0000       FACTOR2=  1.0000       CONSTANT=  0.0000    

 OPERATION SMUL  RESULT= DT        OPERAND1= DF        OPERAND2= RCEN    
  FACTOR1=  1.0000       FACTOR2=  1.0000       CONSTANT=  0.0000    

 SUM ALL THE ACTIVE ENTRIES IN THE ELEMENT TABLE 

 TABLE LABEL     TOTAL
 RCEN       11400.5    
 TAUY      -85239.0    
 AREA       78.5398    
 DF        -2029.54    
 DT        -7500.79    

 *GET  NAME      AS  SSUM       DT         VALUE= -7500.788    

 PARAMETER TORQ1 =     7500.788304    

With the initial guess on the twist we obtained a torque of 7 501 [Nmm] but 10 000 [Nmm] was desired.


5.8) Final load

Return to the solution-module and perform a new analysis with the final load.

To find the final twist the initial twist shall be scaled to compensate for the deviation in torque. Obtained 7 501 [Nmm], desired 10 000 [Nmm]:

θ2 = θ1·10000/7501 = 0.1155549853E-03·10000/7501 = 0.1540571212E-03 [1/mm].

Heat generation pr. unit mass - according to the analogy, as described in § 2, the Heat generation pr. unit mass "q" is simply two times the twist:

q2 = 2·θ2 = 2·0.1540571212E-03 = 0.3081142423E-03.

The heat generation is evenly distributed over cross-section:

fig.5.8.1
Fig. 5.8.1: Final load.


5.9) Finding/checking the final (reaction) torque

The procedure is described in depth in § 4.

Solve the problem and go to the post-processor. The commands are:

 csys,11 $ dsys,11 $ rsys,11          ! cyl. coord. sys. at centre of mass.
 etable,rcen,cent,x                   ! element centre radius.
 etable,tauy,tg,x                     ! therm. gradient, radial component.
 etable,area,volu                     ! volume (depth == 1) eq. area.
 smult,df,area,tauy                   ! tangential force pr. element.
 smult,dt,df,rcen                     ! torque pr. element.
 ssum                                 ! all etable values are summed up.
 *get,name,ssum,dt                    ! torque from sum to parametre "name".
 torq2= -name                         ! torque in the shaft. 

The radial (X-dir.) thermal gradient corresponding to the tangential (Y-dir.) shear stress, that is used in the operations above, look like this:

fig.5.9.1
Fig. 5.9.1: Radial thermal gradient ~ tangential shear stress.

And the output shows:

 ACTIVE COORDINATE SYSTEM SET TO     11  (CYLINDRICAL)

 DISPLAY COORDINATE SYSTEM SET TO     11  (CYLINDRICAL)

 RSYS KEY SET TO  11

 USE COORDINATE SYSTEM  11 FOR SOLUTION RESULTS

 STORE RCEN     FROM ITEM=CENT COMP=X     FOR ALL SELECTED ELEMENTS

 STORE TAUY     FROM ITEM=TG   COMP=X     FOR ALL SELECTED ELEMENTS

 STORE AREA     FROM ITEM=VOLU  FOR ALL SELECTED ELEMENTS

 OPERATION SMUL  RESULT= DF        OPERAND1= AREA      OPERAND2= TAUY    
  FACTOR1=  1.0000       FACTOR2=  1.0000       CONSTANT=  0.0000    

 OPERATION SMUL  RESULT= DT        OPERAND1= DF        OPERAND2= RCEN    
  FACTOR1=  1.0000       FACTOR2=  1.0000       CONSTANT=  0.0000    

 SUM ALL THE ACTIVE ENTRIES IN THE ELEMENT TABLE 

 TABLE LABEL     TOTAL
 RCEN       11400.5    
 TAUY      -113640.    
 AREA       78.5398    
 DF        -2705.77    
 DT        -10000.0    

 *GET  NAME      AS  SSUM       DT         VALUE= -10000.00    

 PARAMETER TORQ2 =     9999.999996    

With the final guess on the twist we are spot-on with the obtained torque!


5.10) Stiffness of the bar

With the general formula for angular deflection of a bar:

θr = (TL)/(KG)

Where:
T ~ Torque.
L ~ Length.
K ~ Stiffness.
G ~ Shear modulus.

By re-arranging the general formula K can be found to:

K = (TL)/(θrG) = T/(θG)

In this example the stiffness "the effective torsional moment of inertia" of the bar of a ¼ circular sector of radius 10 [mm] is calculated to [mm4]:

K = T/(θG) = 10000/(0.1540571212E-03·78846) = 823.261 [mm4]

The analytical solution says K = 824.156 [mm4] a good agreement with the FEA value.


5.11) Stresses

A plot of the temperature will give the "shape" of potential:

fig.5.11.1
Fig. 5.11.1: The potential "temperature".

A plot of the thermal gradient will give the sear stress:

fig.5.11.2
Fig. 5.11.2: The shear stress "thermal gradient".

Max shear stress can be read to: τ = 72.61728 ~ 73. [N/mm2].

The max shear stress can finally be converted to an equivalent stress by the traditionally used relationship, in this case:
σ = τ·√3 = 72.61728·√3 = 125.7768168 ~ 126. [N/mm2].

The analytical solution says τ = 72.6546181092 [N/mm2] a good agreement with the FEA value.

The iteration or better the scaling procedure can be shown as:

Property Symbol Initial Final Formula
Initial stiffness Ip 1097.57
Twist pr. unit length θ 0.115555 E-3 0.154057 E-3
Heat generation pr. unit mass Q 0.231110 E-3 0.308114 E-3
Torque T 7500.78 10000.0 10000.0
Max shear stress τ 54.4687 72.6173 72.6546
Stiffness Keff 823.261 823.261 824.156


Top Up Down Bot.

6) Multiple symmetry

To illustrate the method, I've chosen the somewhat special shape of a straight bar with a cross-section as an equilateral triangle, because there is a nice analytical solution to this problem.


6.1) Geometry

Due to symmetry there is, 3 "teeth" (n=3), but a tooth also have a symmetry line giving an additional 2-fold symmetry. The FEA model will therefore only consist of: ½ "tooth", giving it a 3·2 = 6-fold symmetry:

Dimension: a = 10 [mm].

fig.6.1.1
Fig. 6.1.1: Geometry.

A plot of the area as modelled in ANSYS®:

fig.6.1.2
Fig. 6.1.2: Area.


6.2) Material

The only material property needed is the heat conductivity. According to the analogy, as described in § 2, the heat conductivity "K" is simply the reciprocal of the shear modulus "G".

K = 1/G

If the specific value of the shear modulus isn't known the modulus of elasticity (Young's modulus) together with Poisson's ratio can be used as a basis. This will allow for consistent results if the results are compared with 3-D models.

G = E/(2(1+ν))

In a case with steel and units in N and mm, the values can be:

E = 205.e+3 [N/mm2].

ν = 0.3 [-].

G = E/(2(1+ν)) = (205.e+3)/(2(1+0.3)) = 78846. [N/mm2].

K = 1/G = 1/78846 = 0.1268292683E-04

Now if K is defined to this value in ANSYS®, so that the output shows:

 PARAMETER K =    0.1268292683E-04 

Then the command for the heat conductivity is:

 mp,kxx,1,k 

And the output shows:

 MATERIAL         1     KXX  =  0.1268293E-04 

6.3) Mesh

In this example the area is meshed with the quadrilateral element, no 77 - ANSYS® does a good job in meshing with quad's in 2-D, so why not use this element where ever possible, as quad's in general are better/more efficient than triangles.

fig.6.3.1
Fig. 6.3.1: Mesh.

The shown mesh, with element type 77, has 1703 elements, 5310 nodes and thereby also 5310 degrees of freedom. Though a relatively fine mesh still a small model.


6.4) Preparing load and torque calculation

Before leaving the pre-processor we can as well prepare load and torque calculation by issuing the "asum" command and extract polar moment of inertia. It's natural, and assumed here, that the geometry is modelled so, that the centre of the full geometry is positioned in the centre of the global coordinate system(s). "fold" = 6:

 asum                                 ! sum area properties.
 *get,name,area,,ior,z                ! 1/6 polar moment of inertia.
 ip= name*fold                        ! total polar moment of inertia. 

And the output shows:

 PRINT GEOMETRY ITEMS ASSOCIATED WITH THE CURRENTLY SELECTED AREAS

 *** NOTE ***                            CP =       0.688   TIME= 23:32:39
 Density not associated with all selected areas.  Geometry items are     
 based on a unit density.                                                

   TOTAL NUMBER OF AREAS SELECTED =      1  (OUT OF      1 DEFINED)

   TOTAL SURFACE AREA OF ALL SELECTED AREAS =   7.2169    

   TOTAL VOLUME OF ALL SELECTED AREAS =   7.2169    

   CENTER OF MASS:  XC=  2.4056      YC= 0.83333      ZC=  0.0000    

                     *** MOMENTS OF INERTIA ***
           (BASED ON A UNIT DENSITY AND A UNIT THICKNESS)

         ABOUT ORIGIN     ABOUT CENTER OF MASS    PRINCIPAL
     IXX =    7.5176            2.5059            2.2624    
     IYY =    52.623            10.859            11.102    
     IZZ =    60.141            13.365            13.365    
     IXY =   -13.021            1.4468    
     IYZ =    0.0000            0.0000    
     IZX =    0.0000            0.0000    

     PRINCIPAL ORIENTATION VECTORS (X,Y,Z):
       0.986 -0.166  0.000    0.166  0.986  0.000    0.000  0.000  1.000
       (THXY=   -9.553  THYZ=    0.000  THZX=    0.000)

 *GET  NAME      FROM  AREA  ITEM=IOR  Z     VALUE=  60.1406530    

 PARAMETER IP =     360.8439182    

6.5) Support

The support consists simply of defining a temperature of zero along the outer boundary of the cross-section. There is only this one line on an outer boundary, the two other lines are internal symmetry lines:

fig.6.5.1
Fig. 6.5.1: Support.

The number of active degrees of freedom will be reduced moderately (to 5159) but it's probably without importance.


6.6) Initial load

Let the torque be 1.e+4 = 10 000 [Nmm].

As the formula for angular deflection of a circular bar is:

θr = (TL)/(IpG)

Where:
T ~ Torque.
L ~ Length.
Ip ~ Polar moment of inertia.
G ~ Shear modulus.

Then the expression for angular deflection of a circular bar pr. is unit length "twist" is:

θ = T/(IpG)

With:
T = 1.e+4 = 10000. [Nmm]. Given torque.
Ip = 360.844 [mm4]. Polar moment of inertia as described in § 6.4.
G = 78846. [N/mm2]. Shear modulus of steel as described in § 6.2.

The first guess on twist in the bar is the same as if the bar was circular - this is not the case, but it gives a first guess with some physical sense and a value that can have a reasonable size. Alternatively, a unit value could be used. In the end it should make no difference as the result of the initial load set will be scaled to obtain the correct final result.

The values below are only first guesses and will be designated by index 1.

The angular deflection of a circular steel bar pr. is unit length "twist" is:

θ1 = 10000/(360.844·78846) = 0.3514795785E-03 [1/mm].

Heat generation pr. unit mass - according to the analogy, as described in § 2, the Heat generation pr. unit mass "q" is simply two times the twist:

q1 = 2·θ1 = 2·0.3514795785E-03 = 0.7029591570E-03

The heat generation is evenly distributed over cross-section:

fig.6.6.1
Fig. 6.6.1: Initial load.


6.7) Finding the initial (reaction) torque

The procedure is described in depth in § 4. The coordinate system for calculations is now no 1 - global cylindrical coordinate system in the centre of the full geometry.

Solve the problem and go to the post-processor. The commands are (fold = 6):

 csys,1 $ dsys,1 $ rsys,1        ! cyl. coord. sys. at centre.
 etable,rcen,cent,x              ! element centre radius.
 etable,tauy,tg,x                ! therm. gradient, radial component.
 etable,area,volu                ! volue (depth == 1) eq. area.
 smult,df,area,tauy              ! tangential force pr. element.
 smult,dt,df,rcen                ! torque pr. element.
 ssum                            ! all etable values are summed up.
 *get,name,ssum,dt               ! torque from sum into parametre "name".
 torq1= -name*fold               ! torque for the whole shaft. 

The radial (X-dir.) thermal gradient corresponding to the tangential (Y-dir.) shear stress, that is used in the operations above, look like this:

fig.6.7.1
Fig. 6.7.1: Radial thermal gradient ~ tangential shear stress.

And the output shows:

 ACTIVE COORDINATE SYSTEM SET TO      1  (CYLINDRICAL) Z- CYL. AXIS

 DISPLAY COORDINATE SYSTEM SET TO      1  (CYLINDRICAL)

 RSYS KEY SET TO   1

 USE COORDINATE SYSTEM   1 FOR SOLUTION RESULTS

 STORE RCEN     FROM ITEM=CENT COMP=X     FOR ALL SELECTED ELEMENTS

 STORE TAUY     FROM ITEM=TG   COMP=X     FOR ALL SELECTED ELEMENTS

 STORE AREA     FROM ITEM=VOLU  FOR ALL SELECTED ELEMENTS

 OPERATION SMUL  RESULT= DF        OPERAND1= AREA      OPERAND2= TAUY    
  FACTOR1=  1.0000       FACTOR2=  1.0000       CONSTANT=  0.0000    

 OPERATION SMUL  RESULT= DT        OPERAND1= DF        OPERAND2= RCEN    
  FACTOR1=  1.0000       FACTOR2=  1.0000       CONSTANT=  0.0000    

 SUM ALL THE ACTIVE ENTRIES IN THE ELEMENT TABLE 

 TABLE LABEL     TOTAL
 RCEN       4528.24    
 TAUY      -87181.9    
 AREA       7.21688    
 DF        -372.637    
 DT        -1000.08    

 *GET  NAME      AS  SSUM       DT         VALUE= -1000.077    

 PARAMETER TORQ1 =     6000.461938    

With the initial guess on the twist we obtained a torque of 6 000 [Nmm] but 10 000 [Nmm] was desired.


6.8) Final load

Return to the solution-module and perform a new analysis with the final load.

To find the final twist the initial twist shall be scaled to compensate for the deviation in torque. Obtained 6 000 [Nmm], desired 10 000 [Nmm]:

θ2 = θ1·10000/6 000 = 0.3514795785E-03·10000/6000 = 0.5857542005E-03 [1/mm].

Heat generation pr. unit mass - according to the analogy, as described in § 2, the Heat generation pr. unit mass "q" is simply two times the twist:

q2 = 2·θ2 = 2·0.5857542005E-03 = 0.1171508401E-02.

The heat generation is evenly distributed over cross-section:

fig.6.8.1
Fig. 6.8.1: Final load.


6.9) Finding/checking the final (reaction) torque

The procedure is described in depth in § 4.

Solve the problem and go to the post-processor. The commands are (fold = 6 ):

 csys,1 $ dsys,1 $ rsys,1             ! cyl. coord. sys. at centre.
 etable,rcen,cent,x                   ! element centre radius.
 etable,tauy,tg,x                     ! therm. gradient, radial component.
 etable,area,volu                     ! volue (depth == 1) eq. area.
 smult,df,area,tauy                   ! tangential force pr. element.
 smult,dm,df,rcen                     ! torque pr. element.
 ssum                                 ! all etable values are summed up.
 *get,name,ssum,dt                    ! torque from sum into parametre "name".
 torq2= -name*fold                    ! torque for the whole shaft. 

The radial (X-dir.) thermal gradient corresponding to the tangential (Y-dir.) shear stress, that is used in the operations above, look like this:

fig.6.9.1
Fig. 6.9.1: Radial thermal gradient ~ tangential shear stress.

And the output shows:

 ACTIVE COORDINATE SYSTEM SET TO      1  (CYLINDRICAL) Z- CYL. AXIS

 DISPLAY COORDINATE SYSTEM SET TO      1  (CYLINDRICAL)

 RSYS KEY SET TO   1

 USE COORDINATE SYSTEM   1 FOR SOLUTION RESULTS

 STORE RCEN     FROM ITEM=CENT COMP=X     FOR ALL SELECTED ELEMENTS

 STORE TAUY     FROM ITEM=TG   COMP=X     FOR ALL SELECTED ELEMENTS

 STORE AREA     FROM ITEM=VOLU  FOR ALL SELECTED ELEMENTS

 OPERATION SMUL  RESULT= DF        OPERAND1= AREA      OPERAND2= TAUY    
  FACTOR1=  1.0000       FACTOR2=  1.0000       CONSTANT=  0.0000    

 OPERATION SMUL  RESULT= DT        OPERAND1= DF        OPERAND2= RCEN    
  FACTOR1=  1.0000       FACTOR2=  1.0000       CONSTANT=  0.0000    

 SUM ALL THE ACTIVE ENTRIES IN THE ELEMENT TABLE 

 TABLE LABEL     TOTAL
 RCEN       4528.24    
 TAUY      -145292.    
 AREA       7.21688    
 DF        -621.013    
 DT        -1666.67    

 *GET  NAME      AS  SSUM       DT         VALUE= -1666.667    

 PARAMETER TORQ2 =     9999.999997    

With the final guess on the twist we are spot-on with the obtained torque!


6.10) Stiffness of the bar

With the general formula for angular deflection of a bar:

θr = (TL)/(KG)

Where:
T ~ Torque.
L ~ Length.
K ~ Stiffness.
G ~ Shear modulus.

By re-arranging the general formula K can be found to:

K = (TL)/(θrG) = T/(θG)

In this example the stiffness "the effective torsional moment of inertia" of the bar of a cross-section as an equilateral triangle of side length of 10 [mm] is calculated to [mm4]:

K = T/(θG) = 10000/(0.5857542005E-03·78846) = 216.523 [mm4]

The analytical solution says K = 216.506 [mm4] a very good agreement with the FEA value.


6.11) Stresses

A plot of the temperature will give the "shape" of potential:

fig.6.11.1
Fig. 6.11.1: The potential "temperature".

A plot of the thermal gradient will give the sear stress:

fig.6.11.2
Fig. 6.11.2: The shear stress "thermal gradient".

Max shear stress can be read to: τ = 199.9824 ~ 200. [N/mm2].

The max shear stress can finally be converted to an equivalent stress by the traditionally used relationship, in this case:
σ = τ·√3 = 199.9824·√3 = 346.3797036 ~ 346. [N/mm2].

The analytical solution says τ = 200. [N/mm2] a good agreement with the FEA value.

The stress plot can finally be graphically expanded in ANSYS®, the command used here is:

 /expand,fold,polar,half,0,(360/fold),0 

The stress plot of the full cross-section of the shaft will then look like:

fig.6.11.3
Fig. 6.11.3: The shear stress, full cross-section "thermal gradient".

The iteration or better the scaling procedure can be shown as:

Property Symbol Initial Final Formula
Initial stiffness Ip 360.844
Twist pr. unit length θ 0.351480 E-03 0.585754 E-03
Heat generation pr. unit mass Q 0.702959 E-03 0.117151 E-02
Torque T 6000.46 10000.0 10000.0
Max shear stress τ 119.999 199.982 200.000
Stiffness Keff 216.523 216.523 216.506


Top Up Down Bot.

7) Hollow sections, ex. 1, no symmetry used

To illustrate the method, I've chosen the somewhat special shape of a straight bar with a cross-section as circle with an eccentric circular hole, because there is a nice analytical solution to this problem. The symmetry isn't used to reduce the problem size.


7.1) Geometry

The diameter of the shaft is 4 [mm] and the hole has a diameter of 2 [mm]. The hole is 0.5 [mm] eccentric to the centre of the shaft.

A plot of the areas as modelled in ANSYS®. Both the solid part of the shaft and the hole are modelled - the hole must be meshed too, so it needs to be modelled too (argument in § 7.5):

fig.7.1.1
Fig. 7.1.1: Areas.


7.2) Material

The only material property needed is the heat conductivity. According to the analogy, as described in § 2, the heat conductivity "K" is simply the reciprocal of the shear modulus "G".

K = 1/G

If the specific value of the shear modulus isn't known the modulus of elasticity (Young's modulus) together with Poisson's ratio can be used as a basis. This will allow for consistent results if the results are compared with 3-D models.

G = E/(2(1+ν))

In a case with steel and units in N and mm, the values can be:

E = 205.e+3 [N/mm2].

ν = 0.3 [-].

G = E/(2(1+ν)) = (205.e+3)/(2(1+0.3)) = 78846. [N/mm2].

K = 1/G = 1/78846 = 0.1268292683E-04

Now if K is defined to this value in ANSYS®, so that the output shows:

 PARAMETER K =    0.1268292683E-04 

Then the command for the heat conductivity is:

 mp,kxx,1,k 

And the output shows:

 MATERIAL         1     KXX  =  0.1268293E-04 

7.3) Mesh

In this example the areas are meshed with the quadrilateral element, no 77 - ANSYS® does a good job in meshing with quad's in 2-D, so why not use this element where ever possible, as quad's in general are better/more efficient than triangles. The hole is modelled too, and it must also be meshed (argument in § 7.5):

fig.7.3.1
Fig. 7.3.1: Mesh.

The shown mesh, with element type 77, has 3540 elements, 10719 nodes and thereby also 10719 degrees of freedom. Though a relatively fine mesh still a small model.


7.4) Preparing load and torque calculation

Before leaving the pre-processor we can as well prepare load and torque calculation by issuing the "asum" command and extract polar moment of inertia (about CoG) and the position of CoG and place the coordinate system for (reaction) torque calculation here. Only the solid part of the shaft is analysed here:

 asel,s,area,,1,4
 asum                               ! sum area properties.
 *get,name,area,,imc,z  $ ip= name  ! polar moment of inertia.
 *get,name,area,,cent,x $ xc= name  ! x-pos. of mass centre. 

And the output shows:

 PRINT GEOMETRY ITEMS ASSOCIATED WITH THE CURRENTLY SELECTED AREAS

 *** NOTE ***                            CP =       0.687   TIME= 22:12:13
 Density not associated with all selected areas.  Geometry items are     
 based on a unit density.                                                

   TOTAL NUMBER OF AREAS SELECTED =      4  (OUT OF      5 DEFINED)

   TOTAL SURFACE AREA OF ALL SELECTED AREAS =   9.4248    

   TOTAL VOLUME OF ALL SELECTED AREAS =   9.4248    

   CENTER OF MASS:  XC=-0.16667      YC=-0.91025E-09  ZC=  0.0000    

                     *** MOMENTS OF INERTIA ***
           (BASED ON A UNIT DENSITY AND A UNIT THICKNESS)

         ABOUT ORIGIN     ABOUT CENTER OF MASS    PRINCIPAL
     IXX =    11.781            11.781            11.781    
     IYY =    10.996            10.734            10.734    
     IZZ =    22.776            22.515            22.515    
     IXY =   0.58166E-05       0.58180E-05
     IYZ =    0.0000            0.0000    
     IZX =    0.0000            0.0000    

     PRINCIPAL ORIENTATION VECTORS (X,Y,Z):
       1.000  0.000  0.000    0.000  1.000  0.000    0.000  0.000  1.000
       (THXY=    0.000  THYZ=    0.000  THZX=    0.000)

 *GET  NAME      FROM  AREA  ITEM=IMC  Z     VALUE=  22.5146694    

 PARAMETER IP =     22.51466935    

 *GET  NAME      FROM  AREA  ITEM=CENT X     VALUE=-0.166666669    

 PARAMETER XC =   -0.1666666686    

Afterwards create a local coordinate system in this position - it's a large advantage if this local coordinate system is cylindrical, as shown here, the commands are:

 csys,0 $ dsys,0
 local,12,1,xc,0,0                    ! local coord.sys. at CoG.
 csys,0 $ dsys,0 

And the output shows:

 ACTIVE COORDINATE SYSTEM SET TO      0  (CARTESIAN)  

 DISPLAY COORDINATE SYSTEM SET TO      0  (CARTESIAN)  

 COORDINATE SYSTEM     12 DEFINITION.  TYPE= 1  (CYLINDRICAL)
   XC,YC,ZC=-0.16667      0.0000      0.0000    
   ANGLES=    0.00    0.00    0.00   PARAMETERS=   1.000   1.000
   ORIENTATION=  1.00  0.00  0.00   0.00  1.00  0.00   0.00  0.00  1.00

 ACTIVE COORDINATE SYSTEM SET TO     12  (CYLINDRICAL)

 ACTIVE COORDINATE SYSTEM SET TO      0  (CARTESIAN)  

The centre of the (outer) shaft circle is the global coord.sys. marked by "XYZ". The centre of the hole is the local coord.sys. no 11, to the right. The CoG of (the solid part of) the shaft is the local coord.sys. no 12, to the left:

fig.7.4.1
Fig. 7.4.1: The local coordinate system in CoG.


7.5) Defining the hole boundary condition

While still in the pre-processor, the boundary condition of the hole is defined. It's stated in the collection of formulas that, apart from the thermal analogy, a soap film (or membrane) analogy can be considered. The soap film spans a hole of the contour of the shaft and is subjected to a differential pressure. The height of the bubble corresponds to the potential and the gradient to the stress. The hole in the shaft can then be incorporated as a thin wire bent into the shape of the hole and floating at a constant elevation in the soap film. This explains why the hole must be modelled, meshed and loaded too - otherwise there would be a hole in the soap film. In the thermal analogy this must correspond to all nodes on the hole boundary has the same potential or temperature, thus the temperature must be coupled for these nodes:

fig.7.5.1
Fig. 7.5.1: Coupled degrees of freedom (temperature).


7.6) Support

The support consists simply of defining a temperature of zero along the outer boundary of the cross-section:

fig.7.6.1
Fig. 7.6.1: Support (and cp's).

The number of active degrees of freedom will be reduced moderately, both by this support and the previous coupling of DoFs, (to 10320) but it's probably without importance.


7.7) Initial load

Let the torque be 1.e+3 = 1 000 [Nmm].

As the formula for angular deflection of a circular bar is:

θr = (TL)/(IpG)

Where:
T ~ Torque.
L ~ Length.
Ip ~ Polar moment of inertia.
G ~ Shear modulus.

Then the expression for angular deflection of a circular bar pr. is unit length "twist" is:

θ = T/(IpG)

With:
T = 1.e+3 = 1000. [Nmm]. Given torque.
Ip = 22.5147 [mm4]. Polar moment of inertia as described in § 7.4.
G = 78846. [N/mm2]. Shear modulus of steel as described in § 7.2.

The first guess on twist in the bar is the same as if the bar was circular - this is not the case, but it gives a first guess with some physical sense and a value that can have a reasonable size. Alternatively, a unit value could be used. In the end it should make no difference as the result of the initial load set will be scaled to obtain the correct final result.

The values below are only first guesses and will be designated by index 1.

The angular deflection of a circular steel bar pr. is unit length "twist" is:

θ1 = 1000/(22.5147·78846) = 0.5633183695E-03 [1/mm].

Heat generation pr. unit mass - according to the analogy, as described in § 2, the Heat generation pr. unit mass "q" is simply two times the twist:

q1 = 2·θ1 = 2·0.5633183695E-03 = 0.1126636739E-02

The heat generation is evenly distributed over cross-section:

fig.7.7.1
Fig. 7.7.1: Initial load.


7.8) Finding the initial (reaction) torque

The procedure is described in depth in § 4.

Solve the problem and go to the post-processor.

It's the initial (reaction) torque in the physical cross-section, that should be found, therefore select these elements. It's the thermal gradient in the radial (X) direction that's used:

fig.7.8.1
Fig. 7.8.1: Initial radial gradient.

The commands are:

 csys,12 $ dsys,12 $ rsys,12          ! cyl. coord. sys. at centre of mass.
 etable,rcen,cent,x                   ! element centre radius.
 etable,tauy,tg,x                     ! therm. gradient, radial component.
 etable,area,volu                     ! volue (depth == 1) eq. area.
 smult,df,area,tauy                   ! tangential force pr. element.
 smult,dt,df,rcen                     ! torque pr. element.
 ssum                                 ! all etable values are summed up.
 *get,name,ssum,dt                    ! torque from sum to parametre "name".
 trq1= -name                          ! torque in the shaft. 

And the output shows:

 ACTIVE COORDINATE SYSTEM SET TO     12  (CYLINDRICAL)

 DISPLAY COORDINATE SYSTEM SET TO     12  (CYLINDRICAL)

 RSYS KEY SET TO  12

 USE COORDINATE SYSTEM  12 FOR SOLUTION RESULTS

 STORE RCEN     FROM ITEM=CENT COMP=X     FOR ALL SELECTED ELEMENTS

 STORE TAUY     FROM ITEM=TG   COMP=X     FOR ALL SELECTED ELEMENTS

 STORE AREA     FROM ITEM=VOLU  FOR ALL SELECTED ELEMENTS

 OPERATION SMUL  RESULT= DF        OPERAND1= AREA      OPERAND2= TAUY    
  FACTOR1=  1.0000       FACTOR2=  1.0000       CONSTANT=  0.0000    

 OPERATION SMUL  RESULT= DT        OPERAND1= DF        OPERAND2= RCEN    
  FACTOR1=  1.0000       FACTOR2=  1.0000       CONSTANT=  0.0000    

 SUM ALL THE ACTIVE ENTRIES IN THE ELEMENT TABLE 

 TABLE LABEL     TOTAL
 RCEN       3846.19    
 TAUY      -164990.    
 AREA       9.42478    
 DF        -570.369    
 DT        -948.573    

 *GET  NAME      AS  SSUM       DT         VALUE= -948.5732    

 PARAMETER TRQ1 =     948.5731765    

With the initial guess on the twist we obtained a torque of 948.573 [Nmm] but 1 000 [Nmm] was desired.


7.9) Final load

Return to the solution-module and perform a new analysis with the final load.

To find the final twist the initial twist shall be scaled to compensate for the deviation in torque. Obtained 948.573 [Nmm], desired 1 000 [Nmm]:

θ2 = θ1·1000/948.573 = 0.5633183695E-03·1000/948.573 = 0.5938586326E-03 [1/mm].

Heat generation pr. unit mass - according to the analogy, as described in § 2, the Heat generation pr. unit mass "q" is simply two times the twist:

q2 = 2·θ2 = 2·0.5938586326E-03 = 0.1187717265E-02.

The heat generation is evenly distributed over cross-section:

fig.7.9.1
Fig. 7.9.1: Final load.


7.10) Finding/checking the final (reaction) torque

The procedure is described in depth in § 4.

Solve the problem and go to the post-processor.

It's the final (reaction) torque in the physical cross-section, that should be found, therefore select these elements. It's the thermal gradient in the radial (X) direction that's used:

fig.7.10.1
Fig. 7.10.1: Final radial gradient.

The commands are:

 csys,12 $ dsys,12 $ rsys,12          ! cyl. coord. sys. at centre of mass.
 etable,rcen,cent,x                   ! element centre radius.
 etable,tauy,tg,x                     ! therm. gradient, radial component.
 etable,area,volu                     ! volue (depth == 1) eq. area.
 smult,df,area,tauy                   ! tangential force pr. element.
 smult,dt,df,rcen                     ! torque pr. element.
 ssum                                 ! all etable values are summed up.
 *get,name,ssum,dt                    ! torque from sum to parametre "name".
 trq2= -name                          ! torque in the shaft. 

And the output shows:

 ACTIVE COORDINATE SYSTEM SET TO     12  (CYLINDRICAL)

 DISPLAY COORDINATE SYSTEM SET TO     12  (CYLINDRICAL)

 RSYS KEY SET TO  12

 USE COORDINATE SYSTEM  12 FOR SOLUTION RESULTS

 STORE RCEN     FROM ITEM=CENT COMP=X     FOR ALL SELECTED ELEMENTS

 STORE TAUY     FROM ITEM=TG   COMP=X     FOR ALL SELECTED ELEMENTS

 STORE AREA     FROM ITEM=VOLU  FOR ALL SELECTED ELEMENTS

 OPERATION SMUL  RESULT= DF        OPERAND1= AREA      OPERAND2= TAUY    
  FACTOR1=  1.0000       FACTOR2=  1.0000       CONSTANT=  0.0000    

 OPERATION SMUL  RESULT= DT        OPERAND1= DF        OPERAND2= RCEN    
  FACTOR1=  1.0000       FACTOR2=  1.0000       CONSTANT=  0.0000    

 SUM ALL THE ACTIVE ENTRIES IN THE ELEMENT TABLE 

 TABLE LABEL     TOTAL
 RCEN       3846.19    
 TAUY      -173935.    
 AREA       9.42478    
 DF        -601.291    
 DT        -1000.00    

 *GET  NAME      AS  SSUM       DT         VALUE= -1000.000    

 PARAMETER TRQ2 =     999.9999997    

With the final guess on the twist we are spot-on with the obtained torque!


7.11) Stiffness of the bar

With the general formula for angular deflection of a bar:

θr = (TL)/(KG)

Where:
T ~ Torque.
L ~ Length.
K ~ Stiffness.
G ~ Shear modulus.

By re-arranging the general formula K can be found to:

K = (TL)/(θrG) = T/(θG)

In this example the stiffness "the effective torsional moment of inertia" of the steel bar of a hollow circular cross-section of diameter Ø 4/2 and 0.5 [mm] eccentric is calculated to [mm4]:

K = T/(θG) = 1000/(0.5938586326E-03·78846) = 21.3568 [mm4]

The analytical solution says K = 21.3738 [mm4] a good agreement with the FEA value.


7.12) Stresses

A plot of the temperature will give the "shape" of potential, first the solution of the full problem:

fig.7.12.1
Fig. 7.12.1: The potential "temperature", full problem.

Then a plot of the temperature/potential, of the solid part of the shaft:

fig.7.12.2
Fig. 7.12.2: The potential "temperature", solid part.

A plot of the thermal gradient will give the sear stress:

fig.7.12.3
Fig. 7.12.3: The shear stress "thermal gradient", full problem.

fig.7.12.4
Fig. 7.12.4: The shear stress "thermal gradient", solid part.

Max shear stress can be read to: τ = 124.828 ~ 125. [N/mm2].

The max shear stress can finally be converted to an equivalent stress by the traditionally used relationship, in this case:
σ = τ·√3 = 124.828·√3 = 216.208 ~ 216. [N/mm2].

The analytical solution says τ = 122.287 [N/mm2] an o.k. agreement with the FEA value.

The iteration or better the scaling procedure can be shown as:

Property Symbol Initial Final Formula
Initial stiffness Ip 22.5147
Twist pr. unit length θ 0.563318 E-03 0.593859 E-03
Heat generation pr. unit mass Q 0.112664 E-02 0.118772 E-02
Torque T 948.573 1000.00 1000.00
Max shear stress τ 118.408 124.827 122.287
Stiffness Keff 21.3568 21.3568 21.3738


Top Up Down Bot.

8) Hollow sections, ex. 2, some symmetry used

To illustrate the method, I've also chosen a maybe more common example of a thin-walled rectangular tube. The symmetry is used to ease modelling and reduce problem size.


8.1) Geometry

The rectangular tube has the outer dimension of 200 × 100 [mm] and the wall thickness are respectively 5 and 10 [mm]:

fig.8.1.1
Fig. 8.1.1: Full geometry.

A plot of the areas as modelled in ANSYS®. Both the solid part of the shaft and the hole are modelled - the hole must be meshed too, so it needs to be modelled (argument in § 7.5). The 4-fold symmetry is utilised to limit the model to ¼ of the full geometry:

fig.8.1.2
Fig. 8.1.2: Modelled areas.


8.2) Material

The only material property needed is the heat conductivity. According to the analogy, as described in § 2, the heat conductivity "K" is simply the reciprocal of the shear modulus "G".

K = 1/G

If the specific value of the shear modulus isn't known the modulus of elasticity (Young's modulus) together with Poisson's ratio can be used as a basis. This will allow for consistent results if the results are compared with 3-D models.

G = E/(2(1+ν))

In a case with steel and units in N and mm, the values can be:

E = 205.e+3 [N/mm2].

ν = 0.3 [-].

G = E/(2(1+ν)) = (205.e+3)/(2(1+0.3)) = 78846. [N/mm2].

K = 1/G = 1/78846 = 0.1268292683E-04

Now if K is defined to this value in ANSYS®, so that the output shows:

 PARAMETER K =    0.1268292683E-04 

Then the command for the heat conductivity is:

 mp,kxx,1,k 

And the output shows:

 MATERIAL         1     KXX  =  0.1268293E-04 

8.3) Mesh

In this example the areas are meshed with the quadrilateral element, no 77 - ANSYS® does a good job in meshing with quad's in 2-D, so why not use this element where ever possible, as quad's in general are better/more efficient than triangles. The hole is modelled too, and it must also be meshed (argument in § 7.5):

fig.8.3.1
Fig. 8.3.1: Mesh.

The shown mesh, with element type 77, has 2000 elements, 6161 nodes and thereby also 6161 degrees of freedom. Though a relatively fine mesh still a small model.


8.4) Preparing load and torque calculation

Before leaving the pre-processor we can as well prepare load and torque calculation by issuing the "asum" command and extract polar moment of inertia. It's natural, and assumed here, that the geometry is modelled so, that the centre of the full geometry is positioned in the centre of the global coordinate system(s). "fold" = 4:

 asum                                 ! sum area properties.
 *get,name,area,,ior,z                ! 1/4 polar moment of inertia.
 ip= name*fold                        ! total polar moment of inertia. 

And the output shows:

 PRINT GEOMETRY ITEMS ASSOCIATED WITH THE CURRENTLY SELECTED AREAS

 *** NOTE ***                            CP =       1.856   TIME= 11:06:38
 Density not associated with all selected areas.  Geometry items are     
 based on a unit density.                                                

   TOTAL NUMBER OF AREAS SELECTED =      2  (OUT OF      3 DEFINED)

   TOTAL SURFACE AREA OF ALL SELECTED AREAS =   950.00    

   TOTAL VOLUME OF ALL SELECTED AREAS =   950.00    

   CENTER OF MASS:  XC=  71.316      YC=  35.658      ZC=  0.0000    

                     *** MOMENTS OF INERTIA ***
           (BASED ON A UNIT DENSITY AND A UNIT THICKNESS)

         ABOUT ORIGIN     ABOUT CENTER OF MASS    PRINCIPAL
     IXX =   0.14329E+07       0.22501E+06       0.13251E+06
     IYY =   0.57317E+07       0.90002E+06       0.99252E+06
     IZZ =   0.71646E+07       0.11250E+07       0.11250E+07
     IXY =  -0.21494E+07       0.26645E+06
     IYZ =    0.0000            0.0000    
     IZX =    0.0000            0.0000    

     PRINCIPAL ORIENTATION VECTORS (X,Y,Z):
       0.945 -0.328  0.000    0.328  0.945  0.000    0.000  0.000  1.000
       (THXY=  -19.145  THYZ=    0.000  THZX=    0.000)

 *GET  NAME      FROM  AREA  ITEM=IOR  Z     VALUE=  7164583.33    

 PARAMETER IP =     28658333.33    

8.5) Defining the hole boundary condition

While still in the pre-processor, the boundary condition of the hole is defined. Temperature must be coupled for these nodes (argument in § 7.5):

fig.8.5.1
Fig. 8.5.1: Coupled degrees of freedom (temperature).


8.6) Support

The support consists simply of defining a temperature of zero along the outer boundary of the cross-section:

fig.8.6.1
Fig. 8.6.1: Support (and cp's).

The number of active degrees of freedom will be reduced moderately, both by this support and the previous coupling of DoFs, (to 5920) but it's probably without importance.


8.7) Initial load

Let the torque be 10.e+6 = 10 000 000 [Nmm].

As the formula for angular deflection of a circular bar is:

θr = (TL)/(IpG)

Where:
T ~ Torque.
L ~ Length.
Ip ~ Polar moment of inertia.
G ~ Shear modulus.

Then the expression for angular deflection of a circular bar pr. is unit length "twist" is:

θ = T/(IpG)

With:
T = 10.e+6 = 10 000 000 [Nmm]. Given torque.
Ip = 28 658 333 [mm4]. Polar moment of inertia as described in § 8.4.
G = 78846. [N/mm2]. Shear modulus of steel as described in § 8.2.

The first guess on twist in the bar is the same as if the bar was circular - this is not the case, but it gives a first guess with some physical sense and a value that can have a reasonable size. Alternatively, a unit value could be used. In the end it should make no difference as the result of the initial load set will be scaled to obtain the correct final result.

The values below are only first guesses and will be designated by index 1.

The angular deflection of a circular steel bar pr. is unit length "twist" is:

θ1 = 10000000/(28658333·78846) = 0.4425563302E-05 [1/mm].

Heat generation pr. unit mass - according to the analogy, as described in § 2, the Heat generation pr. unit mass "q" is simply two times the twist:

q1 = 2·θ1 = 2·0.4425563302E-05 = 0.8851126604E-05

The heat generation is evenly distributed over cross-section, incl. the hole:

fig.8.7.1
Fig. 8.7.1: Initial load.


8.8) Finding the initial (reaction) torque

The procedure is described in depth in § 4.

Solve the problem and go to the post-processor.

It's the initial (reaction) torque in the physical cross-section, that should be found, therefore select these elements. It's the radial (X-dir.) thermal gradient corresponding to the tangential (Y-dir.) shear stress that's used:

fig.8.8.1
Fig. 8.8.1: Initial radial gradient.

The commands are (fold = 4):

 csys,1 $ dsys,1 $ rsys,1             ! cyl. coord. sys.
 etable,rcen,cent,x                   ! element centre radius.
 etable,tauy,tg,x                     ! therm. gradient, radial component.
 etable,area,volu                     ! volue (depth == 1) eq. area.
 smult,df,area,tauy                   ! tangential force pr. element.
 smult,dt,df,rcen                     ! torque pr. element.
 ssum                                 ! all etable values are summed up.
 *get,name,ssum,dt                    ! torque from sum to parametre "name".
 trq1= -name*fold                     ! torque in the shaft. 

And the output shows:

 ACTIVE COORDINATE SYSTEM SET TO      1  (CYLINDRICAL) Z- CYL. AXIS

 DISPLAY COORDINATE SYSTEM SET TO      1  (CYLINDRICAL)

 RSYS KEY SET TO   1

 USE COORDINATE SYSTEM   1 FOR SOLUTION RESULTS

 STORE RCEN     FROM ITEM=CENT COMP=X     FOR ALL SELECTED ELEMENTS

 STORE TAUY     FROM ITEM=TG   COMP=X     FOR ALL SELECTED ELEMENTS

 STORE AREA     FROM ITEM=VOLU  FOR ALL SELECTED ELEMENTS

 OPERATION SMUL  RESULT= DF        OPERAND1= AREA      OPERAND2= TAUY    
  FACTOR1=  1.0000       FACTOR2=  1.0000       CONSTANT=  0.0000    

 OPERATION SMUL  RESULT= DT        OPERAND1= DF        OPERAND2= RCEN    
  FACTOR1=  1.0000       FACTOR2=  1.0000       CONSTANT=  0.0000    

 SUM ALL THE ACTIVE ENTRIES IN THE ELEMENT TABLE 

 TABLE LABEL     TOTAL
 RCEN       99006.9    
 TAUY      -20083.6    
 AREA       950.000    
 DF        -15538.7    
 DT       -0.122991E+07

 *GET  NAME      AS  SSUM       DT         VALUE= -1229913.    

 PARAMETER TRQ1 =     4919650.330    

With the initial guess on the twist we obtained a torque of 4 918 274 [Nmm] but 10 000 000 [Nmm] was desired.


8.9) Final load

Return to the solution-module and perform a new analysis with the final load.

To find the final twist the initial twist shall be scaled to compensate for the deviation in torque. Obtained 4 918 274 [Nmm], desired 10 000 000 [Nmm]:

θ2 = θ1·10000000/4918274 = 0.4425563302E-05·10000000/4919650 = 0.8995686695E-05 [1/mm].

Heat generation pr. unit mass - according to the analogy, as described in § 2, the Heat generation pr. unit mass "q" is simply two times the twist:

q2 = 2·θ2 = 2·0.8995686695E-05 = 0.1799137339E-04.

The heat generation is evenly distributed over cross-section:

fig.8.9.1
Fig. 8.9.1: Final load.


8.10) Finding/checking the final (reaction) torque

The procedure is described in depth in § 4.

Solve the problem and go to the post-processor.

It's the initial (reaction) torque in the physical cross-section, that should be found, therefore select these elements. It's the radial (X-dir.) thermal gradient corresponding to the tangential (Y-dir.) shear stress that's used:

fig.8.10.1
Fig. 8.10.1: Final radial gradient.

The commands are (fold = 4):

 csys,1 $ dsys,1 $ rsys,1             ! cyl. coord. sys.
 etable,rcen,cent,x                   ! element centre radius.
 etable,tauy,tg,x                     ! therm. gradient, radial component.
 etable,area,volu                     ! volue (depth == 1) eq. area.
 smult,df,area,tauy                   ! tangential force pr. element.
 smult,dt,df,rcen                     ! torque pr. element.
 ssum                                 ! all etable values are summed up.
 *get,name,ssum,dt                    ! torque from sum to parametre "name".
 trq2= -name*fold                     ! torque in the shaft. 

And the output shows:

 ACTIVE COORDINATE SYSTEM SET TO      1  (CYLINDRICAL) Z- CYL. AXIS

 DISPLAY COORDINATE SYSTEM SET TO      1  (CYLINDRICAL)

 RSYS KEY SET TO   1

 USE COORDINATE SYSTEM   1 FOR SOLUTION RESULTS

 STORE RCEN     FROM ITEM=CENT COMP=X     FOR ALL SELECTED ELEMENTS

 STORE TAUY     FROM ITEM=TG   COMP=X     FOR ALL SELECTED ELEMENTS

 STORE AREA     FROM ITEM=VOLU  FOR ALL SELECTED ELEMENTS

 OPERATION SMUL  RESULT= DF        OPERAND1= AREA      OPERAND2= TAUY    
  FACTOR1=  1.0000       FACTOR2=  1.0000       CONSTANT=  0.0000    

 OPERATION SMUL  RESULT= DT        OPERAND1= DF        OPERAND2= RCEN    
  FACTOR1=  1.0000       FACTOR2=  1.0000       CONSTANT=  0.0000    

 SUM ALL THE ACTIVE ENTRIES IN THE ELEMENT TABLE 

 TABLE LABEL     TOTAL
 RCEN       99006.9    
 TAUY      -40823.2    
 AREA       950.000    
 DF        -31585.1    
 DT       -0.250000E+07

 *GET  NAME      AS  SSUM       DT         VALUE= -2500000.    

 PARAMETER TRQ2 =     9999999.996    

With the final guess on the twist we are spot-on with the obtained torque!


8.11) Stiffness of the bar

With the general formula for angular deflection of a bar:

θr = (TL)/(KG)

Where:
T ~ Torque.
L ~ Length.
K ~ Stiffness.
G ~ Shear modulus.

By re-arranging the general formula K can be found to:

K = (TL)/(θrG) = T/(θG)

In this example the stiffness "the effective torsional moment of inertia" of the rectangular tube with the outer dimension of 200 × 100 [mm] and the wall thickness of respectively 5 and 10 [mm] is calculated to [mm4]:

K = T/(θG) = 10000000/(0.8995686695E-05·78846) = 14 098 898. [mm4]

The analytical solution says K = 13 718 000 [mm4] a quite good agreement with the FEA value.


8.12) Stresses

A plot of the temperature will give the "shape" of potential, first the solution of the full problem:

fig.8.12.1
Fig. 8.12.1: The potential "temperature", full problem.

Then a plot of the temperature/potential, of the solid part of the tube:

fig.8.12.2
Fig. 8.12.2: The potential "temperature", solid part.

A plot of the thermal gradient will give the sear stress:

fig.8.12.3
Fig. 8.12.3: The shear stress "thermal gradient", full problem.

Only stresses in the physical part of the structure is relevant:

fig.8.12.4
Fig. 8.12.4: The shear stress "thermal gradient", solid part.

Max shear stress can be read to: τ = 99.527 ~ 100. [N/mm2] but this value hasn't much meaning as it represents a singularity and will grow to infinite values with increasing mesh density - as predicted by the collection of formulas. The information about the stress concentration located here is however important.

The stress at the centre line of the long, thin side:

fig.8.12.5
Fig. 8.12.5: The shear stress at the centre line of the long, thin side.

The stress at the centre line of the short, thick side:

fig.8.12.6
Fig. 8.12.6: The shear stress "thermal gradient", solid part.

The stress distribution across the centre lines of the two sides are effectively linear and a mean stress can well be found as the average between max and min values:

Long, thin side: τavg = (58.743 + 51.65)/2 = 55.2, formula: 57.0
Short, thick side: τavg = (34.691 + 20.506)/2 = 27.6, formula: 28.5

Good agreement between FEA values and formula.

The iteration or better the scaling procedure can be shown as:

Property Symbol Initial Final Formula
Initial stiffness Ip 28658333.
Twist pr. unit length θ 0.442556 E-05 0.899569 E-05
Heat generation pr. unit mass Q 0.885113 E-05 0.179914 E-04
Torque T 4919650. 10000000. 10000000.
Max shear stress, thin side τa 27.1547 55.1965 56.9801
Max shear stress, thick side τb 13.5775 27.5985 28.4900
Stiffness Keff 14098898. 14098898. 13718000.

The similarity of ideal flow can be recognised in fig. 8.12.1 and fig. 8.12.2 above. As a consequence, the product of mean shear stress (mean flow) and side width must be the same on the two sides (formula values):
τa · ta = 56.9801 · 5 = 284.905
τb · tb = 28.4900 · 10 = 284.900

The same exercise with the FEA values says:
τa · ta = 55.1965 · 5 = 275.983
τb · tb = 27.5985 · 10 = 275.985


8.13) Alternative shape

The singularity in the stress plot in the above section, fig. 8.12.4 can be eliminated by a rather large inner radius. This is only good common sense and mentioned in the collection of formulas.

Shown below an inner corner with a radius of 9 [mm]. A radius of this size moves the location of max stress away from the corner and to the straight, outer side of the thin side:

fig.8.13.1
Fig. 8.13.1: The shear stress "thermal gradient", solid part.

As shown in § 6.11, the plot can be expanded to show the full geometry. The command used for the expansion here is:

 /expand,1,rect,half,0,1.e-4,0,1,rect,half,1.e-4,0,0 

The stress plot of the full cross-section of the tube will then look like:

fig.8.13.2
Fig. 8.13.2: The shear stress "thermal gradient", full cross-section.



Top Up Down Bot.

Overview