Holzapfel model : fibre orientation

[bharathn]

Dear all, this is a trivial question but nevertheless important.

In the Holzapfel model, the orientations specified are in the global coordinate system? Or in the local? I have a cylindrical ring like structure. There are fibres in the wall. I would like to know whether, when I specifiy 1 0 0 as the fibre direction, they point in the radial direction (denoted by the knot number) or whether they point the general 1 0 0 direction.  Note that it is not a perfect cylinder but a wavy one.

My apologies if the question is not worded properly.

Regards,
Bharath

[Prof. R.L. Taylor]

In ver8.4 the direction is a global Cartesian direction. 

For your model you will have to modify the 'pfiber.f' so the fibers in the reference direction for each element are correct.  If the problem is almost cylindrical you may be able to use the coordinates passed in 'elcoor.h' (reference ones) to define the reference orientation.  Be sure to check that the element you are using computes the xref( values before calling constitution.

You will need something like:

        theta = atan2(xref(2) - xv(2),xref(1) - xv(1))

!       Calculate new a0 for each point
        a(1) =  cos(theta)*av(3) - sin(theta)*av(1)
        a(2) =  sin(theta)*av(3) + cos(theta)*av(1)
        a(3) =  av(2)

!       Normalization
        a0( = a(/sqrt(a(1)**2+a(2)**2+a(3)**2)

[bharathn]

Dear professor,

What if it is not cylindrical and I want to generalize it? Can I use the same piece of code you provided?

[Prof. R.L. Taylor]

You will need to use element to project cylinder basis into element ones,  or if not close to a cylinder whatever medical image you might be using.

[bharathn]

Dear Professor,

I have edited the first few lines of code in pfiber.f to the following:

Code: [Select]

     
subroutine pfiber(d, a0, f, j, sig, dd, ftype, xref, xcur) ! xref and xcur have been added as new inputs!!!

c-----[--+---------+---------+---------+---------+---------+---------+-]
c      Purpose: Fiber models for finite deformation constitution
c               N.B. - Use I_4 not \bar{I_4}
c               1.) Holzapfel-Gasser model

c      Psi(I_4) = 0.5*h1/h2*[exp(h2*(I_4 - 1)**2) - 1)]
c      I_4      = A_I*C_IJ*A_J = (F_iI*A_I)*(F_iJ*A_J) = a_i*a_i
c      a_i      = F_iI*A_I

c      Inputs:
c        d(*)    - Material parameters
c        a0(3)   - Structure vector
c        f(3,3)- Deformation and displacement gradient
c        j       - Jacobian determinant
c        ftype   = 1 uses F; else uses G = displacement gradeint

c      Outputs:
c        sig(6)  - Cauchy stress
c        dd(6,6) - Moduli (current configuration)
c-----[--+---------+---------+---------+---------+---------+---------+-]
      implicit   none

      real*8     d(*), a0(3), f(3,3), j, sig(6), dd(6,6), xref(3)
      real*8 xcur(3)

      integer    ftype, i
      real*8     a(3), tau(6), i4, ee, h1,h2, dpsi1, dpsi2, aLocal(3),
      & a0Local(3), theta

c     Compute current structure tensor
! *************** Am I correct in using xcur and xref? ***************     
      theta = atan2(xref(2) - xcur(2),xref(1) - xcur(1))

!       Calculate new a0 for each point
      aLocal(1) =  cos(theta)*a0(3) - sin(theta)*a0(1)
      aLocal(2) =  sin(theta)*a0(3) + cos(theta)*a0(1)
      aLocal(3) =  a0(2)

!       Normalization
      a0Local(:) = aLocal(:)/sqrt(aLocal(1)**2+aLocal(2)**2+ &
      aLocal(3)**2)
     
      if(ftype.eq.1) then             ! Use deformation gradient
        a(:) = 0.0d0
        do i = 1,3
          a(:) = a(:) + f(:,i)*a0Local(i)
        end do ! i
      else                           ! Use displacement gradient
        a(:) = a0Local(:)
        do i = 1,3
          a(:) = a(:) + f(:,i)*a0Local(i)
        end do ! i
      endif
There are a few syntax errors of course but I would like to know if this is what I am supposed to end up doing.

Regards,
Bharath

[Prof. R.L. Taylor]

Sorry, the 'xv' was the origin for the cylinder (assumed in x-y plane with cylinder along z-axis).  You do not use any xcur, the deformation gradient maps the vector to the current orientation.  You just need the alocal to be in the plane of the cylinder surface.

[bharathn]

Dear professor, I made the changes as you have recommended. Is there any problem I could validate this on?

Code: [Select]

c~        *************** ENSURE CORRECT LOCAL FIBER ORIENTATION ***************     
      theta = atan2(xref(2),xref(1))

      aLoc(1) =  cos(theta)*a0(1) - sin(theta)*a0(2)
      aLoc(2) =  sin(theta)*a0(1) + cos(theta)*a0(2)
      aLoc(3) =  a0(3)


c~        Normalization
      a0Loc(:) = aLoc(:)/sqrt(aLoc(1)**2+aLoc(2)**2+ aLoc(3)**2)
     
      if(ftype.eq.1) then             ! Use deformation gradient
        a(:) = 0.0d0
        do i = 1,3
          a(:) = a(:) + f(:,i)*a0Loc(i)
        end do ! i
      else                           ! Use displacement gradient
        a(:) = a0(:)
        do i = 1,3
          a(:) = a(:) + f(:,i)*a0Loc(i)
        end do ! i
      endif

c~        *************** END ***************     


[Prof. R.L. Taylor]

Try a cylinder with just circumferential ones and then with radial. Finally at an angle

[bharathn]

Dear professor,
I have attached the images of fibre orientations in the 100, 010, and 110 directions. I have also attached the input file. Note that the first digit represents the radial direction and the second, the circumferential.
The 100 orientation seems to not have much of an impact on the expansion of the cylinder. The 010 orientation exhibits the greatest resistance to expansion. The 110 is in between. I believe this should be what we are looking for?

One observation is that the 010 orientation seems to not have a uniform displacement distribution as we move radially, unlike the 110 and the 100 orientations.  Any idea whether this is realistic?

[Prof. R.L. Taylor]

I would expect a circumferential one to still give cylindrical response.  I would output the initial orientation of the fibers at the various gauss points to see what they are.  You should be able to see if they are purely circumferential and of unit length.

[bharathn]

I have outputted the fiber orientations along with the coordinates at which they are located. To me they seem alright.

An inspection of the output file for the circumferential system shows a rather poor convergence even though I have allowed for quite a number of tang/utang iterations. Could this be the issue?

[K.Li]

There might be a problem with the plotting of displacement distribution for a cylinder in FEAP (even if you plot it in a cylindrical coordinate system). You should check your displacement with another software such as Paraview, Tecplot. 

By the way, please have a look at the paper where the model was proposed,

https://www.biomech.tugraz.at/index.php/puplications#2000

The fibers are not supposed to resist the compression load meaning they do not work while under compression.

[bharathn]

Thank you for your suggestion! I even checked the output file manually for the displacements and the circumferential orientation gives me non-cylindrical displacements. Let me check it again!

With regards to not working in compression, in the current case, it is only expansion for the circumferential fibres, isn't it? Or am I missing something? For the radial fibres, they are under compression and therefore do not contribute anything.

[bharathn]

Dear professors,

I have plotted the topmost layers of the cylinder with the circumferential fibres. I used the demo in http://nurbscalculator.in/ . I have attached the NURBS control point data along with the images.

I have compared the output I got from the simulation with the output I would get if it were a 'perfect' circle. Do you think the difference between the two is large enough to see the non-cylindrical plots that we saw in the FEAP window earlier?

Regards,
Bharath

[Prof. R.L. Taylor]

The radial plot seemed to have rather large changes in circumference when you added the fibers.  However, it has a lot of symmetry for only a few elements.  Something seems to be not working well.  I will run your mesh to see if I can spot anything.
Is the correct mesh still the one you posted eariler?