Rotation-free Kirchhoff plate model boundary conditions

[iBatistic]

Hi @Seevanibali ,

in our last talk, you mentioned that the second term (the term with M/D) could be considered as some type of second-order correction:

Today I was thinking about this and realised that if we rearrange the equation above, we can see that we are setting the second-order derivative at the boundary.

Our result is the blue curve, and we should get the red curve, but both curves have the same second-order derivative at the end. Maybe our solver converges to the right value of the second derivative but with wrong displacements.

I guess we can confirm this theory if you visualise the second-order derivative for this case in Abaqus. What do you think about this theory?

[Seevanibali]

@iBatistic Yes, I understand your logic of second-order derivatives, and I think your logic is correct. We are indeed setting the second order derivative. 2nd order derivative represents curvature. I got the following explanation from wiki:

'On the graph of a function, the second derivative corresponds to the curvature or concavity of the graph. The graph of a function with a positive second derivative is upwardly concave, while the graph of a function with a negative second derivative curves in the opposite way.'

Also, we know that when f"(x) > 0, we have minima (for the blue graph, the end point looks like a minima), and f"(x) < 0 is maxima - red graph end point is like maxima. Maybe we need to flip the sign of the second-order derivative? I will check this first thing in the morning. You can also try in your code. Let me know what you think about this.

Also, I can run this cantilever case in Abaqus. I will have to install abaqus in my home machine. Please give me a couple days. For abaqus, I need university ethernet access.

[iBatistic]

Also, we know that when f"(x) > 0, we have minima (for the blue graph, the end point looks like a minima), and f"(x) < 0 is maxima - red graph end point is like maxima. Maybe we need to flip the sign of the second-order derivative? I will check this first thing in the morning. You can also try in your code. Let me know what you think about this.

If I'm not mistaken, the point where the curve turns from concave to convex is also the point where the second derivative equals zero. I am guessing this is the case with the red line representing the physical solution.

Also, I can run this cantilever case in Abaqus. I will have to install abaqus in my home machine. Please give me a couple days. For abaqus, I need university ethernet access

The alternative can be the analytical solution from which we can calculate the second derivative. Do you know what the analytical solution would be for plate loaded with pressure and fixed at end?

[Seevanibali]

Plate_Results_SB.zip

Hi @iBatistic ,

I have attached the results for i. Cantilever plate (one side fixed, rest sides free), ii. Symmetric plate with two sides fixed, and two sides free.

The properties that I have used for the simulation are: Geometry: 10x10x0.1 m Material Properties: E = 200 GPa, nu = 0.3 Mesh: 50x50 Uniform Pressure: 1000 N/m^2

I have attached plots for Displacement_z (U_mag or U3), Rotations (UR1 and UR2) as well as section moments (SM1 and SM2). I am trying to figure out what these section moments actually represent. Also note that UR1 is rotation about x-axis, which I think in our case is thetaVf_y, and UR2 is rotation about y-axis. which is equivalent to thetaVf_x. Also note that, I have used a deformation scale factor of 100 for cantilever plate, and 1000 for fixed-free-symmetric plate.

[iBatistic]

Hi @Seevanibali ,

I’m getting better results, but the values are still slightly overpredicted, as you can see from the images below. Additionally, the boundary values aren’t being displayed correctly, which is affecting the visualization at the edges, as you’re already aware.

[Seevanibali]

Thanks @iBatistic ,

So, Are you saying that you could fix the boundary condition for free edge? Or do you still see the S-shape at the free end?

I will go through tour results in detail tomorrow morning and get back to you with my thoughts.

[iBatistic]

freeEdgeDisplacement.zip

@Seevanibali attached is the boundary condition so you can try. The relaxation factor for clampedMomemnt is 0.01 and for freeEdgeDisplacement is 0.1

I'm not getting S shape, but deflection is positive, is that ok or should it be negative?

[Seevanibali]

@iBatistic ,

Thanks Ivan for sharing the zip file. Deflection positive or negative does not matter. It is just a sign convention issue. Its great that you are not getting S-shape. I will try it out as well, and we could discuss our results then.

