Hyperelastic Pipe

Hi @Disla
I fixed 2 rows of nodes and the case definitely runs longer. What else do you think I could do to prevent it from diverging? Incrementation setting maybe?

I don’t know what problem you are trying to solve. ?¿?¿

Hi @Disla
So I am applying a pressure difference to the solid. What else do you think might be causing instabilities in the calculix side?
Does it have a hard time dealing with oscillations perhaps?(even though this case has a displacement that is pretty small)

Hi mishal49,

-try to build a structured mesh.

-You are in dynamics meaning inertia plays a fundamental role here. How do you apply the pressure it’s important. it is not the same a shock pressure wave than a gentle pressure increase.

-Your boundaries are constrained so there will be reflections. You would probably need to damp the system a little to avoid messing the convergence.

-Solve the structural problem first with some pre-estimated pressure distribution (or constant maybe).

-Evaluate the displacements with your external pressure of 319Pa . Is FSI worth known that? Are you looking for some resonant effect?

-Your material is nonlinear. Does it show nonlinearity in that strain range?. Is Neo Hookean worth known that?

-You are talking about pressure. ¿If your result is pressure, why *CLOAD?. That is a punctual load in N on the nodes.

*CLOAD

Ninterface,1,0.0
Ninterface,2,0.0
Ninterface,3,0.0

All that is mainly to identify if your model is failing because of the FSI coupling procedure or due to a wrong structural part set up.
I would also evaluate the need of introducing complexities on a problem before proceeding.

Dear ccx users.
In the last few weeks I have worked together with Leonid to identify the origin of the detected discrepancies and to clear any doubts about the results provided by ccx for the hyperelastic materials, Neo Hokean and Mooney Rivlin. We have both calculated separately and come to the same conclusion that the results provided by ccx conform to the results predicted by the theory and that the discrepancy was simply due to a couple of missprints in the paper that served as the basis for the verification. Leonid has also provided us with the source file that allows us to calculate the expected result for any Neo Hookean material. Sorry for any inconvinience and I take this opportunity again to invite Leo to follow and participate in the forum. I am sure that many users will appreciate your experience and knowledge.

> restart; with(LinearAlgebra):
> H1:=0.1:W1:=0.1:L:=1:
> Id:=Matrix(3,3,[[1,0,0],[0,1,0],[0,0,1]]):
Neo-Hookean
> W:=mu/2*(I1-3):
> mu:=9e5:
Tensors
> F:=Matrix(3,3,[[1/sqrt(lambda),0,0],[0,1/sqrt(lambda),0],[0,0,lambda]]);B:=evalm(F&*transpose(F)):
> C:=evalm(transpose(F)&*F);

                    [    1                           ]
                    [---------        0          0   ]
                    [      1/2                       ]
                    [lambda                          ]
                    [                                ]
               F := [                 1              ]
                    [    0        ---------      0   ]
                    [                   1/2          ]
                    [             lambda             ]
                    [                                ]
                    [    0            0        lambda]


                       [  1                        ]
                       [------      0          0   ]
                       [lambda                     ]
                       [                           ]
                  C := [            1              ]
                       [  0       ------       0   ]
                       [          lambda           ]
                       [                           ]
                       [                          2]
                       [  0         0       lambda ]

Stress tensor
> T:=evalm(-p*Id+2*diff(W,I1)*B);

                   [     900000.0000                  ]
                   [-p + ----------- ,     0. ,     0.]
                   [       lambda                     ]
                   [                                  ]
              T := [              900000.0000         ]
                   [0. ,     -p + ----------- ,     0.]
                   [                lambda            ]
                   [                                  ]
                   [                                 2]
                   [0. , 0. , -p + 900000.0000 lambda ]


> 
Invariant
> I1:=trace(C):
Lopital constant
> p:=expand(solve(T[2,2]=0,p));

                                  900000.
                             p := -------
                                  lambda

Equation
> eq:=9500=int(int(T[3,3],x2=0..H1/sqrt(lambda)),x3=0..W1/sqrt(lambda)):
> Solution:=L*(fsolve(eq,lambda,1..2)-1);

                       Solution := 0.500000000

> 
> 

3 Likes

Hello @Disla
I have made this case work. The errors were due to poor meshes generated from several softwares.
I worked on this case, and it deforms well up to the point where it is close to contact. The pipe is predicted to oscillate which is why I must apply some kind of reverse force to avoid contact between the two sides of the deformed surface.
image

I tried making node sets on the deformed areas and applied an opposing force to the surface however, this only slows down the deformation but does not initiate backward motion.
image
image

Is there a ramp function perhaps that I could use to apply the force? or something with contact pairs?
Any suggestions would be greatly appreciated.

Capturing oscillations of a hyperelastic pipe- OpenFOAM and CalculiX - Using preCICE - preCICE Forum on Discourse

1 Like

Hello
I tried to apply a load using the AMPLITUDE card as shown below:

*INCLUDE, INPUT=all.msh
*INCLUDE, INPUT=inlet.nam
*INCLUDE, INPUT=outlet.nam
*INCLUDE, INPUT=wall.nam
*INCLUDE, INPUT=interface.nam
*INCLUDE, INPUT=contact1.nam
*INCLUDE, INPUT=contact2.nam
*MATERIAL, Name=EL
*HYPERELASTIC,NEO HOOKE
244898, 2.5e-7
**0.2449, 0.25
**ELASTIC
**1440000.0, 0.47
*DENSITY
970.0
*SOLID SECTION, Elset=Eall, Material=EL
*AMPLITUDE,NAME=A1
0.0,-1.E-4,0.004,-5.E-4,0.008,-1.E-3,0.012,-3.E-3,
0.016,-6E-3,0.02,-8E-3
*AMPLITUDE,NAME=A2
0.0,1.E-4,0.004,5.E-4,0.008,1.E-3,0.012,3.E-3,
0.016,6E-3,0.02,8E-3
*STEP, INC=1000000
*DYNAMIC,DIRECT,NLGEOM
0.0001,0.1
*RESTART,WRITE,FREQUENCY=1
*BOUNDARY
Ninlet,1,3,0
Noutlet,1,3,0
*CLOAD
Ninterface,1,0.0
Ninterface,2,0.0
Ninterface,3,0.0
*CLOAD,AMPLITUDE=A1
Ncontact1, 1, 1
Ncontact1, 2, 1
Ncontact1, 3, 1
*CLOAD,AMPLITUDE=A2
Ncontact2, 1, 1
Ncontact2, 2, 1
Ncontact2, 3, 1
*NODE FILE
U
*EL FILE
S, E
*END STEP

However, I get the following error:

Congratulations!!.

Walls touching and so close to a clamped end?
It looks suspicious to me. Have you estimate with non linear how much vacuum you need to do that?

Is your scale, right?
imagen

If I remember right your pipe was like 10mm diameter. According to the picture your walls are far from touching.

Regarding the ramp you could try a sin() or cos() function. They work well : F(t)= 1+cos(pi()*t-pi())

EDIT: Check you are using the same units in both solvers (Pa,m,N…).

Yes, scale is right. This is not the total displacement because I have mapped it several times in the fluid side hence, it is only magnitude of displacement from the most recent mapping.
Pipe is 15 mm, inner diameter.
The displacement location does agree with experimental data but how would I do the nonlinear estimate that you have mentioned?

For the ramp function, do I need to define the variables in a separate file or something? Do I just put this as is under the amplitude card? Sorry, I tried looking in the manual but did not find this.
I do have forces that the fluid applies on the solid in a separate .cvg file.

Nonlinear quasi static. Increasing negative value up to the point where walls touch.
Seems a lot of vacuum to close a pipe so close to the clamped area. I would expect to happen at the center. ¿Is this some kind of implosion?. If that agree with the experimental data then you need to add self-contact on the internal wall. Why do you want to avoid it with an oposite force?

Just an amplitude time table.
This is normalized so you can scale the output. Time goes from 0 to 1. Yours should be lower.
1/2*(1+cos(pi()*t-pi()))

*AMPLITUDE,NAME=A_1
0,0
0.0066,0.0001074761413335
0.0133,0.0004363950856248
0.02,0.0009866357858642
0.0267,0.001757954468813
0.0334,0.002750009417013
.
.
.
.

0.9535,0.9946743431235
0.9613,0.9963091478024
0.969,0.9976307011032
0.9768,0.9986725338368
0.9845,0.999407324011
0.9923,0.9998537149224
1,1

The pipe is supposed to oscillate without contact. Because of the many mappings done, a contact model is necessary to nudge the motion in the right direction otherwise, the pipe just continues to collapse inward.

Ok got it. I tried the table like you suggested and first copied the one you have given just to check if that fixes the issue however, looks like it still gives the same error.

The node that it is giving error for is shared between the interface as well as the contact nodesets. Is it not possible to apply 2 different forces?
I initially tried to apply single values to the contact nodesets and the case ran but diverged in a few timesteps so I would assume that 2 forces is not the real issue like it suggests…

Why don’t you try to load at the center of the pipe, or maybe at one end but relaxing the clamped edge where you are imposing the load.?

Setting a frictionless support free to “expand”,… it would become a symmetry condition and less stiff.

It seems hard to make converge close to that clamped area.

If doing that, keep in mind buckling load depends on the unstiffened length.

If your load is just to initiate the oscillation, give a small vacuum pressure and slightly ovalize the shape to promote the buckling shape you are interested in.

Oh, I see. So, by small vacuum pressure, you mean applying probably a hydrostatic pressure on that area in the opposing direction instead of a force?

Applying a load at the center would be a great way however, the pressure difference from the fluid naturally affects that part of the pipe near the clamped edge so probably that would not be feasible. The shape it gets is exactly as expected. The oscillation part is the one that I need to make work.

What a pity. It definitely helps. I think the problem is probably there.
Maybe some more refinement then on that area.
I would work the uncoupled nonlinear dynamic problem to see how far you can push the vacuum before going to FSI. I have also found that reflection easily mess the convergence after a few oscillations.

ezgif-5-3553cdb91b

1 Like

Okay, let me try this.
Refine the mesh near the deforming edge and then apply an amplitude force you mean?
That looks nice. Could I please get your .inp file for reference?
What exactly do you mean by ‘giving a small vacuum pressure’?

The resulting shape of your pipe after some iterations looks like it has buckle because it is under compression. Your fluid pressure inside the pipe looks like points inwards (vaccuum/sucction).
Am I right? Was that your intention?

Refine and apply a negative pressure. Seems the easier way to apply a fluid load (perpendicular to the surface).
The amplitude will be different depending on the type of analysis. For quasistatic you could gently load with the suggested cos(). The aim is to get an idea of the maximum vacuum/pressure your configuration can support. For Dynamic I would try an initial ramp with final quick release if you want to see the pipe oscillation. If the release is too slow the pipe would accommodate and relax without oscillation.

1 Like