Examples of openVFIFE
GinkGo Lv1
 2021-06-13 00:41:02    

Examples

openVFIFE is quite simple and intuitive. The basic procedures of a calculation is as follows:

  1. create a StructSystem to manage all struct objects
  2. create Particles
  3. create Materials and Sections (necessary in bar structures)
  4. create Elements
  5. add constraints
  6. add external force and solve

the pseudocodes for the procesdure

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
#include "structsystem.h"	// using openVFIFE

StructSystem system = StructSyste(id);	// create a StructSystem to manage all struct objects
// ... system setting
Partice* p = new Partice(id, x, y, z);  // create particle object
system.addPartice(p);				   // add particle into system
BaseMaterial mat = LinearElastic(id);   // create Material
// ... material setting
BaseSection sect = Rectangle(id, width, height);    // create section
BaseElement* e = new Beam3D(id, p1, p2, mat, sect);	// create element
system.addElement(e);							  // add element into system
system.addConstraint(pid, c);					  // add constraints

// time history solve
for (int i = 0; i < nstep; i++)
{
    if (i == 0)
    {
        // solve first time step
		system.addExternalForce(pid, force);		// add external force
		system.solve(dt, zeta, true);
		system.setInternalForce();
    }
    else
    {
		system.solve(h, zeta);
		system.clearParticleForce();
		system.setExternalForce(pid, f);
		system.setInternalForce();
    }

    // save results
    system.saveParticleResult(path);
	system.saveElementResult(path);
	system.saveSupportReact(path);
}

system.releaseContainers();					    // relase system containers

As you can see in the above codes, using openVFIFE to conduct structure analysis is same as other FEM software. To illustrate the using of openVFIFE, several examples are list here. You can learn openVFIFE through these simple examples.

Now, let's ROCK!

Example 1

In the first example, I want to show you the moving trajectory calculation of a single paritcle, which is in horizontal projectile motion. The mass of the particle is \(10kg\), the initial speed is \(v_0 = 10m/s\), and the acceleration of gravity \(g = 9.8 m/s\). The governing equation of particle is \[ Horizontal: m\ddot{x} + c\dot{x} = 0\\\\ Vertical: m\ddot{y} + c\dot{y} = mg \] where, \(c = m \zeta\), \(\zeta\) is the damping parameters in openVFIFE. The solution is

if \(c \neq 0\):

\[ x = \frac{mv_0}{c} [1 - e^{-\frac{c}{m}t}]\\\\ v_x = \dot{x} = v_0 e^{\frac{c}{m}t}\\\\ a_x = \ddot{x} = -\frac{c}{m}v_0 e^{-\frac{c}{m}t}\\\\ y = \frac{mg}{c}t + \frac{m^2 g}{c^2}[e^{-\frac{c}{m}t}-1]\\\\ v_y = \dot{y} = \frac{mg}{c}[1 - e^{-\frac{c}{m}t}]\\\\ a_y = \ddot{y} = g e^{-\frac{c}{m}t} \]

else if \(c=0\): \[ x = v_0 t\\\\ v_x = \dot{x} = v_0\\\\ a_x = \ddot{x} = 0\\\\ y = \frac{1}{2} g t^2\\\\ v_y = \dot{y} = gt\\\\ a_y = \ddot{y} = g \]

In openVFIFE, you can solve the problem by the following code(see ./examples/example1/example1.cpp)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
#include "../../headers/structsystem.h"
using namespace std;