[iBatistic]

@Seevanibali, sure. Great

[iBatistic]

@Seevanibali I have a quick question. When you compared the Abaqus results with the OpenFOAM implementation of the Kirchhoff plate, did the solutions match; whether supported or fixed on all edges?

I’m asking because the solutions I obtain for the cantilever plate closely align with those from Euler-Bernoulli beam theory, which makes me wonder if we’re using the right model in Abaqus.

[Seevanibali]

@iBatistic Hey Ivan,

I did not run results this time. But last time, when Philip had run abaqus model, for simply supported and fixed fixed case, it matched the results with OpenFOAM. If you are suspecting error in the Abaqus results, I will rerun the abaqus models tomorrow in my work machine for fixed and simply supported case and let you know.

[iBatistic]

It doesn't necessarily have to be an error, it could be that a different plate model is being used, so the results are slightly different. Just wanted to double-check before we dig any further...

[Seevanibali]

@iBatistic Hey Ivan,

I was running the free BC condition today, and I found something strange when comparing my equation with yours.

Look at the above screenshot. I have two comments, one as 'incorrect version' - which I had in my code, and the other as the 'correct version' - which is how you coded. I honestly still do not see the difference between the two lines. I know that the clubbing of terms is different, but essentially, in my eyes, it is the same equation. Can you help me recognize my mistake here? It looks like some sort of parsing error. But I cannot figure out what. Nevertheless, both codes compile and run. So, I would never have discovered this if I had not compared it with how you coded it.

Anyway, I get the following results for wVf and thetaVf with the correct version.

As you can see above, the thetaVf values are zero at the free end. Also, when you 'WarpbyVector', use 'D' field instead of 'pointD', and the edges of the plate correctly deform. Let me know if this works for you or not.

Now, the interesting point is that If you compare my wVf field with your D_z value that you shared with me a few days ago, the values are so different, so I don't understand what the mistake is here. Can you tell me what mesh discretization you used?

Anyways, I am still looking into my results and trying to figure out what might have gone wrong. Let me know if you have any comments/suggestions

[iBatistic]

@Seevanibali, It seems that in the incorrect version, everything is in the denominator after the / sign. You’ll need an extra ) after (1 - nu).

I’m using a 50x50 mesh, so I’m confident this is not a mesh-related problem. I’m not sure why you're getting a different solution with the correct version.

I’m using FE41. Which version are you using?

[Seevanibali]

@Seevanibali Ivan, I am using OpenFOAM.com version-2312.

Ahh, I see now what you are saying about parenthesis. I has to do with mag(delta) being in denominator.

I will check again tomorrow with fresh mind, and will let you know. Perhaps a short meeting will help. I will text you.

[iBatistic]

Hi @Seevanibali, the only changes I made were to the boundary conditions; the rest of the code is directly from Philip’s development branch, so we should be able to identify the cause. I’m working with fe41, by the way.

Feel free to commit your latest code to your repository, and I’ll try with OpenFOAM 2312.

I agree with the idea of ​​a short meeting. Feel free to suggest a time when you can.

[iBatistic]

@Seevanibali, below are results:	mesh	max w	nCorr
50x50	0.0762266	6049
100x100	0.0760825	13156
200x200	0.0760239	23036

Euler–Bernoulli beam theory max deflection 0.075

[Seevanibali]

Hi @iBatistic, Thanks for the results. I will let you the abaqus results soon

[Seevanibali]

Hi @iBatistic,

I got some results from Abaqus for the Cantilever Plate. I have used Continuum Shell Element, S4R - shell element for thin or thick plates with reduced integration method. I selected plane stress method, with thickness 0.1 m,as we have in our code. Using this method, the following are the value of maximum deflection, UR2 (equivalent to thetaX) and UR1 (equivalent to thetaY):

mesh	max w	UR2	UR1
50x50	0.070509	9.46645e-3	8.42602e-4
100x100	0.0705082	9.46531e-3	8.43772e-4
200x200	0.0705092	9.46517e-3	8.44219e-4

