bangerth/biharmonic

TODO: Automatically refine the mesh

Closed this issue · 3 comments

TODO: Automatically refine the mesh

From Jason:

		double c1 = jThinPlateBendingSpeedOfSound(freq, h, rho, E.r(), v);
		double c2 = sqrt(T.r() / (rho*h));  // speed of sound in stretched membrane (m/s)
                // T=tension (N/m), rho-density(kg/m^3), h=thickness (m)
		double c;
		if (c1 > c2)
			c = c1;
		else
			c = c2;
		double wavelength = c / freq;

#define MIN_DISCRETIZATION_rect_fea_d 15
// minimum # of discretization in x or y direction for this solver
#define DOF_MAX_rect_fea_d 1000
// # maximum of degrees of freedom for this solver

		int Nx = (int)(10.0 * Lx / wavelength);  // # of discretizations in x direction
		int Ny = (int)(10.0 * Ly / wavelength);  // # of discretizations in y direction
		if (Nx < MIN_DISCRETIZATION_rect_fea_d)
			Nx = MIN_DISCRETIZATION_rect_fea_d;
		if (Ny < MIN_DISCRETIZATION_rect_fea_d)
			Ny = MIN_DISCRETIZATION_rect_fea_d;

		if (Nx*Ny>DOF_MAX_rect_fea_d)// ok, we've got to reduce the size
		{
			double a = Nx*Ny;
			double s = sqrt(a / DOF_MAX_rect_fea_d);

			Nx = (int)(Nx / s);
			Ny = (int)(Ny / s);
			if (Nx < MIN_DISCRETIZATION_rect_fea_d)
			{
				Nx = MIN_DISCRETIZATION_rect_fea_d;
				if (Nx*Ny>DOF_MAX_rect_fea_d)
					Ny = (int)((double)DOF_MAX_rect_fea_d / (double)Nx);
			}
			if (Ny < MIN_DISCRETIZATION_rect_fea_d)
			{
				Ny = MIN_DISCRETIZATION_rect_fea_d;
				if (Nx*Ny>DOF_MAX_rect_fea_d)
					Nx = (int)((double)DOF_MAX_rect_fea_d / (double)Ny);
			}
		}



// computes the speed of sound of a shear or vending wave in a thin plate.
// send:
//		freq		frequency (Hz)
//		h			heigth or thickness of sheet (m)
//		rho			density of material (kg/m^3)
//		E			youngs modulus of material (Pa)
//		v			Poissons ratio of material (dimensionless)
//					(I think v=0 for a sheet of woven material or felt is correct)
// returns:
//		speed of sound in m/s
double jThinPlateBendingSpeedOfSound(double freq, double h, double rho, double E, double v)
{
	double I_w, w2,a;

	I_w = h*h*h / 12; // I/width
	// from http://www.engineeringtoolbox.com/area-moment-inertia-d_1328.html
	w2 = two_pi*freq;
	w2 *= w2;

	// from: http://scholar.lib.vt.edu/theses/available/etd-09172001-144121/unrestricted/Blanc_MsThesis.pdf
	//  speed of sound = (EI*w^2/rho*h*width)^(1/4)
	a = E*I_w*w2 / (rho*h);
	a = pow(a, 0.25);
	return a;
}

Fixed by f342204.

For the record, the formula to compute a in jThinPlateBendingSpeedOfSound is missing the factor 1/(1-nu^2) in the computation of E*I_w when describing the D factor in the equations.