Computatuion of v_i (t+ 0.5 *delta t) and omega_i (t+ 0.5 *delta t).

[Srishti-Siddharth-iitb]

I am trying to simulate a single Cosserat rod in MATLAB using your codes and the numerical scheme discussed in the following 2018 paper: https://doi.org/10.1098/rsos.171628. According to equations (C5) and (C6) of the paper (page 25), the position and orientation at (t+delta t) can be computed using v_i (t+ 0.5 * delta t) and omega_i (t+ 0.5 *delta t). Could you please confirm whether these velocities have been computed using the following formulae in your Python code?

v_i (t+ 0.5 delta t) = v_i (t) +0.5 * delta t dv_dt(t+0.5* delta t),

and

omega_i (t+ 0.5 delta t) = omega_i (t) +0.5 * delta t domega_dt(t+0.5 *delta t),

[armantekinalp]

There is a typo there. You can use the above equations.

[Srishti-Siddharth-iitb]

Thank you.

…

On Thu, Jul 25, 2024 at 2:46 AM Arman Tekinalp ***@***.***> wrote: Screen.Shot.2024-07-24.at.4.15.09.PM.png (view on web) <https://github.com/user-attachments/assets/6d51de7c-61e2-4c0a-b7d1-f112345f58b8> There is a typo there. You can use the above equations. — Reply to this email directly, view it on GitHub <#420 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AXRBXIM4EXRH5V5BJAOWSO3ZOAKRVAVCNFSM6AAAAABLNGA43CVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDENBYHEYTENBTHA> . You are receiving this because you authored the thread.Message ID: ***@***.***>

[Srishti-Siddharth-iitb]

Hello, I hope I am not bothering you with the questions about the PyElastica code. I had a doubt regarding the modified equations which you have sent as a screenshot in this email thread. Specifically, to compute $\Dot{v} (t+0.5 \delta t)$ and $\Dot{w} (t+0.5 \delta t)$, one would need the values of $v ( t+0.5 \delta t )$ and $w ( t+0.5 \delta t )$ to compute the discrete versions of Lagrangian transport and unsteady dilation given in Equation (3.9) in https://doi.org/10.1098/rsos.171628. I was wondering whether $v ( t+0.5 \delta t )$ and $w ( t+0.5 \delta t )$ should've computed as follows: 1. $v ( t+0.5 \delta t ) = v(t) + 0.5 \delta t \Dot{v} (t) , and 2. $w ( t+0.5 \delta t ) = w(t) + 0.5 \delta t \Dot{w} (t). Thank you again for all the help. Regards Srishti On Fri, Jul 26, 2024 at 1:48 PM Srishti siddharth < ***@***.***> wrote:

…

Thank you. On Thu, Jul 25, 2024 at 2:46 AM Arman Tekinalp ***@***.***> wrote: > Screen.Shot.2024-07-24.at.4.15.09.PM.png (view on web) > <https://github.com/user-attachments/assets/6d51de7c-61e2-4c0a-b7d1-f112345f58b8> > > There is a typo there. You can use the above equations. > > — > Reply to this email directly, view it on GitHub > <#420 (comment)>, > or unsubscribe > <https://github.com/notifications/unsubscribe-auth/AXRBXIM4EXRH5V5BJAOWSO3ZOAKRVAVCNFSM6AAAAABLNGA43CVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDENBYHEYTENBTHA> > . > You are receiving this because you authored the thread.Message ID: > ***@***.***> >

[sy-cui]

In our code we use $\omega$ at time $t$, although I do think you can use the formulae you proposed. For the time-stepping method we use (Verlet integration) to be applicable, technically $d\omega/dt$ should not be a function of $\omega$. However, the unsteady dilatation and Lagrangian transport terms are rarely important, and the overall algorithm is first-order accurate in time nonetheless.

[armantekinalp]

I am converting this to a discussion since it more about the methods, than code issue.