I have a meeting with Philip in a short while. I will ask him what element he thinks is best for abaqus modelling of plates.

I also ran the plate with all boundaries fixed in ABAQUS, and my max deflection and UR1 (UR1 = UR2) for different mesh is as follows:

mesh	max w	UR1=UR2
50x50	6.92311e-4	2.15327e-4
100x100	6.92266e-4	2.15086e-4
200x200	6.92255e-4	2.15075e-4

I also ran - all sides fixed - plate case - in OpenFOAM and see the below results for that:	mesh	max w	UR1=UR2	nCorr
50x50	6.93256e-4	2.14852e-4	3502
100x100	6.91463e-4	2.14964e-4	2224
200x200	6.91014e-4	2.15039e-4	1699

We can see from Abaqus and OpenFOAM results comparison that the results are very close but definitely not exact.

Also, I ran the same cantilever test case using the 2312 version of OpenFOAM, with solution tolerance 1e-8 and maximum correctors as 50000. The 200 by 200 mesh did not converge for that tolerance level, even for 50000 iterations. The residual for 'M' reached 1.5 e-7, and for 'w' 1.32e-8. The max deflection for the other two mesh is:

mesh	max w	nCorr
50x50	0.0752151	6302
100x100	0.0755997	11875
200x200	0.0757968	50000 (not converged to 1e-8 tolerance)

[iBatistic]

Hi @Seevanibali ,

The results for both the fully fixed and supported edge cases are quite consistent between OpenFOAM and Abaqus. In my experience, this level of agreement is about as good as we can expect. Given these results, I'm uncertain whether a model that performs well in these cases could still diverge significantly in the cantilever plate scenario.

Regarding the maximum deflection: our values match. I calculate the maximum deflection by checking its value at the edge directly. Paraview, on the other hand, shows the deflection at the center of the cells along the edge, which aligns with the values you included in the table.

The lack of convergence on the finest mesh may be due to slight variations in how OpenFOAM handles gradients across different versions, but I don’t think that should be our focus right now.

Probably 0.07 is indeed the final result. I’ve noticed that our UR distribution isn’t fully aligned with Abaqus’s, which may be the root of the dissimilarity (take a look at the first image in my fifth post). Specifically, at the edge, the gradient in the tangential direction appears to be zero in our cases and we should use zeroGradient for the tangential derivative (we calculate normal but tangential we should use from cell centre). I'm not sure about this, but it seems to me that something like this can be problem.

You can continue discussing this with Philip, and I’ll try to figure it in the meantime.

[Seevanibali]

Hi @iBatistic,

Philip told me that we cannot compare the results of the Euler-Bernoulli theory with the Kirchoff plate like the way we did because for 2-D Kirchhoff theory, as we both discussed in the morning, variables are calculated along the length of the plate and also across the cross-section. So, Philip wants us to set the top and bottom patches of our plate as 'empty' and have only one element in the thickness direction, and then it is, in fact, like an Euler-Bernoulli beam. He believes that our numerical model should give the exact results of Euler-Bernoulli. For that, we will have to look at the mesh discretization error and see whether the error reduces with second-order accuracy as we refine the mesh.

Then, in the next step, we can have an actual 2-D plate, where we could see or have a comparison with Abaqus and OpenFOAM. But, as a first step, according to Philip, it is necessary to see that our plate model exactly goes to the analytical Euler-Bernoulli result. I am not sure how much of it I could clearly explain. But feel free to give your comments/suggestions.

He was showing me one of his test cases of cantilever beam for the block-coupled solver paper that he has (I think it was 2016). And said we could do a similar approach to that paper. Also, he pointed out that in the max deflection formula delta = wL^4/8EI, the E would be different for plane stress and plane strain. Essentially, for nu = 0, plane stress and plane strain are the same; otherwise, for the plane strain method, we will have to substitute the E value with E/(1-nu^2). Well, he mentioned a lot of things like this. I need to personally look into his paper and then process all his comments. At this moment, I am simply dumping here whatever I can remember from the discussions :D