int main()
{
    // create system
    StructSystem system = StructSystem(1);
    system.setJobname("Example1");

    double zeta = 0.0;
    cout << "Please enter a damping coefficient(>=0): " << endl;
    cin >> zeta;

    // solve parameters
    double endTime = 10.0;  // total time, s
    double h = 1.0e-3;      // time step, s
    double v0 = 10;         // initial velocity, m/s
    double mass = 10;       // mass, kg
    double g = 9.8;         // acceleration of gravity, m/s-2

    // create particles
    Particle* p1 = new Particle(1, 0, 0);
    system.addParticle(p1);
    // set mass properties
    StdArray6d m = {mass, mass, 0, 0, 0, 0};
    p1->setLumpedMass(m);
    // activate dof of particle
    p1->activateDof("Ux");
    p1->activateDof("Uy");

    string path = system.workdir() + "/" + system.jobname() + "/model";
    system.saveModel(path);

    int nStep = ceil(endTime / h);
    int interval = ceil(0.1/h);
    for (int i = 0; i <= nStep; i++)
    {
        if (i == 0)
        {
            // set initial velocity and acclerate
            StdArray6d v = {v0, 0, 0, 0, 0, 0};
            p1->setVelocity(v);
            system.setAccelerate(0, 9.8, 0);
            system.solve(h, zeta, true);
        }
        else
        {
            system.solve(h, zeta);
            system.clearParticleForce();
            system.setAccelerate(0, g, 0); // add external force
            system.setInternalForce();
        }

        // save results
        if (i % interval == 0)
        {
            string path = system.workdir() + "/" + system.jobname() + "/" +
                          to_string(i*h);
            system.saveParticleResult(path);
        }
    }

    system.releaseContainers();
    return 0;
}

and compile this code by

1
2
#g++ -static -O3 xx.cpp -L /directory/of/libvfife.a -lvfife -o xx.out
g++ -static -O3 projectile.cpp -L ../../library/static -lvfife -o projectile.out

and run it

1
2
3
4
5
6
7
8
9
#time ./xx.out
time ./projectile.out
# you will see the following information in terminal
########### StructSystem Established! ##########
#Please enter a damping coefficient(>0): 
#0
#10000
#100
#./projectile.out  0.00s user 0.01s system 0% cpu 1.392 total

Finally the results of openVFIFE and theroy solution have been compared in the next figures.

​ (a) (b)

​ (c) (d)

​ (e) (f)

Fig. 1. Results, (a) \(x\), (b) \(y\), (c) \(\dot{x}\), (d) \(\dot{y}\), (e) \(\ddot{x}\) (f) \(\ddot{y}\)

when \(\zeta = \frac{c}{m}\) changes from 0.0 to 1.0, openVFIFE always get the right answer. All these examples are trying to tell you how to use openVFIFE, and illustrate as a verification as well. So there will be no detailed commnets on these examples. If you are interested in these examples, just email me.

Example 2

Consider a single bar structure(as shown in the following figure) , the length \(l = 10m\), the section area is \(1m^2\), Young's modulus \(E=10^6Pa\) , mass density \(\rho = 10 kg/m^3\). In this case, I want to show you how to use openVFIFE to solve static problems, and explain the influence of damping coefficient and time step.

Fig.2. You Can Read It

source code, see ./examples/example2/example2.cpp

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
#include "../../headers/structsystem.h"

using namespace std;


int main()
{
    // create system
    StructSystem system = StructSystem(1);
    system.setJobname("example2");

    // setting damping coefficient
    double zeta = 0.5;
    cout << "Please enter a damping coefficient(>=0): " << endl;
    cin >> zeta;

    // setting time step
    double h = 0.0;
    cout << "Please enter a time step(>=0): " << endl;
    cin >> h;

    // solve parameters
    double endTime = 50.0;
    double p = 100;

    // create particles
    Particle * p1 = new Particle(1, 0, 0);
    Particle * p2 = new Particle(2, 10, 0);
    system.addParticle(p1);
    system.addParticle(p2);

    // create material
    LinearElastic mat = LinearElastic(1);
    mat.setE(1.0E6);
    mat.setDensity(10);

    // create section
    CustomSectionParas paras = {.A=1, .Sy=0, .Sz=0, .SHy=0, .SHz=0,
                                .CGy=0, .CGz=0.0, .Iyy=10, .Izz=10, .Iyz=10};
    CustomSection sect = CustomSection(1, paras);

    // create elements
    Link2D* e = new Link2D(1, p1, p2, &mat, &sect);
    system.addElement(e);

    // constraints, Ux, Uy, Uz, Rotx, Roty, Rotz
    DOF c1 = {.key = {true, true, true, true, true, true},
              .val = {0, 0, 0, 0, 0, 0}};
    DOF c2 = {.key = {false, true, true, true, true, true},
              .val = {0, 0, 0, 0, 0, 0}};
    system.addConstraint(1, c1);
    system.addConstraint(2, c2);

    // save model
    string path = system.workdir() + "/" + system.jobname() + "/model";
    system.saveModel(path);
    system.info();      // print structsystem information

    StdArray6d f1 = {100, 0, 0, 0, 0, 0};	//(x, y, z, mx, my, mz)
    int nStep = ceil(endTime/h);
    cout << nStep << endl;
    int interval = ceil(0.1/h);
    for (int i = 0; i <= nStep; i++)
    {
        if (i == 0)
        {
            system.addExternalForce(2, f1);            // add external force
            system.solve(h, zeta, true);
            system.setInternalForce();
        }
        else
        {
            system.solve(h, zeta);
            system.clearParticleForce();
            system.setExternalForce(2, f1);
            system.setInternalForce();
        }

        // save results
        if (i % interval == 0)
        {
            string path = system.workdir() + "/" + system.jobname() + "/" +
                          to_string(i*h);
            system.saveParticleResult(path);
            system.saveElementResult(path);
            system.saveSupportReact(path);
        }
    }

    system.releaseContainers();
    return 0;
}