[Srishti-Siddharth-iitb]

Thank you. I have simulated $\Dot{v}$ and $\Dot{\omega}$ at $t+0.5 \delta t$ using the following:

1. position $r$ and orientation $Q$ at $t+0.5 \delta t$ using the first half step position update.
2. linear velocity $v$ and angular velocity $\omega$ at time t.

However, my results do not match with the one given in the example in https://docs.cosseratrods.org/en/latest/guide/workflow.html. I have added a PDF containing my results and the results obtained from the link. In particular, the displacement of the tip of the rod and the shape of the rod are different for the MATLAB code and PyElastica. I am happy to share my MATLAB code via email.

[Srishti-Siddharth-iitb]

Uploading Pyelatica_MATLAB_results.pdf…

[Srishti-Siddharth-iitb]

In case the link given above doesn't work, please use the following:

https://1drv.ms/b/s!ApBERl2FKoZlpKYatM0DN666Tu6fVQ?e=KXqbeJ

[sy-cui]

I assume you are trying to simulate the deflection of a cantilever beam (in z direction) under a transverse endpoint force (in x direction). Judging from your plots, the boundary condition is likely not implemented correctly. For a clamped BC, you need to fix both the position and the director. It is possible you did not fix the director, resulting in the final configuration being a straight line (Figure 5).

[Srishti-Siddharth-iitb]

Yes, I have implemented the BC for a clamped base. However, the shape of the beam is a straight line. I have attached the snippets of code where BC have been implemented.

[Srishti-Siddharth-iitb]

A is the total number of spatial intervals, a is the current spatial node and k is the current time instant.

[sy-cui]

It seems that your $R$ matrix is set correctly at the boundary. It is not possible to debug without further looking at the full code. You can send the code (well-documented preferably) over and I will take a look.

[Srishti-Siddharth-iitb]

PFA the MATLAB code for the cosserat rod simulation: the main_code.m simulates the rod and plotting_rod.m generates the plots.

Cosserat_rod_simulations.zip

[sy-cui]

Two problems I can see at the moment:

• In the Rodrigues' rotation function you input the rotation angle and the vector, but you only use the direction of the vector, not its magnitude. This causes problems when, e.g. rodrigues_formula(1, w) which does a rotation about w with angle 1.
• The material shear vector is computed as $Q(et - d_3) = Q(et) - e_3$. I believe you put $e_3$ in the parenthesis.

I seems that resolving these problems does not correct the final solution. I suggest you double check all your inputs, because it is very easy to misplace some considering the number of variables you use. Also, consider testing simpler problems first such as the axial stretching case in PyElastica, which involves only stretches and would help you debug incrementally. For reference, see attached a
elastica_simplified.zip
simplified implementation.

[Srishti-Siddharth-iitb]

Thank you. I will take a look at the simplified example.

…

On Sat, 24 Aug, 2024, 12:24 am Songyuan Cui, ***@***.***> wrote: Two problems I can see at the moment: - In the Rodrigues' rotation function you input the rotation angle and the vector, but you only use the direction of the vector, not its magnitude. This causes problems when, e.g. rodrigues_formula(1, w) which does a rotation about w with angle 1. - The material shear vector is computed as $Q(et - d_3) = Q(et) - e_3$. I believe you put $e_3$ in the parenthesis. I seems that resolving these problems does not correct the final solution. I suggest you double check all your inputs, because it is very easy to misplace some considering the number of variables you use. Also, consider testing simpler problems first such as the axial stretching case in PyElastica, which involves only stretches and would help you debug incrementally. For reference, see attached a elastica_simplified.zip <https://github.com/user-attachments/files/16732657/elastica_simplified.zip> simplified implementation. — Reply to this email directly, view it on GitHub <#429 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AXRBXIL7IENSCB66HJKKTJLZS6APVAVCNFSM6AAAAABM4POPKGVHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTANBTGM3TOMA> . You are receiving this because you authored the thread.Message ID: ***@***.***>

