nickabattista/ib2d

An easy to use immersed boundary method in 2D, with full implementations in MATLAB and Python that contains over 75 built-in examples, including multiple options for fiber-structure models and advection-diffusion, Boussinesq approximations, and/or artificial forcing.

ib2d's Issues

possible bugs in data_analysis_in_python

Hi Nick,

It seems there are some bugs in the data_analysis folder for python:

1. Line 126: fX_Lag,fY_Lag,fLagMag,fLagNorm,fLagTan = import_Lagrangian_Force_Data(pathForce,numSim) in "Standard_Data_Analysis_Script.py", fLagNorm and fLagTan are no longer used but returned here.
2. Line 104: U,V = read_Eulerian_Data_From_vtk(path,numSim,strChoice,getxy) in "import_Eulerian_Data.py". Should the return values be like "U,V,x,y"?

Thanks,
Yang

Connect Springs to Immovable Vertices

Hi,
I want to attach a spring between a part of the simulated object and another vertex that only moves when I update its position. I assume that one would use target points for that, but I do not want these points to interact with the fluid in any way. Is there a way to use vertices that do not interact with the fluid or a method to simply apply a force to a vertex of my shape in an arbitrary direction without using additional vertices?
Best regards,
Fabian

Suppress numpy VisibleDepreciationWarning

Something about a non-integer number appears whenever you run the examples...

Inconsistent dimensions between u.0000.vtk and later indices

Steps to reproduce:
Run the MATLAB version of Example_Channel_Flow/Example_Flow_Around_Cylinder. Compare the dimensions of u.0000.vtk and any of the other u.####.vtk files. The DIMENSIONS field is switched.

For some reason, the first vtk output seems to be the transpose of all the later ones. I'm not sure if this issue is confined to this one example (though I'm inclined to doubt it), or if it is only in the MATLAB version. The problem does not seem to be occurring in Omega.0000.vtk, P.0000.vtk, or uMag.0000.vtk

Regarding the example of Skeleton

Hello @nickabattista
I was performing a simulation of the Skeleton example where a rubber band oscillates from a position of non-equilibrium to a position of equilibrium in space and time as the fluid damps the oscillation.
I closed the end point where fluid was being let in. When I ran the simulation, Ib2d threw an error out saying BLOW UP!, Forces are too large.
I thought this is happening because only springs were in the model. Hence I included beams as well. I toggled the values of almost all the variables such as fluid density, viscosity, elastic constants, etc. Nothing really seems to work.
Why is this happening?

Naturally oscillating 2D-cylinder

Kindly include an example of naturally oscillating (elastically mounted / VIV) 2D cylinder.

I think there might be an issue in the IBM_Driver.py?

If you take a look at line 1014 of the IBM_Driver.py, you should see that this line should belong to the if statement above it, not alongside the if statement. Because the variable springs2 will only be created if the variable m is less than 5. If there are 5 columns of data in structure.spring file, among which the last column is degree of non-linearity, the program will report an error, the reason is that springs2 does not exist. Move the line 1014 into the if statement should be able to avoid this issue.
I am not sure whether my findings are correct or not. Hope the author can give me a reply. Thanks.

Heat tube example under asymmetric deformation

Use the example of heat tube to introduce pumping by changing the spring length. Try to introduce a specific deformation only along one side by establishing target points at other side. However the code cannot run under this situation. Would be great if there is any explanations or insight into this problem. Thanks.

Dose the python version stop update?

Excuse me, dose the python version stop update?

Rigid body motion in the python implementation

Dear @nickabattista

I work on unsteady aerodynamics of flapping airfoils. I would like to know if you could share a tutorial / simple test case for a flapping foil with any prescribed kinematic model.

Thanks
Rahul Sundar

Could you share the case of Turok-Hron example using IB2d

Hi Nick,

I saw the work that you shared on Youtube- 'Smaller Deformation Turek-Hron Simulation using IB2d' . It is really amazing! I am trying to do some similar work using IB2d. I will appreciate that if you could share this case. Thanks in advance.

update_nonInv_Beams

Hi Nick,

I tried to simulate the swimmer case with pyIB2d, and I found that a code snippet in IBM_Driver.py:

if update_nonInv_Beams_Flag and nonInv_beams_Yes:
from update_nonInv_Beams import update_nonInv_Beams
#This function is application specific, located with main2d
nonInv_beams_info = update_nonInv_Beams(dt,current_time,beams_info)

Should I use "nonInv_beams_info" to replace "beams_info" in "update_nonInv_Beams" function since it is the nonInv Beam?

Thanks in advance,
Yang