when \(\Delta t = 0.01\), the axail force of the bar

\(\zeta\) Target (N) openVFIFE
0.0 100.0 --
0.5 100.0 100.0
1.0 100.0 100.0

Fig.3. Axial Force of The Bar

if \(\zeta= 0.0\), the axial force of the bar will not converge (umdamped free vibration); and when \(\zeta\) increased, the answer will finally converge to the static solution, and the larger \(\zeta\) is, the faster convergence is. Besides, \(\zeta\) won't affect the static solution. Hence, for static problems, large \(\zeta\) could be adopted to help converge.

The critical time step \[ \Delta t_{critical} = \frac{2l}{v_c} = \frac{2l}{\sqrt{E/\rho}}=0.0316s \] when \(\zeta= 0.5\), the axail force of the bar

\(\Delta t\) Target (N) openVFIFE
0.001 100.0 100.0
0.01 100.0 100.0
0.1 100.0 nan
0.5 100.0 nan

Fig.4. Axial Force of The Bar, (a) \(\Delta t = 0.001, 0.01\), (b) \(\Delta t = 0.001, 0.01, 0.1, 0.5\),

If \(\Delta t > \Delta t_{critical}\), the solution will unconverge, hence it is important to assign time step. openVFIFE provides a autoTimeStep() function for setting proper time step automaticly.

Example 3

Ref: Ting, E.C., Y. F. Duan, T.Y. Wu. Vector Mechanics of Structure. BeiJing: Science Press, 2012. Chapter 5, example 1

Fig.5. copy from the reference textbook

Parameters \[ a = 10 m, b = 5 m;\\ P_{\alpha} = 20 N, P_{\beta} = 20N;\\ E_{\alpha} = 20000 Pa, E_{\beta} = 10000 Pa. \] Codes: (see ./examples/example3/example3.cpp)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
#include "../../headers/structsystem.h"
using namespace std;

int main()
{
    // create system
    StructSystem system = StructSystem(1);
    system.setJobname("example3");

    // solve parameters
    double endTime = 50.0;
    double zeta = 0.5;
    double h = 0.01;
    double p = 100;

    // create particles
    Particle* p1 = new Particle(1, 0, 0);
    Particle* p2 = new Particle(2, 10, 0);
    Particle* p3 = new Particle(3, 10, 5);
    system.addParticle(p1);
    system.addParticle(p2);
    system.addParticle(p3);
    StdArray6d mass = {10.0, 10.0, 0, 0, 0, 0};
    p2->setLumpedMass(mass);

    // create material
    LinearElastic mat1 = LinearElastic(1);
    mat1.setE(2.0E4);
    LinearElastic mat2 = LinearElastic(2);
    mat2.setE(1.0E4);

    // create section
    CustomSectionParas paras = {.A=1, .Sy=0, .Sz=0, .SHy=0, .SHz=0,
                                .CGy=0, .CGz=0.0, .Iyy=10, .Izz=10, .Iyz=10};
    CustomSection sect = CustomSection(1, paras);

    // create elements
    Link2D* e1 = new Link2D(1, p1, p2, &mat1, &sect);
    Link2D* e2 = new Link2D(2, p2, p3, &mat2, &sect);
    system.addElement(e1);
    system.addElement(e2);

    // constraints
    DOF c1 = {.key = {true, true, true, true, true, true},
              .val = {0, 0, 0, 0, 0, 0}};
    DOF c3 = {.key = {true, true, true, true, true, true},
              .val = {0, 0, 0, 0, 0, 0}};
    system.addConstraint(1, c1);
    system.addConstraint(3, c3);

    // save model
    string path = system.workdir() + "/" + system.jobname() + "/model";
    system.saveModel(path);
    system.info();      // print structsystem information

    StdArray6d f1 = {-20, 0, 0, 0, 0, 0};
    StdArray6d f2 = {-20, 20, 0, 0, 0, 0};
    StdArray6d f3 = {20, 0, 0, 0, 0, 0};

    int nStep = ceil(endTime/h);
    cout << nStep << endl;
    int interval = ceil(0.1/h);
    for (int i = 0; i <= nStep; i++)
    {
        if (i == 0)
        {
            // add external force
            system.addExternalForce(1, f1);
            system.addExternalForce(2, f2);
            system.addExternalForce(3, f3);
            system.solve(h, zeta, true);
            system.setInternalForce();
        }
        else
        {
            system.solve(h, zeta);
            system.clearParticleForce();
            // add external force
            system.addExternalForce(1, f1);
            system.addExternalForce(2, f2);
            system.addExternalForce(3, f3);
            system.setInternalForce();
        }

        // save results
        if (i % interval == 0)
        {
            string path = system.workdir() + "/" + system.jobname() + "/" +
                          to_string(i*h);
            system.saveParticleResult(path);
            system.saveElementResult(path);
            system.saveSupportReact(path);
        }
    }

    system.releaseContainers();
    return 0;
}