[Srishti-Siddharth-iitb]

Hello. I have tried the simplified model which you shared earlier and it's working when axial stretching is applied. However, the code does not run ( the solution is NaN) when I apply the force either in the y-direction or the z- direction. If possible, could you please let me know the reason it is not working? My guess is that is because of the parameters of the beam chosen. I have tried to modify the parameters ( such as increasing the Young's modulus to 1e7) , but have failed to get any valid results. Thank you Srishti On Mon, 26 Aug, 2024, 10:14 am Srishti siddharth, < ***@***.***> wrote:

…

Thank you. I will take a look at the simplified example. On Sat, 24 Aug, 2024, 12:24 am Songyuan Cui, ***@***.***> wrote: > Two problems I can see at the moment: > > - In the Rodrigues' rotation function you input the rotation angle > and the vector, but you only use the direction of the vector, not its > magnitude. This causes problems when, e.g. rodrigues_formula(1, w) > which does a rotation about w with angle 1. > - The material shear vector is computed as $Q(et - d_3) = Q(et) - e_3$. > I believe you put $e_3$ in the parenthesis. > > I seems that resolving these problems does not correct the final > solution. I suggest you double check all your inputs, because it is very > easy to misplace some considering the number of variables you use. Also, > consider testing simpler problems first such as the axial stretching case > in PyElastica, which involves only stretches and would help you debug > incrementally. For reference, see attached a > elastica_simplified.zip > <https://github.com/user-attachments/files/16732657/elastica_simplified.zip> > simplified implementation. > > — > Reply to this email directly, view it on GitHub > <#429 (reply in thread)>, > or unsubscribe > <https://github.com/notifications/unsubscribe-auth/AXRBXIL7IENSCB66HJKKTJLZS6APVAVCNFSM6AAAAABM4POPKGVHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTANBTGM3TOMA> > . > You are receiving this because you authored the thread.Message ID: > ***@***.***> >

[sy-cui]

When bending is involved, the time step restriction is related to the the cross-sectional area. A highly-skewed element aspect ratio (large dx / area) tend to lead to a stiffer problem and a much smaller time step is required. For the parameters used for that particular axial stretching case, you may increase the base radius (or reduce the time steps) to ensure numerical stability.

[Srishti-Siddharth-iitb]

Thank you for the suggestion ! The rod is bending when the radius is increased.

…

On Sun, 15 Sept, 2024, 2:19 am Songyuan Cui, ***@***.***> wrote: When bending is involved, the time step restriction is related to the the cross-sectional area. A highly-skewed element aspect ratio (large dx / area) tend to lead to a stiffer problem and a much smaller time step is required. For the parameters used for that particular axial stretching case, you may increase the base radius (or reduce the time steps) to ensure numerical stability. — Reply to this email directly, view it on GitHub <#429 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AXRBXIJUUK656CR4D4QNGKLZWSOONAVCNFSM6AAAAABM4POPKGVHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTANRUHA4DANQ> . You are receiving this because you authored the thread.Message ID: ***@***.***>

[Srishti-Siddharth-iitb]

Hello, this is Srishti. Sorry for bothering you with another question about the Cosserrat rod equations. I was reading the paper " Forward and inverse problems in the mechanics of soft elements" by Gazzola et. al (2018) to understand the PyElastica codes in detail. Could you please let me know about some references which I could look up to understand how Equations (2.12) - (2.15) were derived? I tried procedure given in the paragraph above these equations and I got some extra terms ( containing \partial e / \ partial t) in (2.14). Thank you Srishti On Mon, 23 Sept, 2024, 2:13 pm Srishti siddharth, < ***@***.***> wrote:

…

