3D ring subjected to internal pressure using CSUR and POLAr

[fethio]

Dear FEAP users,

I am using FEAP to model mechanical behavior of biological tissues. The tissue has fibers in it, and the HOLZapfel model is taken for the passive response. Within the tissue, the fiber orientation at a point, a0, depends on position, x_i. 

In FEAP 8.5, the new STRUcture command in mesh input exactly suits this need (to assign elements with particular fiber orientations). However, I was not able to make it work. Indeed, I verified that the orientation attribute in the HOLZapfel material specification over-ruled the STRUcture assigned directions. The ordering of MATE and STRU commands did not change the result.

In order to show that STRU does not work, I meshed a 3D ring-type structure subjected to internal pressure. The fiber distribution in the STRU command follows a circumferential path (provided in the attached input and include files). In order to create a homogeneous expansion of the ring, I would like to subject it to internal pressure, and look at the deformations. I tried to use CSUR with POLAr option to load the nodes that are located in the inner radius of the ring, but had no luck in that.

So, how can the CSUR command be used to load the ring with internal pressure?

[Prof. R.L. Taylor]

From what I can tell you have some basic mesh errors

1. the nodes on the elements are not numbered by right-hand rule

2. The polar coordinates should be specified in degrees, not radians.  Also do not try to generate more than 90 degrees with any one patch (in 3-d)

[fethio]

Thank you Prof. Taylor,
I will try to fix it and then share the results with you.

[fethio]

Upon following the corrections by Prof.Taylor, I confirm a working copy of the ring subjected to internal pressure loading. The elements are of type HOLZapfel fibers laid onto a NEOHookean matrix.

I also confirm that the STRUcture command is being overriden by the previous fiber directions in the MATErial card. I tried this with two different fiber orientations, one along the x-, and the other along the y- axis. The results for the x- and y- displacements are provided below.

There seems to be a bug related with the STRUcture card sending appropriate directions to elements of the mesh.

A working copy of the input file is attached. The INCLude file struc_v is in attached to the original posting.

[fethio]

The subroutine pfiber in version 8.5 (july 2017 download) calls another subroutine pfiber_vec. We assumed that  pfiber_vec was created to input the fiber orientations from the STRUcture card. These are stored in hr(np(343)). We also assumed that the flag named ngeo in subroutine pfiber was created to reflect STRUcture card input. However, there were two issues in the original source code

• the subroutine pfiber_vec was incomplete, i.e., it contained only a dummy assignment to the fiber orientation vector (as a0 = [1,0,0]).
• the flag ngeo was unresponsive to adding a STRUcture card in the input file. (it's value always remained as 1.)

Thus we had to modify the original source codes of pfiber.f and pfiber_vec.f, in order to transmit the STRUcture card data.

In pfiber_vec.f

Code: [Select]

      include 'pointer.h'
      include 'comblk.h'
 
      a1(1) = hr(np(343)+6*(n-1)+0)
      a1(2) = hr(np(343)+6*(n-1)+1)
      a1(3) = hr(np(343)+6*(n-1)+2)

and in pfiber.f, we added the line

Code: [Select]

        call pfiber_vec(n_el,ii,nv, a0) ! until ngeo=4 works!!
under the condition

Code: [Select]

      if(ngeo.eq.1) then 

Compiling these changes into the executable yielded a perfectly expanding ring of circumferential fiber layout under internal pressure.

[Prof. R.L. Taylor]

 If you check the routine "pmfiber.f" you will find that we provide for (a) Cartesian; (b) cylindrical; and (c) spherical layouts that can be set using the fibers without need to specify individual ones.  However, once you go to general geometry you will need to input individual fibers.

[Prof. R.L. Taylor]

On further look, can you check what happens if you specify your material input for the fibers as

FIBEr HOLZapfel h1 h2 1.0  ! The orientation of 1.0,0,0 is just to pass the error test in inmate.f
FIBEr FILE

That should cause it to read from your fiber files. If you had two fibers the input would be

FIBEr HOLZapfel h1 h2 1.0
FIBEr HOLZapfel h3 h4 1.0
FIBEr FILE

and your file would have 6 columns per row.  I hope this would preclude you having to code the routine too.

Please confirm if it works or not.  These options have not been used on anything like what you are doing!

[fethio]

Dear Professor,

Thank you very much for your feedback. I did not recognize the subroutine pmfiber. I will check it and program our in-house contraction subroutine accordingly.

I would just like to make sure about how the STRUcture card is supposed to be read. I assumed that this card ought to be input during the mesh input phase. However, you suggested that it is input during the MATErial card, following the FIBEr card, as::

Code: [Select]

FIBEr HOLZapfel h1 h2 1.0  ! The orientation of 1.0,0,0 is just to pass the error test in inmate.f
FIBEr FILE
So, does the FILE point to the filename in which the STRUcture card exists?

Another issue is how this would work for a user material fiber. After the user material, one is not supposed to enter any additional input cards, so in that case would it work like this?

Code: [Select]

FIBEr HOLZapfel h1 h2 1.0  ! The orientation of 1.0,0,0 is just to pass the error test in inmate.f
FIBEr FILE
UCONstitutive CONTratile ...
Best regards,

[fethio]

I confirm that the suggestion to input general fiber orientations as

Code: [Select]

FIBEr HOLZapfel h1 h2 1.0  ! The orientation of 1.0,0,0 is just to pass the error test in inmate.f
FIBEr struc_v
where struc_v is the name of the file which contains the element fiber orientations did not produce the expected result (symmetric inflation of the ring in our case).
The result did not change whether struc_v started with the STRUcture command in its first line, or not.

[Prof. R.L. Taylor]

The name is inconsistent probably,  you just say FIBEr FILE then add the STRUct data.  See the attached file