Answer comparison

Descriptions Items Target openVFIFE
support reactions \(f_{1x} (N)\) -20 -20
\(f_{1y} (N)\) 20 20
\(f_{3x} (N)\) -20 -20
\(f_{3y} (N)\) 20 20
displacement of particle 2 \(U_x (m)\) 0.01 0.00999503
\(U_y (m)\) -0.01 -0.00999004
element axial force \(F_{Na} (N)\) 20 20
\(F_{Nb} (N)\) 20 20

Fig.6. time history of displacement of particle 2

You can compare the results with the refence. I don't want to do it.

Example 4

Ref: Ting, E.C., Y. F. Duan, T.Y. Wu. Vector Mechanics of Structure. BeiJing: Science Press, 2012. Chapter 6, example 2

Fig.7. coordinate system and discrete scheme

Parameters: \[ b = 10m,\\ m_2 = 6.25 kg, m_4 = m_5 = m_6 = 2.5 kg,\\ E_{1} = 1e5Pa, E_2 = 1e4Pa,\\ P_2 = 1.0N, P_v = 0.1N. \] \(P_2\) is always perpendicular to bar 1, and \(P_v\) is always keep vertical.

codes, (see ./examples/example4/example4.cpp)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
#include "../../headers/structsystem.h"

using namespace std;