Calculating the non Invariant beam stiffness

Hey @nickabattista
How do you know what beam stiffness to use given the Re input, fluid viscosity and density?

The force on torsional spring doesn't add up?

from the please_find_lagrangnian_force_on_eulerian_grid.py, from line 821-864 is below
id_1 = pts_1[ii] # 1ST Node index
id_2 = pts_2[ii] # (MIDDLE) 2nd Node index -> index that gets force applied to it!
id_3 = pts_3[ii] # 3RD Node index
k_Beam = K_Vec[ii] # Beam stiffness of i-th spring
C = C_Vec[ii] # Curvature of the beam between these three nodes

    Xp = xLag[id_1]          # xPt of 1ST Node Pt. in beam
    Xq = xLag[id_2]          # xPt of 2ND (MIDDLE) Node Pt. in beam
    Xr = xLag[id_3]          # xPt of 3RD Node Pt. in beam

    Yp = yLag[id_1]          # yPt of 1ST Node Pt. in beam
    Yq = yLag[id_2]          # yPt of 2ND (MIDDLE) Node Pt. in beam
    Yr = yLag[id_3]          # yPt of 3RD Node Pt. in beam

    # Checks if Lag. Pts. have passed through the boundary and translates appropriately
    Xp,Xq,Xr = check_If_Beam_Points_Pass_Through_Boundary(ds,Lx,Xp,Xq,Xr)
    Yp,Yq,Yr = check_If_Beam_Points_Pass_Through_Boundary(ds,Ly,Yp,Yq,Yr)

    # Compute Cross-Product
    cross_prod = (Xr-Xq)*(Yq-Yp) - (Yr-Yq)*(Xq-Xp)

    # FORCES FOR LEFT NODE
    bF_x_L = -k_Beam * ( cross_prod - C ) * ( Yr-Yq )
    bF_y_L =  k_Beam * ( cross_prod - C ) * ( Xr-Xq )

    # FORCES FOR MIDDLE NODE
    bF_x_M =  k_Beam * ( cross_prod - C ) * (  (Yq-Yp) + (Yr-Yq) )
    bF_y_M = -k_Beam * ( cross_prod - C ) * (  (Xr-Xq) + (Xq-Xp) )

    # FORCES FOR RIGHT NODE
    bF_x_R = -k_Beam * ( cross_prod - C ) * ( Yq-Yp )
    bF_y_R =  k_Beam * ( cross_prod - C ) * ( Xq-Xp )

    fx[id_1] -= bF_x_L  # Sum total forces for left node,
                                 # in x-direction (this is LEFT node for this beam)
    fy[id_1] -= bF_y_L  # Sum total forces for left node, 
                                 # in y-direction (this is LEFT node for this beam)
    fx[id_2] += bF_x_M  # Sum total forces for middle node,
                                 # in x-direction (this is MIDDLE node for this beam)
    fy[id_2] += bF_y_M  # Sum total forces for middle node, 
                                 # in y-direction (this is MIDDLE node for this beam)
    fx[id_3] -= bF_x_R  # Sum total forces for right node,
                                 # in x-direction (this is RIGHT node for this beam)
    fy[id_3] -= bF_y_R  # Sum total forces for right node, 
                                 # in y-direction (this is RIGHT node for this beam)

it can be seen that bF_x_L + bF_x_M + bF_x_R=0, but it then adds the force to node as -bF_x_L, bF_x_M, -bF_x_R ?
The TOTAL force to add on the three nodes should be zero, but it doesn't here, is it a bug?

Anguilliform Swimmer

Hey @nickabattista
I have looked at the anguilliform swimmer and went through your publications which are based on the parameter subspaces of the simple Anguilliform swimmer. The Re-in and Re-out depend on the length of the swimmer's body in direct proportion. However, upon closely looking, I noticed that the swimmer's length in certain simulations when I change the Re-in is not exactly the same, instead, It changes quite a bit. [I kept a thread on my screen and actually measured]

The code for testing the swimmer's Lagrangian motion, if studied carefully, reveals that all the endpoints of the head lie in the same straight vertical line instead of moving in a rather circular way. This is changing the length significantly. The interpolation is applied to y only. What about x?. I have attached my revised code in my fork which you may please find in my git hub account. The uploading of simulation results is tough because of the size...however the code can be run.

failed to import shared library

Hi Nick,

I am trying to generate shared library to use the C write method on Mac OS using:

python setup.py build_ext —-inplace

A write.so file can be generated, however, when I tried to import write in ipython, an error occurs:

I tried to find a solution online, but it seems a bit complicated. Have you ever met this problem before?