Thank you for the suggestion ! The rod is bending when the radius is increased. On Sun, 15 Sept, 2024, 2:19 am Songyuan Cui, ***@***.***> wrote: > When bending is involved, the time step restriction is related to the the > cross-sectional area. A highly-skewed element aspect ratio (large dx / > area) tend to lead to a stiffer problem and a much smaller time step is > required. For the parameters used for that particular axial stretching > case, you may increase the base radius (or reduce the time steps) to ensure > numerical stability. > > — > Reply to this email directly, view it on GitHub > <#429 (reply in thread)>, > or unsubscribe > <https://github.com/notifications/unsubscribe-auth/AXRBXIJUUK656CR4D4QNGKLZWSOONAVCNFSM6AAAAABM4POPKGVHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTANRUHA4DANQ> > . > You are receiving this because you authored the thread.Message ID: > ***@***.***> >

[armantekinalp]

Hi @Srishti-Siddharth-iitb

Thanks for pointing out. So the equation 2.6 is not correct as it is written in the paper. But final equation 2.14 is correct.

You can start your derivation by writing the linear momentum balance for an infinitesimal element.

If you cannot get it let us know.

[Srishti-Siddharth-iitb]

Thank you for the suggestion. I have tried to derive the linear momentum balance equation for an infinitesimal element. I have attached my proof below. The only term in my final equation which does not match Equation (2.14) in the paper is the one containing the external force. Could you please let me know where I went wrong? Also, I have a bit of a confusion about the arc-length parameter. Should we consider the current arc length or the reference arc length while deriving the equations of motion given in the paper? I have read a few papers and textbooks which consider the strain potential energy per unit undeformed length to derive the equations of motion (by taking the variation of the Lagrangian). Thank you Srishti

…

On Fri, Oct 18, 2024 at 2:53 AM Arman Tekinalp ***@***.***> wrote: Hi @Srishti-Siddharth-iitb <https://github.com/Srishti-Siddharth-iitb> Thanks for pointing out. So the equation 2.6 is not correct as it is written in the paper. But final equation 2.14 is correct. You can start your derivation by writing the linear momentum balance for an infinitesimal element. If you cannot get it let us know. — Reply to this email directly, view it on GitHub <#429 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AXRBXIJL67ZQQ6BIB2XY5VLZ4ATGPAVCNFSM6AAAAABM4POPKGVHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTAOJXGY2DAMI> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[Srishti-Siddharth-iitb]

Hello. I have tried to obtain the equations (2.14) and (2.15) given in the paper and have not been able to do so. If possible, could you please share the proof of the above equations? Thank you Srishti On Mon, 21 Oct, 2024, 9:29 pm Srishti siddharth, < ***@***.***> wrote:

…

Thank you for the suggestion. I have tried to derive the linear momentum balance equation for an infinitesimal element. I have attached my proof below. The only term in my final equation which does not match Equation (2.14) in the paper is the one containing the external force. Could you please let me know where I went wrong? Also, I have a bit of a confusion about the arc-length parameter. Should we consider the current arc length or the reference arc length while deriving the equations of motion given in the paper? I have read a few papers and textbooks which consider the strain potential energy per unit undeformed length to derive the equations of motion (by taking the variation of the Lagrangian). Thank you Srishti On Fri, Oct 18, 2024 at 2:53 AM Arman Tekinalp ***@***.***> wrote: > Hi @Srishti-Siddharth-iitb <https://github.com/Srishti-Siddharth-iitb> > > Thanks for pointing out. So the equation 2.6 is not correct as it is > written in the paper. But final equation 2.14 is correct. > > You can start your derivation by writing the linear momentum balance for > an infinitesimal element. > > If you cannot get it let us know. > > — > Reply to this email directly, view it on GitHub > <#429 (reply in thread)>, > or unsubscribe > <https://github.com/notifications/unsubscribe-auth/AXRBXIJL67ZQQ6BIB2XY5VLZ4ATGPAVCNFSM6AAAAABM4POPKGVHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTAOJXGY2DAMI> > . > You are receiving this because you were mentioned.Message ID: > ***@***.***> >