int main()
{
    // create system
    StructSystem system = StructSystem(1);
    system.setJobname("example4");

    // solve parameters
    double endTime = 50.0;  // total time, s
    double zeta = 0.5;      // damping coefficient
    double h = 0.01;        // time step, s

    // create particles
    Particle* p1 = new Particle(1, 0.0, 10.0);
    Particle* p2 = new Particle(2, 0.0, 0.0);
    Particle* p3 = new Particle(3, 10.0, 0.0);
    Particle* p4 = new Particle(4, 2.5, 0.0);
    Particle* p5 = new Particle(5, 5.0, 0.0);
    Particle* p6 = new Particle(6, 7.5, 0.0);
    system.addParticle(p1);
    system.addParticle(p2);
    system.addParticle(p3);
    system.addParticle(p4);
    system.addParticle(p5);
    system.addParticle(p6);

    StdArray6d mass2 = {6.25, 6.25, 0, 0, 0, 0};
    p2->setLumpedMass(mass2);
    StdArray6d mass4 = {2.5, 2.5, 0, 0, 0, 0};
    p4->setLumpedMass(mass4);
    p5->setLumpedMass(mass4);
    p6->setLumpedMass(mass4);

    // create material
    LinearElastic mat1 = LinearElastic(1);
    mat1.setE(1.0E5);
    LinearElastic mat2 = LinearElastic(2);
    mat2.setE(1.0E4);

    // create section
    CustomSectionParas paras = {.A=1, .Sy=0, .Sz=0, .SHy=0, .SHz=0,
                                .CGy=0, .CGz=0.0, .Iyy=10, .Izz=10, .Iyz=10};
    CustomSection sect = CustomSection(1, paras);

    // create elements, flexible link
    Link2DLD* e1 = new Link2DLD(1, p1, p2, &mat1, &sect);
    Link2DLD* e2 = new Link2DLD(2, p2, p4, &mat2, &sect);
    Link2DLD* e3 = new Link2DLD(3, p4, p5, &mat2, &sect);
    Link2DLD* e4 = new Link2DLD(4, p5, p6, &mat2, &sect);
    Link2DLD* e5 = new Link2DLD(5, p6, p3, &mat2, &sect);
    system.addElement(e1);
    system.addElement(e2);
    system.addElement(e3);
    system.addElement(e4);
    system.addElement(e5);

    // constraints
    DOF c1 = {.key = {true, true, true, true, true, true},
              .val = {0, 0, 0, 0, 0, 0}};
    DOF c3 = {.key = {true, true, true, true, true, true},
              .val = {0, 0, 0, 0, 0, 0}};
    system.addConstraint(1, c1);
    system.addConstraint(3, c3);

    // save model
    string path = system.workdir() + "/" + system.jobname() + "/model";
    system.saveModel(path);
    system.info();      // print structsystem information

    // get the direction vector of bar1
    const Eigen::Vector3d* ex = e1->ex();
    StdArray6d P2 = {0, 0, 0, 0, 0, 0};
    StdArray6d Pv = {0, -0.1, 0, 0, 0, 0};

    int nStep = ceil(endTime/h);
    cout << nStep << endl;
    int interval = ceil(0.1/h);
    for (int i = 0; i <= nStep; i++)
    {
        // update P2
        P2[0] = 10 * (*ex)(0);
        P2[1] = -10 * (*ex)(1);
        if (i == 0)
        {
            // add external force
            system.addExternalForce(2, P2);
            system.addExternalForce(4, Pv);
            system.addExternalForce(5, Pv);
            system.addExternalForce(6, Pv);
            system.solve(h, zeta, true);
            system.setInternalForce();
        }
        else
        {
            system.solve(h, zeta);
            system.clearParticleForce();
            // add external force
            system.addExternalForce(2, P2);
            system.addExternalForce(4, Pv);
            system.addExternalForce(5, Pv);
            system.addExternalForce(6, Pv);
            system.setInternalForce();
        }

        // save results
        if (i % interval == 0)
        {
            string path = system.workdir() + "/" + system.jobname() + "/" +
                          to_string(i);
            system.saveModel(path);
            system.saveParticleResult(path);
            system.saveElementResult(path);
            system.saveSupportReact(path);
        }
    }

    system.releaseContainers();
    return 0;
}

Fig.8. movement of the structure

Example 5

