CFD Online Discussion Forums

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   Fluent UDF and Scheme Programming (https://www.cfd-online.com/Forums/fluent-udf/)
-   -   Understanding of contact detection UDF (https://www.cfd-online.com/Forums/fluent-udf/236691-understanding-contact-detection-udf.html)

Silence June 10, 2021 04:30

Understanding of contact detection UDF
 
Hello, everyone:
I am now using ‘contact detection’ function to keep a distance between two faces(one is stationary and the other is moving), some part of the UDF for contact detection is a little hard for me to understand since some macros cannot be found in UDF manual, so I am here asking for help to figure the whole things out.
UDF code and my personal understanding of it have been presented in Figure 1 and Figure 2.
Any advice will be appreciated, please help me understand those macros and share your point on the UDF code and my corresponding comments, thanks a lot!
Regards,
Silence

UDF macros which cannot be found in UDF manual:
1.What is ‘Objp’? (Part 1, line 6)
2. Loop(o, contacts) (Par2 line 33)
Loop the elements which pointer 'o' points to in contacts?
3.face = O_F(o) (Par2 line 35)
4.thread = O_F_THREAD(o) (Par2 line 36)
5.REAL_MIN (Par3 line 79)
6. SDOF_Get_Motion(dt, vel0, omega0, theta0) (Par4 line 113)
Get the motion property of the structure which thread dt points to, velocity(vel0), revolution speed(omega0) and degree produced by revolution(theta0), am I right?
7. SDOF_Overwrite_Motion(ndt, vel1, omega1, theta1) (Par5 line 144)
Overwrite the motion property of thread ‘ndt’, the new property of motion is: velocity_new = vel1, revolution speed_new = omega1 and degree_new = theta1, am I right?

Logic problem
1./* 'Fetch current thread ID'---UDF manual, */, what does 'current thread' mean? (part 2, line 29),
2.tid = THREAD_ID(DT_THREAD(dt)), 'dt' was not defined in previous statement, why it can be used here? (part 2, line 30)
3.nid = -1; what is nid? why nid = -1? (part 2, line 31)
4.n_faces++, Obtain the total number of faces on face ‘nid’? (part 2, line 56)
5. N3V_S(xc, /=, n_faces), Is this operation to obtain the centroid of whole face? if so, is it an axiom? (part, line 81)
6. ndt = THREAD_DT(thread), convert thread to dynamic thread, thread is the pointer of face ‘nid’, why do this operation? (part 4, line 110)
7./* Reflect velocity across the normal */
v1dotn0 = 2 * N3V_DOT(vel1, norm0);
N3V_S(norm0, *=, v1dotn0);
N3V_V(vel1, -=, norm0);

why not just N3V_V(vel1, -=, 2*vel1), why obtain the final vel1 through normal vector? (part 5, line 139-142)?

https://i.loli.net/2021/06/10/WX7BUNv8EkQheAP.png
Figure 1a

https://i.loli.net/2021/06/10/bmRT3lfEe5jh26t.png
Figure 1b

https://i.loli.net/2021/06/10/ax3ZS1nJsDrjBE8.png
Figure 1c

https://i.loli.net/2021/06/10/y4KHAUx21PEaX5h.png
Figure 1d

https://i.loli.net/2021/06/10/bz4U3hluafcn7xZ.png
Figure 1e


Code:

#include "udf.h"

DEFINE_CONTACT(contact_props, dt, contacts)
{
        /* Define variables */
        Objp* o;  /* What is Objp? */
        face_t face;
        Thread* thread; 
        Domain* domain = NULL;
        Dynamic_Thread* ndt = NULL;
        int tid, nid, n_faces;
        real v0dotn1, v1dotn0;
        real nc_mag, norm0_mag, norm1_mag;
        /* Define 3D vectors */
        real N3V_VEC(vel_rel);
        real N3V_VEC(nc), N3V_VEC(nctmp);
        real N3V_VEC(xc), N3V_VEC(xctmp);
        real N3V_VEC(vel0), N3V_VEC(omega0), N3V_VEC(theta0), N3V_VEC(norm0);
        real N3V_VEC(vel1), N3V_VEC(omega1), N3V_VEC(theta1), N3V_VEC(norm1);
        /* Check accessibility of cell values of variable in UDF */
        if (!Data_Valid_P())
        {
                return;
        }
        /* Define a common contact point / plane */
        N3V_S(nc, =, 0.0);
        N3V_S(xc, =, 0.0);

        /* 'Fetch current thread ID'---UDF manual, What does 'current thread' mean? */
        tid = THREAD_ID(DT_THREAD(dt));  /* 'dt' was not defined in previous statement, why it can be used here? */
        nid = -1;  /* What is nid? why nid = -1??? */
        n_faces = 0;
        loop(o, contacts)  /*  What does this macro mean?  */
        {
                face = O_F(o);  /* What does this macro mean? */
                thread = O_F_THREAD(o);  /* What does this macro mean? */

                /* Skip faces on current thread */
                if (THREAD_ID(thread) == tid)
                {
                        continue;
                }
                /* Note ID of the other thread for posterity */
                if (nid == -1)
                {
                        nid = THREAD_ID(thread);
                }
                /* Initialize to zero */
                N3V_S(nctmp, =, 0.0); 
                N3V_S(xctmp, =, 0.0); 
                F_AREA(nctmp, face, thread);
                F_CENTROID(xctmp, face, thread);
                /* Negative sum because wall normals point out of the fluid domain */
                N3V_V(nc, -=, nctmp);
                N3V_V(xc, +=, xctmp);
                n_faces++;  /* Obtain the total number of faces on face nid? */
        }

        /* Assign calculation task to node in case of some errors, such as division of zero */
# if RP_NODE
        {
                /* Reduce in parallel */
                nid = PRF_GIHIGH1(nid);  /* Return the maximum value of nid among all of the compute nodes,
                                                                    why do this operation? */
                n_faces = PRF_GISUM1(n_faces);
                PRF_GRSUM3(nc[0], nc[1], nc[2]);
                PRF_GRSUM3(xc[0], xc[1], xc[2]);
        }
# endif
        /* Propagate to host, sending data between host and node */
        node_to_host_int_2(nid, n_faces);
        node_to_host_real_3(nc[0], nc[1], nc[2]);
        node_to_host_real_3(xc[0], xc[1], xc[2]);
        /* Obtain parameters for nid face:
          (1) unit normal vector;
          (2) face gravity center. */
        if (n_faces > 0)
        {
                nc_mag = N3V_MAG(nc) + REAL_MIN;  /* What is the 'REAL_MIN'? */
                N3V_S(nc, /=, nc_mag);  /* Normalization for nc */
                N3V_S(xc, /=, n_faces);  /* Is this operation to obtain the centroid of whole face? if so, is it an axiom? */
        }
        else
        {
                return;
        }

#if RP_HOST
        Message
        (
                "\nContact:: tid: %d nid: %d n_faces: %d "
                "Point: (%f %f %f) Normal: (%f %f %f)",
                tid, nid, n_faces,
                xc[0], xc[1], xc[2],
                nc[0], nc[1], nc[2]
        );
#endif

        /* Fetch thread for opposite body */
        domain = THREAD_DOMAIN(DT_THREAD(dt));  /* What does this macro mean? */
        thread = Lookup_Thread(domain, nid);
        /* Check whether thread is null */
        if (NULLP(thread))
        {
                Message("\nWarning: No thread for nid: %d ", nid);
                return;
        }
        else
        {
                ndt = THREAD_DT(thread);  /* Convert thread to dynamic thread, why do this operation? */
        }
        /* Fetch body parameters */
        SDOF_Get_Motion(dt, vel0, omega0, theta0);
        /* Compute difference vectors and normalize */
        N3V_VV(norm0, =, xc, -, DT_CG(dt));  /* Centroid - center of gravity, why do this operation? */
        norm0_mag = N3V_MAG(norm0) + REAL_MIN;
        N3V_S(norm0, /=, norm0_mag); 
        if (NULLP(ndt))
        {
                /* Stationary body / wall. Use contact normal */
                N3V_V(norm1, =, nc);
                N3V_S(vel1, =, 0.0);
                N3V_V(vel_rel, =, vel0);  /* Compute relative velocity */
        }
        else
        {
                /* Fetch body parameters */
                SDOF_Get_Motion(ndt, vel1, omega1, theta1);
                /* Compute relative velocity */
                N3V_VV(vel_rel, =, vel0, -, vel1);
                N3V_VV(norm1, =, xc, -, DT_CG(ndt)); 
                norm1_mag = N3V_MAG(norm1) + REAL_MIN;
                N3V_S(norm1, /=, norm1_mag);


                /* Check if velocity needs to be reversed */
                if (N3V_DOT(vel_rel, nc) < 0.0)
                {
                        /* Reflect velocity across the normal, why not just N3V_V(vel1, -=, 2*vel1)? */
                        v1dotn0 = 2 * N3V_DOT(vel1, norm0);
                        N3V_S(norm0, *=, v1dotn0);
                        N3V_V(vel1, -=, norm0);
                        /* Override body velocity */
                        SDOF_Overwrite_Motion(ndt, vel1, omega1, theta1);
                }
        }

        /* Check if velocity needs to be reversed */
        if (N3V_DOT(vel_rel, nc) < 0.0)
        {
                /* Reflect velocity across the normal */
                v0dotn1 = 2 * N3V_DOT(vel0, norm1);
                N3V_S(norm1, *=, v0dotn1);
                N3V_V(vel0, -=, norm1);
                /* Override body velocity */
                SDOF_Overwrite_Motion(dt, vel0, omega0, theta0);
        }

#if RP_HOST
        Message
        (
                "\ncontact_props: Updated :: vel0 = (%f %f %f) vel1 = (%f %f %f)",
                vel0[0], vel0[1], vel0[2], vel1[0], vel1[1], vel1[2]
        );
#endif
}



All times are GMT -4. The time now is 18:16.