Compile and run
In order to run the first example go into the ulambator folder and type `make tutorial6a`. An executable will be created. Then type `./tutorial6a.o` to execute the program. It needs a subfolder called output to write output data.
Tutorial 6a – Fiber motion
In this tutorial the flow around and transport of solid objects is explained. In contrast to liquid droplets solid particles do not change their shape. As an example we study elongated, fiber like, particles.
/* Ulambator configuration file Date: 04.07.2017 Version 0.71 */File Header
#include <iostream> #include <stdio.h> #include "../source/matblc1.h" #include "../source/bndblock.h" #ifndef _WIN32 #include <sys/time.h> #endif #ifdef _WIN32 #define _USE_MATH_DEFINES // for C++ #include <cmath> #endif #include <time.h> #include <stdio.h> #include <stdlib.h> #include <math.h> int main (int argc, char * const argv[]) { std::cout<< "Welcome to ULAMBATOR++ v.8\n"; std::cout<<"Unsteady LAminar Microfluidic Boundary element Algorithm To Obtain Results\n";Reading problem parameters
We specify a working directory, a folder id for sub-directories and some geometric parameters and the particles y-position and orientation angle.
char* savedir; savedir = new char[120]; sprintf(savedir,"./output");Here beta is the fiber height compared to the channel height, ksi is the fiber length compared to the channel width and L_h is the fiber aspect ratio height over length. Please note that the channel width is 2 as it spans from -1 to 1. Hence ksi = 1 means the fiber length is 2.
{ int folderid = 61; double beta = 0.8; // h/H double ksi = 0.8; // L/L_c double L_h = 8.0; // L/hThe fiber initial angle is given in degrees For to steep inclinations, depending on the confinement, time step and boundary resolution collisions between particle and wall may occur, which will result in an erroneous trajectory. The y-position has to be set such that collision with the lateral walls is avoided.
double angle = 15.0; double ypos = 0;With a special boundary condition the outer flow can be imposed without outer walls, resulting in a truely infinite domain (this is only the switch, see how it is done below).
int switch_to_infinitedomain = 0;Optional reading of input parameters from the command line.
if (argc>4) { folderid = atoi(argv[1]); beta = atof(argv[2]); ksi = atof(argv[3]); ypos = atof(argv[4]); if (argc>5){ angle = atof(argv[5]); } }Computing the channel aspect ratio Lc_H, fiber width, reduced fiber length and gap size. Lc_H is the aspect ratio of the channel and relates a length “1” (channel half width) to the channel height.
double Lc_H = beta*L_h/ksi; // L_c/H double RH = Lc_H/2.0; double Bf = beta/RH; double Lf = 2.0*ksi-Bf; double gap = 0.5*(1.0-beta); std::cout<<"Fiber Simulation, RH "<<RH<<" gap "<<gap<<" length "<<Lf+Bf<<" angle "<<angle<<".\n";Geometry generation
We initialize two boundaries. One for the wall and another for the particle.
BndBlock bem = BndBlock(2);Variables for boundary conditions. Here inflow in the x-direction.
double bcs[] = {1.0, 0.0}; double pos[6]; double W=1.0; if (switch_to_infinitedomain==1) {In Infinite media an empty panel is set and the boundary condition type is ‘5’ – velocity at infinity.
bem.Elements[0].Init(1,10,0); bem.PanelInit(0,0,0,5,bcs,bcs); } else {In a finite channel from x=-5 to 5 and y =-W to W
bem.Elements[0].Init(4, 2000, 0); pos[0] = -5.0; pos[1] = W; pos[2] =-5.0; pos[3]= -W; bem.PanelInit(0,150, 0, 3, bcs, pos); bcs[0] = 0.0; pos[0] = -5.0; pos[1]= -W; pos[2]= 5.0; pos[3] = -W; bem.PanelInit(0,850, 0, 0, bcs, pos); pos[0] = 5.0; pos[1]= -W; pos[2]= 5.0; pos[3] = W; bem.PanelInit(0,150, 0, 2, bcs, pos); pos[0] = 5.0; pos[1]= W; pos[2]= -5.0; pos[3] = W; bem.PanelInit(0,850, 0, 0, bcs, pos); }Option 1 As an example a particle could also be created from panels
/* bem.Elements[1].Init(4, 500, 1); pos[0] = -Lf; pos[1] = -Bf/2; pos[2] =Lf; pos[3]= -Bf/2; bem.PanelInit(1,100, 0, 11, bcs, pos); bcs[0] = 0.0; pos[0] = Lf; pos[1]= -Bf/2; pos[2]= Lf; pos[3] = Bf/2; bem.PanelInit(1,30, 0, 11, bcs, pos); pos[0] = Lf; pos[1]= Bf/2; pos[2]= -Lf; pos[3] = Bf/2; bem.PanelInit(1,100, 0, 11, bcs, pos); pos[0] = -Lf; pos[1]= Bf/2; pos[2]= -Lf; pos[3] = -Bf/2; bem.PanelInit(1,30, 0, 11, bcs, pos); // Reduction of many panels to a single panel bem.Elements[1].PanelSize[0] = bem.Elements[1].getFullSize(); bem.Elements[1].nbPanels = 1; bem.Elements[1].nbPanRes = 1; */Option 2 Particle created from a predefined function
bem.Elements[1].Init(1,810,1); pos[0] = 0.0; pos[1]= ypos; pos[2]= -Lf; pos[3] = Bf/2; pos[4] = 1.0; bem.PanelInit(1,800, 4, 21, bcs, pos);Option 3 Instead of initializing with PanelInit the fber can be read from file bem.EnableWrite(savedir); bem.DropRead(1,10,savedir);
Also a second fiber could be created, which would require to initialize bem with 3 panel objects (line 73).
/* bem.Elements[2].Init(1,300,1); bem.DropRead(2,111,savedir); pos[0] = -0.5; pos[1]= 0.0; pos[2]= -0.125; pos[3] = 0.0; pos[4] = 1.0; */Problem generation
Matrix initialization and creation of a directory in a sub-directory. The matrix initialization takes as arguments the boundary object, the aspect ratio (channel halfwidth over channel height), viscosity ratio (dispersed phase over continuous phase) and capillary number (viscous stress over surface tension). The last two will have no effect and are only there because the code also works for two-phase flows.
MatBlC1 mat(&bem, RH, 0.0, 0.0); mat.EnableWrite(savedir); bem.UpdateCourbure(1); if (folderid>0) { sprintf(savedir,"case%d", folderid); mat.SubFolder(savedir); }The fiber orientation cannot be set at initialization. Instead the fiber is turned from its initially vertical orientation. According to our definition 0 degrees is a horizontal fiber. The Stock() function is called because changes are made in a volatile memory and need to be copied into the actual boundary position memory.
bem.Translate(1, 0.0, 0.0, (90.0-angle)/180.0*M_PI); bem.Stock(1);Log file 0 and 1 is prepared and erased if it exists. The initial geometry is written and the initial time logged to log file 0.
mat.LogPrepare(0, ""); mat.LogPrepare(1, ""); bem.allWrite(); mat.LogRuntime(0);The time step is set and the number of time steps before writting output data is set to 10
double deltaT = 0.1; int stepsT = 10; mat.setTimeStep(deltaT,stepsT);For faster computation the sub-matrix of the outer (fixed) problem is inverted once, which reduces the final equation system.
mat.PreCondense();Memory is allocated to stroe the displacement velocity and the x,y position.
double move[3]; double xcm, ycm; xcm = 0; ycm = 0;A customized log routine (type 2) is called and the data stored in log file 1.
mat.LogCustom(1,2); /* double gc = (1.0-gap)/(1.0-2*gap+4.0/3.0*gap*gap); // gap flow correction */Iterative solution
Running the loop of 10 iterations 20 times.
for (int j=1; j<=20; j++) {This illustrates how the fiber evolution is done. Build matrix and add the force free constraint to the particle. Solve the matrix and advance the particle. Copy volatile memory to fix memory. The particle movement is corrected by a prefactor gc that takes into account the fact the the velocity in the channel center plane is higher than the mean velocity. Details are in the reference at the bottom of this page.
/* mat.Build(); mat.ForceFree(1, 1.0, gap); mat.Solve(); bem.Translate(1, gc*mat.rhs[2*mat.fixsize]*deltaT, gc*mat.rhs[2*mat.fixsize+1]*deltaT, gc*mat.rhs[2*mat.fixsize+2]*deltaT); bem.Stock(1); */This wraps the scheme above into more generic form. A second order in time ‘RunRigidRK2’ scheme is also implemented.
mat.RunRigidRK1(gap, move);Manual Resetting the fiber in the center
bem.getCoM(1, xcm, ycm); bem.Translate(1, -xcm, 0.0, 0.0); bem.Stock(1); bem.Elements[1].Xrelative += xcm;Then write the geometry and print the velocities and y position to the terminal.
bem.allWrite(); std::cout<<"U, V, omega, Y\n"<<move[0]<<" "<<move[1]<<" "<<move[2]<<" "<<ycm<<" "<<xcm<<"\n"; mat.LogCustom(1,2); mat.LogRuntime(0); } mat.LogRuntime(0); mat.DisableWrite(); } End of program
std::cout << "Finished and Exiting\n"; return 0; }Results
The output folder contains one or several case folders. In the case folder of this example (‘case62’) there are log files and geometry files.
‘log0.txt’ contains the computing time (roughly run time times number of cores).
‘log1.txt’ contains the position and orientation and simulated time.
‘bnd0pos.txt’ contains the outer geometry as a collect of point wise x,y coordinates.
‘drop1tsNNN.txt’ contains the particle geometry as a collect of point wise x,y coordinates and the curvature. the file is called drop because it may also represent a drop. The first line contains the time, relative X and relative Y displacement.
The image above shows the evolution of the particle/fiber as it is streamed from left to right. As the fiber approaches the wall it aligns with the wall before turning away from it. Its trajectory describes an oscillatory motion. One notes that the fiber moves out of the channel domain after a few time steps.
This is only apparently so. We reset the fiber to the center after each iteration cycle and store the relative motion. The plotting routine displaces the fiber by the relative X position.