In example 1, a planar three-bar truss subjected to a load \(P\) in the \(y\) is analyzed, as shown in Fig.9. The cross-sectional area (\(A\)) of bars is \(1m^2\); the length (\(L_{BD}\) is \(1m\); \(\ang ADB = \ang CDB = 45^ \circ\). different plastic material models are considered: an ideal elastic-plastic model and an elastic linear hardening model, as shown in Fig.11. The density (\(\rho\)) of the bars is \(7850kg\). Young's modulus \(E\) is \(206GPa\) in the elastic state; the tangent modulus $Et $ is \(20.6GPa\) in the plastic state; and the yield stress (\(\sigma _y\)) is \(235MPa\).

Fig.9. example 5

codes, (see ./examples/example5/example5.cpp)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
#include "../../headers/structsystem.h"
using namespace std;


int main()
{
    // command-line parameters
    cout << "Please enter P: ";
    double load_factor;
    cin >> load_factor;

    // create system
    StructSystem system = StructSystem(1);
    system.setJobname("BilinearTruss");

    // solve parameters
    double endTime = 100.0;
    double h = 1.0e-3;
    double zeta = 1;
    double length = 1;
    double theta = 45.0 / 180 * PI;

    // create particles
    Particle* p1 = new Particle(1, -length, 0);
    Particle* p2 = new Particle(2, 0, 0);
    Particle* p3 = new Particle(3, length, 0);
    Particle* p4 = new Particle(4, 0, length);
    system.addParticle(p1);
    system.addParticle(p2);
    system.addParticle(p3);
    system.addParticle(p4);

    // create material
    double E = 2.06e11;
    double yield_stress = 2.35e8;
    double Et = 2.06e10;
    double m = 0;
    UniBilinear mat = UniBilinear(1, E, Et, yield_stress, m);
    mat.setDensity(7850);
    // LinearElastic mat = LinearElastic(1);
    // mat.setE(E);
    // mat.setDensity(7850);

    // create section
    CustomSectionParas paras = {.A=1.0, .Sy=0, .Sz=0, .SHy=0, .SHz=0, .CGy=0,
                                .CGz=0.0, .Iyy=0.01, .Izz=0.01, .Iyz=0};
    CustomSection sect = CustomSection(1, paras);

    // create elements
    Link2DLD* e1 = new Link2DLD(1, p1, p4, &mat, &sect);
    Link2DLD* e2 = new Link2DLD(2, p2, p4, &mat, &sect);
    Link2DLD* e3 = new Link2DLD(3, p3, p4, &mat, &sect);
    system.addElement(e1);
    system.addElement(e2);
    system.addElement(e3);

    // constraints
    DOF c = {.key = {true, true, false, false, false, false},
              .val = {0, 0, 0, 0, 0, 0}};
    for (int i = 1; i <= 3; i++)
    {
        system.addConstraint(i, c);
    }
    system.info();

    // save model information
    string path = system.workdir() + "/" + system.jobname() + "/model";
    system.saveModel(path);

    StdArray6d f {};
    // f[1] = load_factor;
    system.addExternalForce(4, f);
    system.solve(h, zeta, true);
    system.setInternalForce();

    int nStep = ceil(endTime / h);
    int interval = ceil(1 / h);
    for (int j = 0; j <= nStep; j++)
    {
        system.solve(h, zeta, false);
        system.clearParticleForce();
        f[1] =  load_factor * (j * 1.0 / nStep);
        system.setExternalForce(4, f);
        system.setInternalForce();

        // save results
        if (j % interval == 0)
        {
            string path = system.workdir() + "/" + system.jobname() + "/" +
                          to_string(f[1]);
            system.saveParticleResult(path);
            system.saveElementResult(path);
            system.saveSupportReact(path);
        }
    }

    system.releaseContainers();
    return 0;
}

Fig.10. results comparison

Example 6

To illustrate the capability of VFIFE in geometric nonlinear analysis, and to test the Link3DLD element, a 24-member shallow dome (as shown in Fig.13) is analyzed by openVFIFE. The topological relationship between elements and particles is depicted in Fig.11, as well as geometry information. The density \(\rho\) of the bars is \(20lb/in^3\); Young’s modulus \(E\) is \(10^6ksi\); and the cross-sectional area (A) of the bars is \(0.1in^2\). A concentrated force \(P\) is imposed on node 1 in the \(z\) direction.

Fig.11. structures

codes, (see ./examples/example6/example6.cpp)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
#include "../../headers/structsystem.h"

using namespace std;

struct Point
{
    int id;
    double x, y, z;
};

struct Element
{
    int id;
    int mat;
    int typ;
    int sec;
    int p1, p2;
};

void importNodes(const string &fname, vector<Point> &v1)
{
    ifstream fin;
    fin.open(fname.c_str());
    if (!fin.is_open())
    {
        cerr << fname << ": open file failed" << endl;
        exit(-1);
    }

    string line;
    getline(fin, line);

    char comma;
    while (true)
    {
        Point p;
        fin >> p.id >> comma >> p.x >> comma >> p.y >> comma >> p.z;
        if(!fin.good()) break;
        v1.push_back(p);
    }
    fin.close();
}

void importElements(const string &fname, vector<Element> &v1)
{
    ifstream fin;
    fin.open(fname.c_str());
    if (!fin.is_open())
    {
        cerr << fname << ": open file failed" << endl;
        exit(-1);
    }

    string line;
    getline(fin, line);

    char comma;
    while (true)
    {
        Element e;
        fin >> e.id >> comma >> e.mat >> comma >> e.typ >> comma >> e.sec
            >> comma >> e.p1 >> comma >> e.p2;
        if(!fin.good()) break;
        v1.push_back(e);
    }
    fin.close();
}

void saveExternalForce(const std::string &path, Particle* p)
{
    // create file
    string fname = path + "/particle_" + to_string(p->id())+ ".csv";
    fstream fout;
    fout.open(fname, ios::out | ios::app);

    // write header
    string header = "PID, Fx, Fy, Fz, Mx, My, Mz";
    fout << header << endl;

    // write data
    p->outputReactionForce(fout);

    // close file
    fout.close();
}

int main()
{
    // create system
    StructSystem system = StructSystem(1);
    system.setJobname("Star");

    // solve parameters
    double endTime = 100.0;
    double zeta = 1;
    double d = 6.0;
    double rtime = 50.0;

    // import structure parameters
    string fnode = system.workdir() + "/" + "nodes.csv";
    vector<Point> points;
    importNodes(fnode, points);

    string felems = system.workdir() + "/" + "elements.csv";
    vector<Element> elements;
    importElements(felems, elements);

    // create particles
    map<int, Particle*> particles;
    for (int i = 0; i < points.size(); i++)
    {
        Particle* p = new Particle(points[i].id, points[i].x, points[i].y, points[i].z);
        system.addParticle(p);
        particles[points[i].id] = p;
    }

    // create material
    LinearElastic mat = LinearElastic(1);
    mat.setE(1.0E6);
    mat.setDensity(20);

    // create section
    CustomSectionParas paras = {.A=0.1, .Sy=0, .Sz=0, .SHy=0, .SHz=0,
                          .CGy=0, .CGz=0.0, .Iyy=0.01, .Izz=0.01, .Iyz=10};
    CustomSection sect = CustomSection(1, paras);

    // create elements
    for (int i = 0; i < elements.size(); i++)
    {
        int id = elements[i].id;
        int p1 = elements[i].p1;
        int p2 = elements[i].p2;
        Link3DLD* e = new Link3DLD(id, particles[p1], particles[p2],
                           &mat, &sect);
        system.addElement(e);
    }

    // constraints
    DOF v1 = {.key={true, true, true, true, true, true},
              .val={0, 0, 0, 0, 0, 0}};
    system.addConstraint(8, v1);
    system.addConstraint(9, v1);
    system.addConstraint(10, v1);
    system.addConstraint(11, v1);
    system.addConstraint(12, v1);
    system.addConstraint(13, v1);
    system.info();

    string path = system.workdir() + "/" + system.jobname() + "/model";
    system.saveModel(path);
    StdArray6d d1 {};

    // setting parameters
    double h = system.autoTimeStep();
    system.setDampCoeff(zeta);
    int nStep = ceil(endTime / h);

    int interval = ceil(0.1 / h);
    cout << "##### start calculating ######" << endl;
    for (int i = 0; i <= nStep; i++)
    {
        d1[2] = -d * (double)i / nStep;
        if (i == 0)
        {
            system.setParticleDisplace(1, d1);
            system.solve(h, zeta, true);
            system.setInternalForce();
        }
        else
        {
            system.solve(h, zeta, false);
            system.clearParticleForce();
            system.setParticleDisplace(1, d1);
            system.setInternalForce();
        }
        // save results
        if (i % interval == 0)
        {
            string path = system.workdir() + "/" + system.jobname() + "/" +
                          to_string(i * h);
            system.saveParticleResult(path);
            system.saveElementResult(path);
            system.saveSupportReact(path);
            saveExternalForce(path, particles[1]);
            saveExternalForce(path, particles[2]);
        }
    }
    cout << "##### finished ######" << endl;

    system.releaseContainers();
    return 0;
}

​ Fig. 12 force-displace curve, (a) \(F-U_z\) at particle 1 (b) \(F-U_x\) at particle 2, (c) \(F-U_z\) at particle 2

More examples will be added in the near future, I believe you can get the basic concepts of openVFIFE.

Hope you can have a good experience in using openVFIFE.

I will write a detailed documentation for users and developers once I graduate. And I am looking forward to your advices and suggestions.

Best wishes!

  • Post title：Examples of openVFIFE
  • Post author：GinkGo
  • Create time：2021-06-13 00:41:02
  • Post link：https://ginkgoltd.github.io/2021/06/13/openVFIFE_examples/
  • Copyright Notice：All articles in this blog are licensed under BY-NC-SA unless stating additionally.