I actually ran a test case of the cantilever now with 50 elements along the length direction and one element along the thickness direction. I set nu = 0 so that we do not get any thickness effects. I made the top and bottom patches empty for all initial conditions (p, M, and w) and in the `faMeshDefinition'. See below my max deflection image (the max deflection at the face is 0.0750276, error of 0.0368% when compared to 0.075 - so it matches):

No, I have to run the mesh discretization error and check its accuracy.

Regarding your other comments, I also think they are pretty close. But, after discussions with Philip and his comments, I am now trying to think of other alternative ways to compare the two results.

[iBatistic]

Hi @Seevanibali, I ran the case in Ansys and got the same results as you in Abaqus.

In the case of nu = 0, I got 0.075, which is the same as what you got in the 1D case. I also tried setting nu=0 in solids4foam and got 0.075 without reducing it to the 1D problem.

Our problem can be observed in the thetaVf Y field. If you compare this field with the UR field from Abaqus, you’ll notice that thetaVf Y goes to zero along the right edge, rather than displaying a distribution as expected. I'm not sure why this is happening.

[Seevanibali]

Hi @iBatistic,

Can you please check the overleaf you shared for the free BC? I was attempting to understand and re-derive the BC for the w Laplace equation. I have added what I was attempting to derive in Section 3 of the overleaf. Please let me know your comments on it.

[iBatistic]

I’ll take a look over the weekend or by Monday at the latest. I have an idea: we should try running the case with three edges simply supported and one free edge. This might help us rule out the possibility that the issue is related to the clamped boundary condition. Although I doubt that’s the cause, it’s worth testing to be sure.

What puzzles me is that we’re working with a scalar Laplace equation, and for discretization, we only need the normal gradient at the boundary, which is also scalar. I’m not sure how the tangential gradient, which shouldn't be zero in our case, has come into play.

[Seevanibali]

Hi @iBatistic ,

Its a good idea. I will trying running three side simply supported and one side free plate case in both OpenFOAM and Abaqus today, and let you know. I know it's a weekend, so please take your time.

[Seevanibali]

Hi @iBatistic,

To be sure that I am simulating the right model in ABAQUS (I use S4R elements), I ran a fully simple supported plate for 50 by 50 mesh, and the values compare as follows:	variable	ABAQUS	solids4foam
max w	2.233e-3	2.21794e-3
thetaX or thetaY	7.426e-4	7.35443e-4

Now, for our case:

I have attached the ABAQUS results for a plate with simply supported edges for three sides and one side free.

SS_plate_3sides_1side_free.zip

The corresponding plots in OpenFOAM for the same test case (I used a mesh of 50by50 in both cases) are given below:

SS_plate_3sides_1side_free.zip

When you compare the results, you will notice that the values in solids4Foam are overpredicted compared to ABAQUS. Here are some comparisons of max values:		ABAQUS	solids4Foam
w	7.055e-3	9.249e-3
thetaX	1.079e-3	1.1827e-3
thetaY	2.269e-3	2.963e-3

However, the trend for displacements and rotations is similar. For example, see UR1 of ABAQUS

As compared to solids4Foam below:

Maybe we could have a short meeting tomorrow. Let me know if some time suits you.

[iBatistic]

A short meeting is a good idea, later today I will let you know for timings.

[Seevanibali]

Hi @iBatistic,

If you remember that for freeEdgeDisplacement, we used a relaxation factor of 0.1, and for clamped edge 0.001. I was trying to run the simply supported plate with one free edge with different relaxation factors for the free edge, and as you can see from the table below, the results are getting closer to ABAQUS results. In fact, for a factor of 0.001 for free edge, the max values of w, thetaX, and thetaY have dropped lower than what was observed through ABAQUS.

Relaxation Factor	w	thetaX	thetaY	nCorr
0.1	9.249e-3	1.1827e-3	2.963e-3	1601
0.01	8.991e-3	1.17e-3	2.8823e-3	8484
0.001	6.813e-3	1.114e-3	2.1933e-3	13871
ABAQUS results	7.055e-3	1.079e-3	2.269e-3

We can discuss more on this during the meeting.