gibbonCode / GIBBON

The Geometry and Image-Based Bioengineering add-On for MATLAB
http://gibboncode.org/
GNU Affero General Public License v3.0
188 stars 60 forks source link

How to define more than 1 set of mat_axis on each element with fibers #35

Closed danielcpnkimo closed 5 years ago

danielcpnkimo commented 5 years ago

Hi, I'm trying to simulate biaxial extension of orthotropic material matrix with exponential-power law fiber distributed in random direction. When I try the compressible neo-Hookean material with fibers by modifying the DEMO_febio_0015_cube_fibers_transiso, it works. However, when I switch to 'Fung-ortho-compressible' material matrix and assign febio_spec.Material.material{1}.mat_axis.ATTR.type='vector'; febio_spec.Material.material{1}.mat_axis.a=[1 0 0]; febio_spec.Material.material{1}.mat_axis.d=[0 1 0]; , I cannot run the simulation. I've check the usermanual of solid mixture but I didn't see any part that preventing this to work. I'm wondering how to set multiple mat_axis's on each element because my fianl goal is to have two kinds of fibers at each orthogonal element. I need triple mat_axis's on each element. I'd like to know how to do that with GIBBON, thanks.

Kevin-Mattheus-Moerman commented 5 years ago

Thanks for posting this. I think @belgiansteve should update the FEBio user manual on this (e.g. to clarify exactly how to set material axes for material models). You need to do it like below:

<MeshData>
     <ElementData var="mat_axis" elem_set="part3">
          <elem lid="1">
               <a>1,0,0</a>
               <d>0,1,0</d>
          </elem>
          <elem lid="2">
               <a>1,1,0</a>
               <d>0,1,1</d>
          </elem>
     </ElementData>
</MeshData>

I think the demo DEMO_febio_0015_cube_fibers_transiso was already using that. Did you test it with the Fung law?

You may also be interested in the demo: DEMO_febio_0035_blob_shear_contact_hex8. It uses the Fung orthotropic formulation as well (but for a shell element layer). It is still a work in progress but check it out to see if it is useful. Also I just added a function for inverting the Lame constants to get the Fung parameters needed:

[E1,E2,E3,G12,G23,G31,v12,v23,v31]=lameInvertHookeOrthotropic(mu1,mu2,mu3,lambda11,lambda22,lambda33,lambda12,lambda23,lambda13);

Let me know how you get on.

danielcpnkimo commented 5 years ago

I did try to modify DEMO_febio_0015_cube_fibers_transiso with Fung-ortho-compressible material but I think I did it the wrong way by simply changing the material section in the following, %Material section febio_spec.Material.material{1}.ATTR.type='solid mixture'; % febio_spec.Material.material{1}.ATTR.id=1; febio_spec.Material.material{1}.ATTR.id=3;

febio_spec.Material.material{1}.mat_axis.ATTR.type='vector'; febio_spec.Material.material{1}.mat_axis.a=[1 0 0]; febio_spec.Material.material{1}.mat_axis.d=[0 1 0];

febio_spec.Material.material{1}.ATTR.type='Fung-ortho-compressible'; febio_spec.Material.material{1}.ATTR.Name='Plant CMT'; febio_spec.Material.material{1}.E1=124; % inplane strech direction febio_spec.Material.material{1}.E2=124; % inplane side direction febio_spec.Material.material{1}.E3=36; % thickness direction febio_spec.Material.material{1}.G1=67; febio_spec.Material.material{1}.G23=40; febio_spec.Material.material{1}.G31=40; febio_spec.Material.material{1}.v12=-0.075; febio_spec.Material.material{1}.v23=0.87; febio_spec.Material.material{1}.v31=0.26; febio_spec.Material.material{1}.c=1; febio_spec.Material.material{1}.k=120; % Bulk Modulus

but it doesn't work at all.

What is the 'elem_set="part3"' in the above xml code As I said I already have randomly assign fiber direction by mat_axis, and in my file, not the DEMO_febio_0015_cube_fibers_transiso, my output .feb file now look like this.

<MeshData>
    <ElementData elem_set="elementSetTransiso" var="mat_axis">
        <elem lid="1">
            <a>0, 1, 0</a>
            <d>-5.0000000e-01, -0.0000000e+00, 8.6602540e-01</d>
        </elem>
        <elem lid="2">
            <a>0, 1, 0</a>
            <d>-5.0000000e-01, -0.0000000e+00, 8.6602540e-01</d>
        </elem>
        <elem lid="3">
            <a>0, 1, 0</a>
            <d>-5.0000000e-01, -0.0000000e+00, 8.6602540e-01</d>
        </elem>
        .............
        <elem lid="1000">
            <a>0, 1, 0</a>
            <d>-5.0000000e-01, -0.0000000e+00, 8.6602540e-01</d>
        </elem>
    </ElementData>
</MeshData>

What I want is each element has three mat_axis: orthotropic matrix, fiber1, fiber2 So it seemed to me that for my case, each lid's number should have three mat_axis's right? Why is the code above look like each elem lid just has one mat_axis defined by vectors?

danielcpnkimo commented 5 years ago

I checked the DEMO_febio_0035_blob_shear_contact_hex8. It use Ogden material, not Fung material. Is it been updated to use Fung material now? What I see in DEMO_febio_0035_blob_shear_contact_hex8 is the following, %Material section febio_spec.Material.material{1}.ATTR.type='Ogden'; febio_spec.Material.material{1}.ATTR.id=1; febio_spec.Material.material{1}.c1=c1; febio_spec.Material.material{1}.m1=m1; febio_spec.Material.material{1}.c2=c1; febio_spec.Material.material{1}.m2=-m1; febio_spec.Material.material{1}.k=k;

febio_spec.Material.material{2}.ATTR.type='rigid body'; febio_spec.Material.material{2}.ATTR.id=2; febio_spec.Material.material{2}.density=1; febio_spec.Material.material{2}.center_of_mass=center_of_mass_plate;

febio_spec.Material.material{3}.ATTR.type='rigid body'; febio_spec.Material.material{3}.ATTR.id=3; febio_spec.Material.material{3}.density=1; febio_spec.Material.material{3}.center_of_mass=center_of_mass_probe;

danielcpnkimo commented 5 years ago

Hi, any update on this question? I've tried to use the Fung orthotropic in the bar bending case and it works around 2 weeks ago. However, I still don't know how to have more three set of user-defined material axes for the same element parts. Could you tell me how to do that? the example DEMO_febio_0035_blob_shear_contact_hex8 you provide use Ogden material which is isotropic material so it doesn't work for my case.

danielcpnkimo commented 5 years ago

Hi, I've solved this issue now after asking Professor Maas on FEBio Forum. Here is the thread [https://forums.febio.org/showthread.php?1721-How-to-define-more-three-sets-of-mat_axis-on-each-element&p=7476#post7476]