Thanks in advance,
Yang

fluid density problem in input2d file

I am using your ib2d software package for related research work, and I admire you very much for writing such an outstanding program for the study of fluid behavior and making it public. In the process of using, I have a problem and would like to consult you: The fluid density provided in the input2d file is two-dimensional density. Is its value consistent with the density of the material in the three-dimensional state, or does it need to be converted accordingly? Hope to get your reply, thank you

Concentration Dynamics

• Need to reincorporate infrastructure for multiple concentrations, user specified source terms, etc. (circa Oct 2021)
• Current examples are compatible with 2024 IB2d engine; they are the examples from a Aug 2021 commit

Update needed for MATLAB data analysis toolbox after VTK Update from June 2, 2021

The Data Analysis toolbox needs to update some of the "read in" VTK scripts to account for the update to how VTK data is now stored with time information since there is more header information in each .vtk file.

I will fix this by the end of the week (7/16/21).

Turek-hron case Eulerian grid

Hi, the Eulerian grid description in channel_channel.m file is different from that in input2d

mesh point inconsistency

I'm only going to speak to this from the VTK output perspective, not the internal plotting.

Let's suppose you run a simulation with grid parameters as follows (see Flow_Around_Cylinder):
Nx = 512
Ny = 128
Lx = 1.0
Ly = 0.25

For simplicity, let's consider only the x-direction. The important thing here is that the fluid domain is of unit length (periodic) in the x-direction, from 0 to 1.0. However, when the VTK files are written, if you look at the mesh points in the x-dimension, you'll find something surprising: they begin at x=0.0 but end at 0.9980341. That last number is not roundoff error. Here's what I think happened:

• IB2d solved for the fluid velocity at the center of each grid cell.
• The lower-left data point was placed at the Euclidean origin, because that is the lower left corner of all IB2d simulations
• However, because the lower left data point is actually the center of the lower left cell and not it's lower left corner, this actually was a spatial translation of the original domain. The result is that all data points are now misplaced spatially by dx/2 and dy/2.
• The end x-mesh value of 0.9980341 is (within roundoff error) given by 1.0 - 1/512 (dx=1/512). That's because the right-most center cell is (originally) located at 1.0-dx/2, and the shift placed it a further dx/2 to the left.

This is a problem because it means that the fluid data is incorrectly placed in the spatial domain with respect to the original parameters of the simulation. The error is slight and amounts to a straight translation (the domain is now [-0.00098, 0.9990]x[-0.00098,0.246] in this example, instead of [0,1]x[0,0.25]), but in some applications it could matter, particularly where the results are being used for something other than just visualization of the fluid flow.

My suggestion for solving this is two part:

• The VTK writer should be updated to print out the correct spatial mesh information.
• Verify that the correct mesh coordinates are used for internal plotting and fix as necessary. As long as we explicitly specify the spatial limits of the plots, everything should still display with (0,0) in the lower left corner, even if the lower left mesh point is placed at (dx,dy). For quiver plots, everything is straightforward. For pcolormesh etc., some extrapolation to the edge of the domain may be needed, though this is technically interpolation if the periodic boundary condition can be used. It may not be worth fooling with.

Until this is addressed, Planktos will expect this inconsistency and auto-correct it. So please notify me in LOUD terms if this issue is closed!!

issue of convergence study

Hi Nick,

I am trying to reproduce the convergence analysis in your Bioinspiration & Biomimetics paper. I tried three sets of meshes: 256256, 512512, 1024*1024. However, the discrepancies of the results on these three meshes are large. I just changed Nx and Ny for each mesh. I don't know if there are any other parameters need to be changed for different cases, and I hope you can help me on this.

Thanks in advance,
Yang

Trouble with new IBM_Driver

@nickabattista @mountaindust

I had gotten a few simulations with multiple cylinders based off the VIV cylinder example running with the version of the IBM driver from
this commit:

fdbf3f0ffeb2df6d6a8950b9df320999c13c7c5b

But then I merged back in your changes from recently and when I attempt to run the examples, I get this error:

 |****** Prepping Immersed Boundary Simulation ******|


--> Reading input data for simulation...

Index exceeds matrix dimensions.

Error in IBM_Driver (line 98)
boussinesq_Yes = model_Info(20);       % Boussinesq Approx.: 0 (for no) or 1 (for yes)

Error in main2d (line 100)
[X, Y, U, V, xLags, yLags] = IBM_Driver(struct_name, mu, rho, grid_Info, dt, T_final, model_Info);

What more information would you like from me about this?
Here are the examples which run with the old IBM_Driver (on a branch using the old IBM_Driver).