[Fenil1341]

Hello,

I am new to UDF and I need help me with a problem. Thanks in advance.

So I want to simulate an acoustic wave on a flat plate. but there is no equation that can be used so I have to use UDF. the equation are in the picture and with respect to those equations I want to give value of the 4 variable purturbation values and solve (u, v, P, and rho) on each mesh nodes for the solution. But I dont know how can I implement them on the whole mesh. the code for UDF I have written is only for the boundary condition.(code is below attached for the reference.)

also H_c_real, dH_c_real, H_c_imag, and dH_c_imag variables are heun function which varies with y coordinate. I have written only 3 for only the code to work but originallly and the array values of them are calculated for each mesh node with the help of MATLAB code.

********EQUATION**********
Screenshot 2024-07-15 181737.png



*********UDF CODE**********


#include "udf.h"
#include <math.h>
#include <complex.h>

double M = 0.5;
double omega = 5.118277;
double P = 101325;

static double H_c_real[] = {0.24509769, 0.245235439, 0.245373191};
static double H_c_imag[] = {-0.381717036, -0.381763624, -0.381810167};
static double dH_c_real[] = {1.073832121, 1.074074661, 1.074317076};
static double dH_c_imag[] = {-0.363362371, -0.363084281, -0.362806137};
static int array_size = sizeof(H_c_real) / sizeof(double); // Assuming all arrays have the same size


/******************** X - VELOCITY PROFILE ********************/



DEFINE_PROFILE(x_velocity, thread, i)
{
real x[ND_ND];
double y;
double theta;
double eta;
face_t f;
double alpha = 2;
double current_time;
double eta_prime;
current_time = CURRENT_TIME;
double x_vel = 1700;

begin_f_loop(f, thread)
{

F_CENTROID(x, f, thread);
y = x[1];
theta = -pow(M, 2) * pow(2 - omega, 2) + pow(2, 2);
eta = creal(_Cmulcc(_Cbuild(0.0, M * 2 * exp(-y)), _Cbuild(0.0, 1.0))) * cos(pow(theta, 0.5) * y);
eta_prime = -M * exp(-y) - sqrt(theta); /* Calculate the derivative of eta with respect to y */

double denominator = M * M * (omega - alpha + alpha * exp(-y));
double inner_expression = (alpha / denominator) - (eta_prime * exp(-y) / (denominator * denominator));

// int index = (int)(y / (1.0 / (array_size - 1))); // Assuming y ranges from 0 to 1
// if (index < 0) index = 0;
// if (index >= array_size) index = array_size - 1;

// _Dcomplex H_c = _Cbuild(H_c_real[index], H_c_imag[index]);
// _Dcomplex dH_c = _Cbuild(dH_c_real[index], dH_c_imag[index]);
_Dcomplex H_c = _Cbuild(H_c_real[1], H_c_imag[1]);
_Dcomplex dH_c = _Cbuild(dH_c_real[1], dH_c_imag[1]);

_Dcomplex complex_value = cexp(_Cbuild(0.0, -omega * current_time));
_Dcomplex result = _Cmulcc(_Cbuild(eta, 0.0), complex_value);
F_PROFILE(f, thread, i) = x_vel + creal(result) * (creal(_Cmulcc(H_c, _Cbuild(inner_expression, 0.0))) + creal(_Cmulcc(dH_c, _Cbuild((alpha * exp(-2 * y)) / (M * M * denominator * denominator * (alpha - omega)), 0.0))));
}
end_f_loop(f, thread)
}


/******************** Y - VELOCITY PROFILE ********************/


DEFINE_PROFILE(y_velocity, thread, i)
{

real x[ND_ND];
double y;
double theta;
double eta;
face_t f;
double alpha = 2;
double OMEGA = 5.118277;
double t;
double eta_prime;
t = CURRENT_TIME;
double y_vel = 0;

begin_f_loop(f, thread)
{

F_CENTROID(x, f, thread);
y = x[1];

theta = -pow(M, 2) * pow(2 - omega, 2) + pow(2, 2);
eta = creal(_Cmulcc(_Cbuild(0.0, M * 2 * exp(-y)), _Cbuild(0.0, 1.0))) * cos(pow(theta, 0.5) * y);
eta_prime = -M * exp(-y) - sqrt(theta);

double denominator = M * M * (omega - alpha + alpha * exp(-y));

_Dcomplex Hc = _Cbuild(H_c_real[1], H_c_imag[1]); ////Hc = Hc_re[cell_index];
_Dcomplex Hc_prime = _Cbuild(dH_c_real[1], dH_c_imag[1]); ////Hc_prime = Hc_prime_re[cell_index];

_Dcomplex complex_value = _Cmulcc(_Cmulcc(cexp(_Cbuild(eta, 0.0)), cexp(_Cbuild(0.0, alpha * x[0] - OMEGA * t))), _Cbuild((eta_prime * creal(Hc) - (alpha * exp(-y) / (alpha - OMEGA)) * creal(Hc_prime)) / (denominator), (eta_prime * cimag(Hc) - (alpha * exp(-y) / (alpha - OMEGA)) * cimag(Hc_prime)) / denominator));
F_PROFILE(f, thread, i) = y_vel + creal(complex_value);


//F_PROFILE(f, thread, i) = cexp(eta) * cexp(I * (alpha * x - OMEGA * t)) * (eta_prime * Hc - (alpha * exp(-y) / (alpha - OMEGA)) * Hc_prime) / (I * M * M * (OMEGA - alpha + alpha * exp(-y)));
}
end_f_loop(f, thread)
}



/******************** DENSITY PROFILE ********************/



DEFINE_PROPERTY(Density_profile, thread, i)
{
real x[ND_ND]; /* this will hold the position vector */
real y;
double theta;
double eta;
face_t f;
double alpha = 2;
double rho = 1.125;
double C = 1;
_Dcomplex rho_prime;
real t = CURRENT_TIME;

begin_f_loop(f, thread)
{
F_CENTROID(x, f, thread);
y = x[1];

theta = -pow(M, 2) * pow(2 - omega, 2) + pow(2, 2);
eta = creal(_Cmulcc(_Cbuild(0.0, M * 2 * exp(-y)), _Cbuild(0.0, 1.0))) * cos(pow(theta, 0.5) * y);

_Dcomplex Hc = _Cbuild(H_c_real[1], H_c_imag[1]);

rho_prime = _Cmulcc(_Cmulcc(_Cbuild(C, 0.0), cexp(_Cbuild(eta, 0.0))), _Cmulcc(cexp(_Cbuild(0.0, alpha * x[0] - omega * t)), Hc));
F_PROFILE(f, thread, i) = rho + creal(rho_prime);
}
end_f_loop(f, thread)
}




/******************** PRESSURE PROFILE ********************/




DEFINE_PROFILE(Pressure_profile, thread, i)
{


real x[ND_ND]; /* this will hold the position vector */
real y;
double theta;
double eta;
face_t f;
double alpha = 2;
double rho = 1.125;
double C = 1;
_Dcomplex rho_prime;
real t = CURRENT_TIME;
double Pressure = 101325;

begin_f_loop(f, thread)

{
F_CENTROID(x, f, thread);
y = x[1];

theta = -pow(M, 2) * pow(2 - omega, 2) + pow(2, 2);
eta = creal(_Cmulcc(_Cbuild(0.0, M * 2 * exp(-y)), _Cbuild(0.0, 1.0))) * cos(pow(theta, 0.5) * y);

_Dcomplex Hc = _Cbuild(H_c_real[1], H_c_imag[1]);

rho_prime = _Cmulcc(_Cmulcc(_Cbuild(C, 0.0), cexp(_Cbuild(eta, 0.0))), _Cmulcc(cexp(_Cbuild(0.0, alpha * x[0] - omega * t)), Hc));
F_PROFILE(f, thread, i) = Pressure + 1.4 * creal(rho_prime);
}

end_f_loop(f, thread)


}