diff --git a/.github/workflows/test_pvade.yaml b/.github/workflows/test_pvade.yaml index 897d46c..477ec3d 100644 --- a/.github/workflows/test_pvade.yaml +++ b/.github/workflows/test_pvade.yaml @@ -33,6 +33,8 @@ jobs: shell: bash -l {0} run: PYTHONPATH=. pytest -sv pvade/tests/ - name: test all inputs + env: + OMPI_MCA_regx: "naive" shell: bash -l {0} run: pytest -sv test_all_inputs.py diff --git a/input/duramat_case_study.yaml b/input/duramat_case_study.yaml index fb5939e..1f40a19 100644 --- a/input/duramat_case_study.yaml +++ b/input/duramat_case_study.yaml @@ -23,6 +23,13 @@ pv_array: elevation: 2.1 tracker_angle: 0.0 span_fixation_pts: [13.2] + torque_tube_separation: 0.2 # gap between panel and tube center + torque_tube_outer_radius: 0.1 # radius of the torque tube + torque_tube_inner_radius: 0.09 # radius of the torque tube + modules_per_span: 10 + fixed_location: 5 # fixed location of the panel along the span (from 0 (fixed at left) to modules_per_span (fixed at right)) + block_chord_div_by_panel_chord: 0.02 + solver: dt: 0.01 t_final: 20.0 @@ -58,4 +65,7 @@ structure: bc_list: [] motor_connection: true tube_connection: true - beta_relaxation: 0.5 \ No newline at end of file + beta_relaxation: 0.5 + elasticity_modulus_tube: 2.0e+11 + poissons_ratio_tube: 0.3 + rho_tube: 7800.0 \ No newline at end of file diff --git a/pvade/IO/input_schema.yaml b/pvade/IO/input_schema.yaml index 4aefbfa..18f1beb 100644 --- a/pvade/IO/input_schema.yaml +++ b/pvade/IO/input_schema.yaml @@ -136,7 +136,7 @@ properties: minimum: 0.0 maximum: 10.0 type: "number" - description: "The vertical distance between the center of the panel and the ground." + description: "The vertical distance between the center of the torque tube and the ground." units: "meter" stream_spacing: default: 7.0 @@ -198,6 +198,39 @@ properties: be unrolled such that tracking_angle_0 is applied to the panel at (x_min, y_min), tracking_angle_1 is applied to the panel at (x_min + x_spacing, y_min), etc., i.e., x-major direction. If argument is a single value, that tracking angle will be applied to the every panel." units: "degree" + torque_tube_separation: + default: 0.0 + minimum: 0.0 + maximum: 1.0 + type: "number" + description: "The distance between the bottom of the panel structure and the centerline of the torque tube, in meters." + torque_tube_outer_radius: + default: 0.0 + minimum: 0.0 + maximum: 1.0 + type: "number" + description: "The radius of the torque tube, in meters." + torque_tube_inner_radius: + default: 0.0 + minimum: 0.0 + maximum: 1.0 + type: "number" + description: "The radius of the torque tube, in meters." + modules_per_span: + default: 1 + minimum: 1 + type: "integer" + description: "The number of modules/panels per every panel_span distance. This can be thought of as the number of modules/panels per each table, counted in the spanwise direction only (i.e., 2-in-portait systems should report only the number in the spanwise direction, vs 2x that value)." + fixed_location: + default: 5 + type: "integer" + description: "If an integer between 0 and modules_per_span, indicates the fixed location" + block_chord_div_by_panel_chord: + default: 0.1 + minimum: 0.01 + type: "number" + description: "block length or width divided by panel chord" + solver: additionalProperties: false type: "object" @@ -648,6 +681,27 @@ properties: type: "number" description: "poissons ratio of the panel structure." units: "None" + elasticity_modulus_tube: + default: 200.0e+9 + minimum: 1.0e+2 + maximum: 500.0e+9 + type: "number" + description: "The effective Young's modulus of the torque tube." + units: "Pascal" + poissons_ratio_tube: + default: 0.3 + minimum: 0.1 + maximum: 0.9 + type: "number" + description: "poissons ratio of the torque tube." + units: "None" + rho_tube: + default: 1.0 + # minimum: 0.001 + # maximum: 1000.0 + type: "number" + description: "The density of the structure." + units: "kg/meter^3" body_force_x: default: 100 type: "number" diff --git a/pvade/fluid/FlowManager.py b/pvade/fluid/FlowManager.py index bd2c5c4..cdb97ce 100644 --- a/pvade/fluid/FlowManager.py +++ b/pvade/fluid/FlowManager.py @@ -165,6 +165,17 @@ def build_forms(self, domain, params): """ + # initialize total torque on each panel to zero + for panel_id in range( + int(params.pv_array.stream_rows * params.pv_array.span_rows) + ): + attr_name = f"total_torque_panel_{panel_id:.0f}" + setattr(self, attr_name, dolfinx.fem.Constant(domain.fluid.msh, 0.0)) + + for module_id in range(params.pv_array.modules_per_span+1): + attr_d_name = f"double_integral_total_torque_panel_{panel_id:.0f}_{module_id:.0f}" + setattr(self, attr_d_name, dolfinx.fem.Constant(domain.fluid.msh, 0.0)) + # Define fluid properties if self.ndim == 2: self.dpdx = dolfinx.fem.Constant(domain.fluid.msh, (0.0, 0.0)) @@ -242,6 +253,54 @@ def build_forms(self, domain, params): self.p_k = dolfinx.fem.Function(self.Q, name="pressure") self.p_k1 = dolfinx.fem.Function(self.Q) + # define function to save x, y, z coordinates + self.spatial_X_coords = dolfinx.fem.Function(self.Q, name="X_coords") + self.spatial_Y_coords = dolfinx.fem.Function(self.Q, name="Y_coords") + self.spatial_Z_coords = dolfinx.fem.Function(self.Q, name="Z_coords") + + class FillFunctionWithXCoords: + def __init__(self): + pass + + def __call__(self, x): + coords = np.zeros((1, x.shape[1]), dtype=PETSc.ScalarType) + + + coords[0] = x[0] + + return coords + + self.spatial_X_coords.interpolate(FillFunctionWithXCoords()) + + class FillFunctionWithYCoords: + def __init__(self): + pass + + def __call__(self, x): + coords = np.zeros((1, x.shape[1]), dtype=PETSc.ScalarType) + + + coords[0] = x[1] + + return coords + + self.spatial_Y_coords.interpolate(FillFunctionWithYCoords()) + + class FillFunctionWithZCoords: + def __init__(self): + pass + + def __call__(self, x): + coords = np.zeros((1, x.shape[1]), dtype=PETSc.ScalarType) + + + coords[0] = x[2] + + + return coords + + self.spatial_Z_coords.interpolate(FillFunctionWithZCoords()) + # initial conditions if params.fluid.initialize_with_inflow_bc: # self.inflow_profile = dolfinx.fem.Function(self.V) @@ -534,7 +593,38 @@ def _all_interior_surfaces(x): ds_fluid = ufl.Measure( "ds", domain=domain.fluid.msh, subdomain_data=domain.fluid.facet_tags ) - + + # def compute_total_torques_on_each_panel(self, domain, params): + + # ds_fluid = ufl.Measure( + # "ds", domain=domain.fluid.msh, subdomain_data=domain.fluid.facet_tags + # ) + + # for panel_id in range( + # int(params.pv_array.stream_rows * params.pv_array.span_rows) + # ): + # total_torque = 0 + + # for module_id in range(params.pv_array.modules_per_span): + + # total_torque += dolfinx.fem.assemble_scalar(dolfinx.fem.form(ufl.inner(self.spatial_X_coords, self.traction[2]) * ds_fluid( + # domain.domain_markers[f"panel_bottom_{panel_id:.0f}_{module_id:.0f}"]["idx"])) + # ) + # total_torque += dolfinx.fem.assemble_scalar(dolfinx.fem.form(ufl.inner(self.spatial_X_coords, self.traction[2]) * ds_fluid( + # domain.domain_markers[f"panel_top_{panel_id:.0f}_{module_id:.0f}"]["idx"])) + # ) + # total_torque += dolfinx.fem.assemble_scalar(dolfinx.fem.form(ufl.inner(self.spatial_Z_coords, self.traction[0]) * ds_fluid( + # domain.domain_markers[f"panel_top_{panel_id:.0f}_{module_id:.0f}"]["idx"])) + # ) + # total_torque += dolfinx.fem.assemble_scalar(dolfinx.fem.form(ufl.inner(self.spatial_Z_coords, self.traction[0]) * ds_fluid( + # domain.domain_markers[f"panel_bottom_{panel_id:.0f}_{module_id:.0f}"]["idx"])) + # ) + + # attr_name = f"total_torque_panel_{panel_id:.0f}" + + # setattr(self, attr_name, total_torque) + + # define integrated force forms self.integrated_force_x_form = [] self.integrated_force_y_form = [] self.integrated_force_z_form = [] @@ -548,79 +638,291 @@ def _all_interior_surfaces(x): # print(f"loc = {loc}, idx = {idx}, s = {s}") self.integrated_force_x_form.append(0) + self.integrated_force_y_form.append(0) + self.integrated_force_z_form.append(0) + self.integrated_force_x_form[-1] += self.traction[0] * ds_fluid( - domain.domain_markers[f"left_{panel_id:.0f}"]["idx"] - ) - self.integrated_force_x_form[-1] += self.traction[0] * ds_fluid( - domain.domain_markers[f"top_{panel_id:.0f}"]["idx"] - ) - self.integrated_force_x_form[-1] += self.traction[0] * ds_fluid( - domain.domain_markers[f"right_{panel_id:.0f}"]["idx"] - ) - self.integrated_force_x_form[-1] += self.traction[0] * ds_fluid( - domain.domain_markers[f"bottom_{panel_id:.0f}"]["idx"] - ) - if self.ndim == 3: - self.integrated_force_x_form[-1] += self.traction[0] * ds_fluid( - domain.domain_markers[f"front_{panel_id:.0f}"]["idx"] + domain.domain_markers[f"panel_left_{panel_id:.0f}"]["idx"] ) - self.integrated_force_x_form[-1] += self.traction[0] * ds_fluid( - domain.domain_markers[f"back_{panel_id:.0f}"]["idx"] + self.integrated_force_x_form[-1] += self.traction[0] * ds_fluid( + domain.domain_markers[f"panel_right_{panel_id:.0f}"]["idx"] ) - - self.integrated_force_x_form[-1] = dolfinx.fem.form( - self.integrated_force_x_form[-1] - ) - - self.integrated_force_y_form.append(0) - self.integrated_force_y_form[-1] += self.traction[1] * ds_fluid( - domain.domain_markers[f"left_{panel_id:.0f}"]["idx"] - ) - self.integrated_force_y_form[-1] += self.traction[1] * ds_fluid( - domain.domain_markers[f"top_{panel_id:.0f}"]["idx"] - ) self.integrated_force_y_form[-1] += self.traction[1] * ds_fluid( - domain.domain_markers[f"right_{panel_id:.0f}"]["idx"] - ) + domain.domain_markers[f"panel_left_{panel_id:.0f}"]["idx"] + ) self.integrated_force_y_form[-1] += self.traction[1] * ds_fluid( - domain.domain_markers[f"bottom_{panel_id:.0f}"]["idx"] - ) - if self.ndim == 3: - self.integrated_force_y_form[-1] += self.traction[1] * ds_fluid( - domain.domain_markers[f"front_{panel_id:.0f}"]["idx"] + domain.domain_markers[f"panel_right_{panel_id:.0f}"]["idx"] ) - self.integrated_force_y_form[-1] += self.traction[1] * ds_fluid( - domain.domain_markers[f"back_{panel_id:.0f}"]["idx"] + self.integrated_force_z_form[-1] += self.traction[2] * ds_fluid( + domain.domain_markers[f"panel_left_{panel_id:.0f}"]["idx"] ) - - self.integrated_force_y_form[-1] = dolfinx.fem.form( - self.integrated_force_y_form[-1] - ) - - self.integrated_force_z_form.append(0) + self.integrated_force_z_form[-1] += self.traction[2] * ds_fluid( + domain.domain_markers[f"panel_right_{panel_id:.0f}"]["idx"] + ) + if self.ndim == 3: - self.integrated_force_z_form[-1] += self.traction[2] * ds_fluid( - domain.domain_markers[f"left_{panel_id:.0f}"]["idx"] + self.integrated_force_x_form[-1] += self.traction[0] * ds_fluid( + domain.domain_markers[f"panel_front_{panel_id:.0f}"]["idx"] + ) + self.integrated_force_x_form[-1] += self.traction[0] * ds_fluid( + domain.domain_markers[f"panel_back_{panel_id:.0f}"]["idx"] + ) + self.integrated_force_y_form[-1] += self.traction[1] * ds_fluid( + domain.domain_markers[f"panel_front_{panel_id:.0f}"]["idx"] + ) + self.integrated_force_y_form[-1] += self.traction[1] * ds_fluid( + domain.domain_markers[f"panel_back_{panel_id:.0f}"]["idx"] + ) + self.integrated_force_z_form[-1] += self.traction[2] * ds_fluid( + domain.domain_markers[f"panel_front_{panel_id:.0f}"]["idx"] + ) + self.integrated_force_z_form[-1] += self.traction[2] * ds_fluid( + domain.domain_markers[f"panel_back_{panel_id:.0f}"]["idx"] + ) + for module_id in range(params.pv_array.modules_per_span): + + self.integrated_force_x_form[-1] += self.traction[0] * ds_fluid( + domain.domain_markers[f"panel_bottom_{panel_id:.0f}_{module_id:.0f}"]["idx"] ) - self.integrated_force_z_form[-1] += self.traction[2] * ds_fluid( - domain.domain_markers[f"top_{panel_id:.0f}"]["idx"] + self.integrated_force_x_form[-1] += self.traction[0] * ds_fluid( + domain.domain_markers[f"panel_bottom_{panel_id:.0f}_{module_id:.0f}"]["idx"] ) - self.integrated_force_z_form[-1] += self.traction[2] * ds_fluid( - domain.domain_markers[f"right_{panel_id:.0f}"]["idx"] + self.integrated_force_y_form[-1] += self.traction[1] * ds_fluid( + domain.domain_markers[f"panel_bottom_{panel_id:.0f}_{module_id:.0f}"]["idx"] ) - self.integrated_force_z_form[-1] += self.traction[2] * ds_fluid( - domain.domain_markers[f"bottom_{panel_id:.0f}"]["idx"] + self.integrated_force_y_form[-1] += self.traction[1] * ds_fluid( + domain.domain_markers[f"panel_bottom_{panel_id:.0f}_{module_id:.0f}"]["idx"] ) self.integrated_force_z_form[-1] += self.traction[2] * ds_fluid( - domain.domain_markers[f"front_{panel_id:.0f}"]["idx"] + domain.domain_markers[f"panel_bottom_{panel_id:.0f}_{module_id:.0f}"]["idx"] ) self.integrated_force_z_form[-1] += self.traction[2] * ds_fluid( - domain.domain_markers[f"back_{panel_id:.0f}"]["idx"] + domain.domain_markers[f"panel_bottom_{panel_id:.0f}_{module_id:.0f}"]["idx"] ) + - self.integrated_force_z_form[-1] = dolfinx.fem.form( - self.integrated_force_z_form[-1] - ) + self.integrated_force_x_form[-1] = dolfinx.fem.form( + self.integrated_force_x_form[-1] + ) + self.integrated_force_y_form[-1] = dolfinx.fem.form( + self.integrated_force_y_form[-1] + ) + self.integrated_force_z_form[-1] = dolfinx.fem.form( + self.integrated_force_z_form[-1] + ) + + def compute_double_integral_panel_torques(self, domain, params): + + ds_fluid = ufl.Measure( + "ds", domain=domain.fluid.msh, subdomain_data=domain.fluid.facet_tags + ) + + for panel_id in range( + int(params.pv_array.stream_rows * params.pv_array.span_rows) + ): + whole_top_surface_facet_index = [] # facet index of the whole panel top surface in a row + whole_bot_surface_facet_index = [] # facet index of the whole panel bottom surface in a row + + for module_id in range(params.pv_array.modules_per_span): + whole_top_surface_facet_index += domain.fluid.facet_tags.find(domain.domain_markers[f"panel_top_{panel_id:.0f}_{module_id:.0f}"]["idx"]).tolist() + whole_bot_surface_facet_index += domain.fluid.facet_tags.find(domain.domain_markers[f"panel_bottom_{panel_id:.0f}_{module_id:.0f}"]["idx"]).tolist() + + whole_top_surf_submesh, entity_map, vertex_map, geom_map = dolfinx.mesh.create_submesh(domain.fluid.msh, 2, np.array(whole_top_surface_facet_index, dtype=np.int32)) + whole_bot_surf_submesh, entity_map, vertex_map, geom_map = dolfinx.mesh.create_submesh(domain.fluid.msh, 2, np.array(whole_bot_surface_facet_index, dtype=np.int32)) + + spatial_X_coords_top = dolfinx.fem.Function(dolfinx.fem.FunctionSpace(whole_top_surf_submesh, ('CG', 1)), name="X_coords_top") + spatial_X_coords_bot = dolfinx.fem.Function(dolfinx.fem.FunctionSpace(whole_bot_surf_submesh, ('CG', 1)), name="X_coords_bot") + spatial_Z_coords_top = dolfinx.fem.Function(dolfinx.fem.FunctionSpace(whole_top_surf_submesh, ('CG', 1)), name="Z_coords_top") + spatial_Z_coords_bot = dolfinx.fem.Function(dolfinx.fem.FunctionSpace(whole_bot_surf_submesh, ('CG', 1)), name="Z_coords_bot") + spatial_X_coords_top.interpolate(self.spatial_X_coords) + spatial_X_coords_bot.interpolate(self.spatial_X_coords) + spatial_Z_coords_top.interpolate(self.spatial_Z_coords) + spatial_Z_coords_bot.interpolate(self.spatial_Z_coords) + """ + the torque should be relative to central of tube, traction x times z (z=0 at center tube) + traction z times x, (x=0 at center tube) + for z, how we should put the tube center at zero? It is not minus elevation + """ + + whole_top_submesh_function_space = dolfinx.fem.FunctionSpace(whole_top_surf_submesh, ('CG', 1)) + whole_bot_submesh_function_space = dolfinx.fem.FunctionSpace(whole_bot_surf_submesh, ('CG', 1)) + + traction_z_array_top = dolfinx.fem.Function(whole_top_submesh_function_space, name="traction_z_top") + traction_x_array_top = dolfinx.fem.Function(whole_top_submesh_function_space, name="traction_x_top") + traction_z_array_bot = dolfinx.fem.Function(whole_bot_submesh_function_space, name="traction_z_bot") + traction_x_array_bot = dolfinx.fem.Function(whole_bot_submesh_function_space, name="traction_x_bot") + + traction_z_fluid = dolfinx.fem.Function(self.Q, name="traction_z_fluid") + traction_x_fluid = dolfinx.fem.Function(self.Q, name="traction_x_fluid") + + # traction is expression on whole fluid mesh, we need to + # project it to a function defined on whole fluid mesh, then + # interpolate it to the top and bot submesh function space. + + u_inter_fluid = ufl.TrialFunction(self.Q) + v_inter_fluid = ufl.TestFunction(self.Q) + + a_inter_fluid = ufl.inner(u_inter_fluid, v_inter_fluid) * ufl.ds + L_inter_fluid_z = ufl.inner(self.traction[2], v_inter_fluid) * ufl.ds + L_inter_fluid_x = ufl.inner(self.traction[0], v_inter_fluid) * ufl.ds + + # Assemble and solve + A_inter_fluid = dolfinx.fem.petsc.assemble_matrix(dolfinx.fem.form(a_inter_fluid)) + A_inter_fluid.assemble() + + b_inter_fluid_z = dolfinx.fem.petsc.assemble_vector(dolfinx.fem.form(L_inter_fluid_z)) + b_inter_fluid_x = dolfinx.fem.petsc.assemble_vector(dolfinx.fem.form(L_inter_fluid_x)) + + solver = PETSc.KSP().create(traction_z_array_top.function_space.mesh.comm) + solver.setOperators(A_inter_fluid) + solver.setType("cg") + solver.getPC().setType("jacobi") + solver.setFromOptions() + solver.solve(b_inter_fluid_z, traction_z_fluid.vector) + traction_z_fluid.x.scatter_forward() + + solver = PETSc.KSP().create(traction_x_array_top.function_space.mesh.comm) + solver.setOperators(A_inter_fluid) + solver.setType("cg") + solver.getPC().setType("jacobi") + solver.setFromOptions() + solver.solve(b_inter_fluid_x, traction_x_fluid.vector) + traction_x_fluid.x.scatter_forward() + + traction_z_array_top.interpolate(traction_z_fluid) + traction_x_array_top.interpolate(traction_x_fluid) + traction_z_array_bot.interpolate(traction_z_fluid) + traction_x_array_bot.interpolate(traction_x_fluid) + + # traction_z_array_top.x.array[:] = 10*np.sin(params.pv_array.tracker_angle[0]*np.pi/180) + # traction_z_array_bot.x.array[:] = 10*np.sin(params.pv_array.tracker_angle[0]*np.pi/180) + # traction_x_array_top.x.array[:] = 10*np.cos(params.pv_array.tracker_angle[0]*np.pi/180) + # traction_x_array_bot.x.array[:] = 10*np.cos(params.pv_array.tracker_angle[0]*np.pi/180) + torque_function_on_panels_top = dolfinx.fem.Function(whole_top_submesh_function_space) + torque_function_on_panels_top.x.array[:] = ( + (spatial_X_coords_top.x.array[:] - np.linspace( + 0.0, + params.pv_array.stream_spacing * (params.pv_array.stream_rows - 1), + params.pv_array.stream_rows, + )[panel_id]) * traction_z_array_top.x.array[:] + - (spatial_Z_coords_top.x.array[:] - params.pv_array.elevation) * traction_x_array_top.x.array[:] + ) + + torque_function_on_panels_bot = dolfinx.fem.Function(whole_bot_submesh_function_space) + torque_function_on_panels_bot.x.array[:] = ( + (spatial_X_coords_bot.x.array[:] - np.linspace( + 0.0, + params.pv_array.stream_spacing * (params.pv_array.stream_rows - 1), + params.pv_array.stream_rows, + )[panel_id]) * traction_z_array_bot.x.array[:] + - (spatial_Z_coords_bot.x.array[:] - params.pv_array.elevation) * traction_x_array_bot.x.array[:] + ) + dx_top_whole = ufl.Measure("dx", domain=whole_top_surf_submesh) + dx_bot_whole = ufl.Measure("dx", domain=whole_bot_surf_submesh) + + total_torque_on_panel_array = ( + dolfinx.fem.assemble_scalar(dolfinx.fem.form(torque_function_on_panels_top * dx_top_whole)) + + dolfinx.fem.assemble_scalar(dolfinx.fem.form(torque_function_on_panels_bot * dx_bot_whole)) + ) + + attr_name = f"total_torque_panel_{panel_id:.0f}" + total_torque_constant = getattr(self, attr_name) + + total_torque_constant.value = total_torque_on_panel_array + + coords_top = whole_top_submesh_function_space.tabulate_dof_coordinates() + coords_bot = whole_bot_submesh_function_space.tabulate_dof_coordinates() + + top_integrated_func_along_y = np.zeros_like(torque_function_on_panels_top.x.array[:]) + bot_integrated_func_along_y = np.zeros_like(torque_function_on_panels_bot.x.array[:]) + + # in PVade, the front has the minimum y coordinates, if we integral from large y to small y, it is negative, + # to keep it consistant as the theoretical model, we need to flip the sign of the integral + + # Sort by y first, then z and z (group by y-levels) + sort_idx_top = np.lexsort((coords_top[:, 0], -coords_top[:, 1])) # y is sorted, from large to small + sort_idx_bot = np.lexsort((coords_bot[:, 0], -coords_bot[:, 1])) # y is sorted, from large to small + + top_coords_sorted = coords_top[sort_idx_top] #sorted from largest y to smallest y + bot_coords_sorted = coords_bot[sort_idx_bot] #sorted from largest y to smallest y + + ### integral of top surface: + # Unique x-values (up to numerical tolerance) + ys_top = np.unique(np.round(top_coords_sorted[:, 1], decimals=6)) + + # get f(x) = integral of torque on top surface from left end to x (x is all unique x coordinates on top surface) + indicator_top = 0 + for y_value in ys_top: + # Mask for points at this y-level + y_mask_top = np.isclose(top_coords_sorted[:, 1], y_value, atol=1e-6) + + facet_centers_top = dolfinx.mesh.compute_midpoints(whole_top_surf_submesh, 2, np.array(np.arange(whole_top_surf_submesh.topology.index_map(2).size_local), dtype=np.int32)) + # Identify cells where the midpoint y-coordinate is higher + marked_facet_top = np.where(facet_centers_top[:, 1] >= y_value)[0] + + # Create a submesh from those cells + submesh_top_surface_part, entity_map, vertex_map, geom_map = dolfinx.mesh.create_submesh(whole_top_surf_submesh, 2, marked_facet_top.astype(np.int32)) + top_sub_function_space = dolfinx.fem.FunctionSpace(submesh_top_surface_part, ('CG', 1)) + top_sub_func = dolfinx.fem.Function(top_sub_function_space) + top_sub_func.interpolate(torque_function_on_panels_top) + + dx_top = ufl.Measure("dx", domain=submesh_top_surface_part) + integrated_top = dolfinx.fem.assemble_scalar(dolfinx.fem.form(top_sub_func * dx_top)) + top_integrated_func_along_y[sort_idx_top[y_mask_top]] = integrated_top + indicator_top += 1 + ### integral of bot surface: + # Unique y-values (up to numerical tolerance) + ys_bot = np.unique(np.round(bot_coords_sorted[:, 1], decimals=6)) + indicator_bot = 0 + for y_value in ys_bot: + + # Mask for points at this y-level + y_mask_bot = np.isclose(bot_coords_sorted[:, 1], y_value, atol=1e-6) + facet_centers_bot = dolfinx.mesh.compute_midpoints(whole_bot_surf_submesh, 2, np.array(np.arange(whole_bot_surf_submesh.topology.index_map(2).size_local), dtype=np.int32)) + # Identify cells where the midpoint y-coordinate is less than 1 + marked_facet_bot = np.where(facet_centers_bot[:, 1] >= y_value)[0] + # Create a submesh from those cells + submesh_bot_surface_part, entity_map, vertex_map, geom_map = dolfinx.mesh.create_submesh(whole_bot_surf_submesh, 2, marked_facet_bot.astype(np.int32)) + bot_sub_function_space = dolfinx.fem.FunctionSpace(submesh_bot_surface_part, ('CG', 1)) + bot_sub_func = dolfinx.fem.Function(bot_sub_function_space) + bot_sub_func.interpolate(torque_function_on_panels_bot) + + dx_bot = ufl.Measure("dx", domain=submesh_bot_surface_part) + integrated_bot = dolfinx.fem.assemble_scalar(dolfinx.fem.form(bot_sub_func * dx_bot)) + bot_integrated_func_along_y[sort_idx_bot[y_mask_bot]] = integrated_bot + + indicator_bot += 1 + single_integral_function_top = dolfinx.fem.Function(whole_top_submesh_function_space, name='single_integral_function_top') + single_integral_function_top.x.array[:] = top_integrated_func_along_y + + single_integral_function_bot = dolfinx.fem.Function(whole_bot_submesh_function_space, name='single_integral_function_bot') + single_integral_function_bot.x.array[:] = bot_integrated_func_along_y + + torque_double_integral = 0.0 + + + attr_name = f"double_integral_total_torque_panel_{panel_id:.0f}_0" + # setattr(self, attr_name, torque_double_integral) + + double_integral_total_torque_panel_constant = getattr(self, attr_name) + double_integral_total_torque_panel_constant.value = torque_double_integral + + single_integral_function_top_fluid = dolfinx.fem.Function(self.Q) + single_integral_function_top_fluid.interpolate(single_integral_function_top) + single_integral_function_bot_fluid = dolfinx.fem.Function(self.Q) + single_integral_function_bot_fluid.interpolate(single_integral_function_bot) + + for connector_id in range(params.pv_array.modules_per_span): + torque_double_integral += dolfinx.fem.assemble_scalar(dolfinx.fem.form(single_integral_function_top_fluid * ds_fluid( + domain.domain_markers[f"panel_top_{panel_id:.0f}_{params.pv_array.modules_per_span-1-connector_id:.0f}"]["idx"])))/params.pv_array.panel_chord + torque_double_integral += dolfinx.fem.assemble_scalar(dolfinx.fem.form(single_integral_function_bot_fluid * ds_fluid( + domain.domain_markers[f"panel_bottom_{panel_id:.0f}_{params.pv_array.modules_per_span-1-connector_id:.0f}"]["idx"])))/params.pv_array.panel_chord + attr_name = f"double_integral_total_torque_panel_{panel_id:.0f}_{connector_id+1:.0f}" + + double_integral_total_torque_panel_constant = getattr(self, attr_name) + double_integral_total_torque_panel_constant.value = torque_double_integral + + # double_integral_total_torque_panel_0_0 is the double integral of first panel array at the most back connector (maximum y), + # double_integral_total_torque_panel_0_10 is the double integral of first panel array at the most front connector (smallest y) def _assemble_system(self, params): """Pre-assemble all LHS matrices and RHS vectors @@ -750,9 +1052,10 @@ def solve(self, domain, params, current_time): self.mesh_vel.vector.array[:] = self.mesh_vel.vector.array / params.solver.dt self.mesh_vel.x.scatter_forward() + # self.compute_double_integral_panel_torques(domain, params) + # Calculate the tentative velocity self._solver_step_1(params) - # Calculate the change in pressure to enforce incompressibility self._solver_step_2(params) @@ -770,10 +1073,13 @@ def solve(self, domain, params, current_time): self.compute_lift_and_drag(params, current_time) + # self.compute_panel_torques(domain, params) + self.compute_double_integral_panel_torques(domain, params) + # Compute the pressure drop between the inlet and outlet if params.pv_array.stream_rows > 0: self.compute_pressure_drop_between_points(domain, params) - + # Update new -> old variables self.u_k2.x.array[:] = self.u_k1.x.array self.u_k1.x.array[:] = self.u_k.x.array diff --git a/pvade/geometry/MeshManager.py b/pvade/geometry/MeshManager.py index 13f1764..f35a7aa 100644 --- a/pvade/geometry/MeshManager.py +++ b/pvade/geometry/MeshManager.py @@ -207,6 +207,8 @@ def build(self, params): gdim=self.ndim, ) + + self.msh.topology.create_connectivity(self.ndim, self.ndim - 1) # Specify names for the mesh elements @@ -214,6 +216,11 @@ def build(self, params): self.cell_tags.name = "cell_tags" self.facet_tags.name = "facet_tags" + with dolfinx.io.XDMFFile(self.comm, f"parent_mesh_only.xdmf", "w") as fp: + fp.write_mesh(self.msh) + fp.write_meshtags(self.cell_tags) + fp.write_meshtags(self.facet_tags) + # if ( # params.general.geometry_module == "panels3d" # or params.general.geometry_module == "heliostats3d" @@ -276,10 +283,19 @@ def _create_submeshes_from_parent(self, params): print(f"Creating {sub_domain_name} submesh") # Get the idx associated with either "fluid" or "structure" - marker_id = self.domain_markers[sub_domain_name]["idx"] - # Find all cells where cell tag = marker_id - submesh_cells = self.cell_tags.find(marker_id) + # if structure includes modules and connectors + if sub_domain_name == "structure" and "structure" not in self.domain_markers: + marker_id = self.domain_markers["modules"]["idx"] + # Find all cells where cell tag = marker_id + submesh_cells_modules = self.cell_tags.find(marker_id) + marker_id = self.domain_markers["connectors"]["idx"] + submesh_cells = np.hstack((self.cell_tags.find(marker_id), submesh_cells_modules)) + + else: + marker_id = self.domain_markers[sub_domain_name]["idx"] + # Find all cells where cell tag = marker_id + submesh_cells = self.cell_tags.find(marker_id) # Use those found cells to construct a new mesh submesh, entity_map, vertex_map, geom_map = dolfinx.mesh.create_submesh( @@ -309,7 +325,7 @@ def _transfer_mesh_tags_to_submeshes(self, params): num_facets = f_map.size_local + f_map.num_ghosts all_values = np.zeros(num_facets, dtype=np.int32) - # Assign non-zero facet tags using the facet tag indices + # Assign non-zero facet tags using the facet tag indices, save the facet marker to all_values all_values[self.facet_tags.indices] = self.facet_tags.values cell_to_facet = self.msh.topology.connectivity(self.ndim, facet_dim) @@ -345,14 +361,32 @@ def _transfer_mesh_tags_to_submeshes(self, params): for child, parent in zip(child_facets, parent_facets): sub_values[child] = all_values[parent] + # sub_cell_map = sub_domain.msh.topology.index_map(self.ndim) + f_map_cell = self.msh.topology.index_map(self.ndim) + + # Get the total number of cells in the parent mesh + num_cells = f_map_cell.size_local + f_map_cell.num_ghosts + all_cell_values = np.zeros(num_cells, dtype=np.int32) + all_cell_values[self.cell_tags.indices] = self.cell_tags.values + + sub_cell_map = sub_domain.msh.topology.index_map(self.ndim) sub_num_cells = sub_cell_map.size_local + sub_cell_map.num_ghosts + sub_cell_values = np.empty(sub_num_cells, dtype=np.int32) + + + for k, entity in enumerate(sub_domain.entity_map): + sub_cell_values[k] = all_cell_values[entity] + + + # sub_num_cells = sub_cell_map.size_local + sub_cell_map.num_ghosts + sub_domain.cell_tags = dolfinx.mesh.meshtags( sub_domain.msh, sub_domain.msh.topology.dim, np.arange(sub_num_cells, dtype=np.int32), - np.ones(sub_num_cells, dtype=np.int32), + sub_cell_values, ) sub_domain.cell_tags.name = "cell_tags" @@ -713,7 +747,7 @@ def test_mesh_functionspace(self): print(f"Rank {self.rank} owns {num_nodes_owned_by_proc} nodes\n{coords}") - def test_submesh_transfer(self, params): + def test_submesh_transfer(self, params, domain, elasticity): P2 = ufl.VectorElement("Lagrange", self.msh.ufl_cell(), 2) # P2 = ufl.FiniteElement("Lagrange", self.fluid.msh.ufl_cell(), 1) diff --git a/pvade/geometry/panels3d/DomainCreation.py b/pvade/geometry/panels3d/DomainCreation.py index 5ead1e4..fdd46af 100644 --- a/pvade/geometry/panels3d/DomainCreation.py +++ b/pvade/geometry/panels3d/DomainCreation.py @@ -149,11 +149,32 @@ def Rz(theta): domain_tag_list = [] domain_tag_list.append(domain_tag) + # panel_tag_list includes all structure, panel+connector panel_tag_list = [] panel_ct = 0 + + # only include panels + structure_panel_only_list = [] + # only include connectors + structure_connector_only_list = [] + + module_distances = np.linspace( + -half_span, half_span, params.pv_array.modules_per_span + 1 + ) - prev_surf_tag = [] + transformed_com = {} # all panels + + if ( + params.pv_array.torque_tube_separation > 0.0 + and params.pv_array.torque_tube_outer_radius > 0.0 + ): + modeling_torque_tube = True + else: + modeling_torque_tube = False + vol_tags_modules = [] + vol_tags_connectors = [] + # start to add panel for panel_id_y, yy in enumerate(y_centers): for panel_id_x, xx in enumerate(x_centers): if isinstance(params.pv_array.tracker_angle, list): @@ -169,186 +190,465 @@ def Rz(theta): else: tracker_angle_rad = np.radians(params.pv_array.tracker_angle) - # Create an 0-tracking-degree panel centered at (x, y, z) = (0, 0, 0) - panel_id = self.gmsh_model.occ.addBox( - -half_chord, - -half_span, - -half_thickness, - params.pv_array.panel_chord, - params.pv_array.panel_span, - params.pv_array.panel_thickness, - ) + this_panel_tag_list = [] # for each row + this_panel_transformed_com = {} # for each row + + numpy_pt_list = [] # for each row + embedded_lines_tag_list = [] # for each row - panel_tag = (self.ndim, panel_id) - panel_tag_list.append(panel_tag) + for module_id in range(params.pv_array.modules_per_span): - numpy_pt_list = [] - embedded_lines_tag_list = [] + module_span = ( + module_distances[module_id + 1] - module_distances[module_id] + ) - # Add a bisecting line to the bottom of the panel in the spanwise direction - pt_1 = self.gmsh_model.occ.addPoint(0, -half_span, -half_thickness) - pt_2 = self.gmsh_model.occ.addPoint(0, half_span, -half_thickness) + if modeling_torque_tube: + # Create an 0-tracking-degree panel centered at (x, y, z) = (0, 0, 0) + this_module = self.gmsh_model.occ.addBox( + -half_chord, + module_distances[module_id], + params.pv_array.torque_tube_separation, + params.pv_array.panel_chord, + module_span, + params.pv_array.panel_thickness, + ) + + panel_tag_list.append((self.ndim, this_module)) # all panels + this_panel_tag_list.append((self.ndim, this_module)) # this panel array/row + + structure_panel_only_list.append((self.ndim, this_module)) + + this_standoff = self.gmsh_model.occ.addBox( + -params.pv_array.block_chord_div_by_panel_chord * half_chord, + module_distances[module_id], + 0.0, + 2.0 * params.pv_array.block_chord_div_by_panel_chord * half_chord, + params.pv_array.block_chord_div_by_panel_chord * half_chord, + params.pv_array.torque_tube_separation, + ) + panel_tag_list.append((self.ndim, this_standoff)) + this_panel_tag_list.append((self.ndim, this_standoff)) - numpy_pt_list.append( - [0, -half_span, -half_thickness, 0, half_span, -half_thickness] - ) + structure_connector_only_list.append((self.ndim, this_standoff)) - torque_tube_id = self.gmsh_model.occ.addLine(pt_1, pt_2) - torque_tube_tag = (1, torque_tube_id) - embedded_lines_tag_list.append(torque_tube_tag) + # Add a bisecting line to the bottom of the connector in the spanwise direction + pt_1 = self.gmsh_model.occ.addPoint(0, module_distances[module_id], 0.0) + pt_2 = self.gmsh_model.occ.addPoint(0, module_distances[module_id] + params.pv_array.block_chord_div_by_panel_chord * half_chord, 0.0) + numpy_pt_list.append([0, module_distances[module_id], 0.0, 0, module_distances[module_id] + params.pv_array.block_chord_div_by_panel_chord * half_chord, 0.0]) # for this row + torque_tube_id = self.gmsh_model.occ.addLine(pt_1, pt_2) + torque_tube_tag = (1, torque_tube_id) + embedded_lines_tag_list.append(torque_tube_tag) # for this panel row - # Add lines in the streamwise direction to mimic sections of panel held rigid by motor - if params.pv_array.span_fixation_pts is not None: - if not isinstance(params.pv_array.span_fixation_pts, list): - num_fixation_pts = int( - np.floor( - params.pv_array.panel_span - / params.pv_array.span_fixation_pts - ) + else: + this_module = self.gmsh_model.occ.addBox( + -half_chord, + module_distances[module_id], + 0.0, + params.pv_array.panel_chord, + module_span, + params.pv_array.panel_thickness, ) + panel_tag_list.append((self.ndim, this_module)) + this_panel_tag_list.append((self.ndim, this_module)) + + structure_panel_only_list.append((self.ndim, this_module)) + + pt_1 = self.gmsh_model.occ.addPoint(0, module_distances[module_id], 0.0) + pt_2 = self.gmsh_model.occ.addPoint(0, module_distances[module_id+1], 0.0) + numpy_pt_list.append([0, module_distances[module_id], 0.0, 0, module_distances[module_id+1], 0.0]) + torque_tube_id = self.gmsh_model.occ.addLine(pt_1, pt_2) + torque_tube_tag = (1, torque_tube_id) + embedded_lines_tag_list.append(torque_tube_tag) # for this panel row + + if modeling_torque_tube: # add the last connector + last_standoff = self.gmsh_model.occ.addBox( + -params.pv_array.block_chord_div_by_panel_chord * half_chord, + module_distances[module_id+1] - params.pv_array.block_chord_div_by_panel_chord * half_chord, + 0.0, + 2.0 * params.pv_array.block_chord_div_by_panel_chord * half_chord, + params.pv_array.block_chord_div_by_panel_chord * half_chord, + params.pv_array.torque_tube_separation, + ) + this_panel_tag_list.append((self.ndim, last_standoff)) + panel_tag_list.append((self.ndim, last_standoff)) + + structure_connector_only_list.append((self.ndim, last_standoff)) + + pt_1 = self.gmsh_model.occ.addPoint(0, module_distances[module_id+1] - params.pv_array.block_chord_div_by_panel_chord * half_chord, 0.0) + pt_2 = self.gmsh_model.occ.addPoint(0, module_distances[module_id+1], 0.0) + numpy_pt_list.append([0, module_distances[module_id+1] - params.pv_array.block_chord_div_by_panel_chord * half_chord, 0.0, 0, module_distances[module_id+1], 0.0]) + torque_tube_id = self.gmsh_model.occ.addLine(pt_1, pt_2) + torque_tube_tag = (1, torque_tube_id) + embedded_lines_tag_list.append(torque_tube_tag) + + """ this is not used as we are using block surface to simulate the fixed bc applied by motor """ + # # Add lines in the streamwise direction to mimic sections of panel held rigid by motor + # if params.pv_array.span_fixation_pts is not None: + # if not isinstance(params.pv_array.span_fixation_pts, list): + # num_fixation_pts = int( + # np.floor( + # params.pv_array.panel_span + # / params.pv_array.span_fixation_pts + # ) + # ) + + # fixation_pts_list = [] + + # for k in range(1, num_fixation_pts + 1): + # next_pt = k * params.pv_array.span_fixation_pts + + # eps = 1e-9 + + # if ( + # next_pt > eps + # and next_pt < params.pv_array.panel_span - eps + # ): + # fixation_pts_list.append(next_pt) + + # else: + # fixation_pts_list = params.pv_array.span_fixation_pts + + # for fp in fixation_pts_list: + # # FIXME: don't add the fixation points into the numpy tagging for now + # if modeling_torque_tube: + # numpy_pt_list.append( + # [ + # -params.pv_array.torque_tube_radius, + # -half_span + fp, + # 0.0, + # params.pv_array.torque_tube_radius, + # -half_span + fp, + # 0.0, + # ] + # ) + + # else: + # numpy_pt_list.append( + # [ + # -half_chord, + # -half_span + fp, + # 0.0, + # half_chord, + # -half_span + fp, + # 0.0, + # ] + # ) + + # # If the separation line already exists at a module division, + # # there's no need to redraw it, but otherwise, add the gmsh points + # # to the model for embedding + # if not np.any(np.isclose(-half_span + fp, module_distances)): + # print( + # f"Embedding a line for a no-deformation boundary condition." + # ) + + # if modeling_torque_tube: + # pt_1 = self.gmsh_model.occ.addPoint( + # -params.pv_array.torque_tube_radius, + # -half_span + fp, + # 0.0, + # ) + # pt_2 = self.gmsh_model.occ.addPoint( + # params.pv_array.torque_tube_radius, + # -half_span + fp, + # 0.0, + # ) + + # else: + # pt_1 = self.gmsh_model.occ.addPoint( + # -half_chord, -half_span + fp, 0.0 + # ) + # pt_2 = self.gmsh_model.occ.addPoint( + # half_chord, -half_span + fp, 0.0 + # ) + + # fixation_line_id = self.gmsh_model.occ.addLine(pt_1, pt_2) + # fixation_line_tag = (1, fixation_line_id) + # embedded_lines_tag_list.append(fixation_line_tag) + # else: + # print( + # f"Applying no-deformation boundary condition at {fp}." + # ) + + # Fragment the lines into the surfaces (equivalent to embedding these lines) + self.gmsh_model.occ.fragment( + this_panel_tag_list, embedded_lines_tag_list + ) - fixation_pts_list = [] - - for k in range(1, num_fixation_pts + 1): - next_pt = k * params.pv_array.span_fixation_pts + for panel_tag in this_panel_tag_list: # 3d cell domain + self.gmsh_model.occ.synchronize() - eps = 1e-9 + # Get the list of 2D surfaces (surfaces) that make up this panel + surf_tags_for_this_panel = self.gmsh_model.getBoundary( + [panel_tag], oriented=False + ) + vol_com = self.gmsh_model.occ.getCenterOfMass(self.ndim, panel_tag[1]) + if np.isclose(vol_com[2], params.pv_array.torque_tube_separation + params.pv_array.panel_thickness / 2.0): + target_key = f"modules" + elif np.isclose(vol_com[2], params.pv_array.torque_tube_separation / 2.0): + target_key = f"connectors" + + + + if target_key is not None: + if target_key in this_panel_transformed_com: + this_panel_transformed_com[target_key].append(vol_com) + else: + this_panel_transformed_com[target_key] = [vol_com] + + for surf_tag in surf_tags_for_this_panel: + surf_dim = surf_tag[0] + surf_id = surf_tag[1] + com = self.gmsh_model.occ.getCenterOfMass(surf_dim, surf_id) + + target_key = None + + surface_located_or_not = False + + # sturctures tagging + if ( + np.isclose(com[0], -half_chord) + ): + target_key = f"panel_left_{panel_ct:.0f}" + surface_located_or_not = True + + if ( + np.isclose(com[0], half_chord) + ): + target_key = f"panel_right_{panel_ct:.0f}" + surface_located_or_not = True + + if ( + np.isclose(com[1], module_distances[0]) + and np.isclose(com[0], 0) and np.isclose(com[2], params.pv_array.torque_tube_separation + params.pv_array.panel_thickness/2.0) + ): + target_key = f"panel_front_{panel_ct:.0f}" + surface_located_or_not = True + + if ( + np.isclose(com[1], module_distances[params.pv_array.modules_per_span]) + and np.isclose(com[0], 0) and np.isclose(com[2], params.pv_array.torque_tube_separation + params.pv_array.panel_thickness/2.0) + ): + target_key = f"panel_back_{panel_ct:.0f}" + surface_located_or_not = True + + for module_id in range(params.pv_array.modules_per_span): if ( - next_pt > eps - and next_pt < params.pv_array.panel_span - eps + com[1] >= module_distances[module_id]+module_span/2.0-params.pv_array.block_chord_div_by_panel_chord * half_chord + and com[1] <= module_distances[module_id]+module_span/2.0+params.pv_array.block_chord_div_by_panel_chord * half_chord + and np.isclose(com[0], 0) and np.isclose(com[2], params.pv_array.torque_tube_separation) ): - fixation_pts_list.append(next_pt) - - else: - fixation_pts_list = params.pv_array.span_fixation_pts - - for fp in fixation_pts_list: - pt_1 = self.gmsh_model.occ.addPoint( - -half_chord, -half_span + fp, -half_thickness - ) - pt_2 = self.gmsh_model.occ.addPoint( - half_chord, -half_span + fp, -half_thickness - ) - - # FIXME: don't add the fixation points into the numpy tagging for now - numpy_pt_list.append( - [ - -half_chord, - -half_span + fp, - -half_thickness, - half_chord, - -half_span + fp, - -half_thickness, - ] - ) - - fixed_pt_id = self.gmsh_model.occ.addLine(pt_1, pt_2) - fixed_pt_tag = (1, fixed_pt_id) + target_key = f"panel_bottom_{panel_ct:.0f}_{module_id:.0f}" + surface_located_or_not = True + break - embedded_lines_tag_list.append(fixed_pt_tag) + if ( + com[1] >= module_distances[module_id]+module_span/2.0-params.pv_array.block_chord_div_by_panel_chord * half_chord + and com[1] <= module_distances[module_id]+module_span/2.0+params.pv_array.block_chord_div_by_panel_chord * half_chord + and np.isclose(com[0], 0) and np.isclose(com[2], params.pv_array.torque_tube_separation + params.pv_array.panel_thickness) + ): + target_key = f"panel_top_{panel_ct:.0f}_{module_id:.0f}" + surface_located_or_not = True + break + # if ( + # np.isclose( + # com[2], + # params.pv_array.torque_tube_separation + # + params.pv_array.panel_thickness + # ) + + # ): + # target_key = f"panel_top_{panel_ct:.0f}" + # surface_located_or_not = True + + # mark the last block in each row + if ( + np.isclose(com[0], params.pv_array.block_chord_div_by_panel_chord * half_chord) + and modeling_torque_tube and np.isclose(com[1], module_distances[params.pv_array.modules_per_span] - params.pv_array.block_chord_div_by_panel_chord * half_chord/2.0) + ): + target_key = f"block_right_{panel_ct:.0f}_{params.pv_array.modules_per_span:.0f}" + surface_located_or_not = True + + if ( + np.isclose(com[0], -params.pv_array.block_chord_div_by_panel_chord * half_chord) + and modeling_torque_tube and np.isclose(com[1], module_distances[params.pv_array.modules_per_span] - params.pv_array.block_chord_div_by_panel_chord * half_chord/2.0) + ): + target_key = f"block_left_{panel_ct:.0f}_{params.pv_array.modules_per_span:.0f}" + surface_located_or_not = True + + if ( + np.isclose(com[2], params.pv_array.torque_tube_separation/2.0) and modeling_torque_tube and np.isclose(com[1], module_distances[params.pv_array.modules_per_span] - params.pv_array.block_chord_div_by_panel_chord * half_chord) + ): + target_key = f"block_front_{panel_ct:.0f}_{params.pv_array.modules_per_span:.0f}" + surface_located_or_not = True + + if ( + np.isclose(com[2], params.pv_array.torque_tube_separation/2.0) and modeling_torque_tube and np.isclose(com[1], module_distances[params.pv_array.modules_per_span]) + ): + target_key = f"block_back_{panel_ct:.0f}_{params.pv_array.modules_per_span:.0f}" + surface_located_or_not = True + + if ( + np.isclose(com[2], 0.0) and modeling_torque_tube and np.isclose(com[1], module_distances[params.pv_array.modules_per_span] - params.pv_array.block_chord_div_by_panel_chord * half_chord/2.0) + ): + target_key = f"block_bottom_{panel_ct:.0f}_{params.pv_array.modules_per_span:.0f}" + surface_located_or_not = True + + if ( + np.isclose(com[2], params.pv_array.torque_tube_separation) and modeling_torque_tube and np.isclose(com[1], module_distances[params.pv_array.modules_per_span] - params.pv_array.block_chord_div_by_panel_chord * half_chord/2.0) + ): + target_key = f"interior_surface_{panel_ct:.0f}" # block/panel interface + surface_located_or_not = True + + # if not the most left/right panel boundary, it is panel/panel interface + if not surface_located_or_not: + for module_id in range(params.pv_array.modules_per_span): + + # sturctures tagging + + if ( + np.isclose(com[1], module_distances[module_id]) + and np.isclose(com[0], 0) and np.isclose(com[2], params.pv_array.torque_tube_separation + params.pv_array.panel_thickness/2.0) + ): + target_key = f"interior_surface_{panel_ct:.0f}" + surface_located_or_not = True + break + + if ( + np.isclose(com[1], module_distances[module_id+1]) + and np.isclose(com[0], 0) and np.isclose(com[2], params.pv_array.torque_tube_separation + params.pv_array.panel_thickness/2.0) + ): + target_key = f"interior_surface_{panel_ct:.0f}" + surface_located_or_not = True + break + + # if ( + # np.isclose(com[2], params.pv_array.torque_tube_separation) + # and np.isclose(com[0], 0.0) and com[1] >= module_distances[module_id]+module_span/2.0-params.pv_array.block_chord_div_by_panel_chord * half_chord/2.0 + # and com[1] <= module_distances[module_id]+module_span/2.0+params.pv_array.block_chord_div_by_panel_chord * half_chord/2.0 + # ): + # target_key = f"panel_bottom_{panel_ct:.0f}" + # surface_located_or_not = True + # break + + if ( + (np.isclose(com[2], params.pv_array.torque_tube_separation) + and np.isclose(com[0], 0.0) and np.isclose(com[1], module_distances[module_id] + params.pv_array.block_chord_div_by_panel_chord * half_chord/2.0) + and modeling_torque_tube) + ): + target_key = f"interior_surface_{panel_ct:.0f}" + surface_located_or_not = True + break + + if ( + np.isclose(com[0], params.pv_array.block_chord_div_by_panel_chord * half_chord) + and modeling_torque_tube and np.isclose(com[1], module_distances[module_id] + params.pv_array.block_chord_div_by_panel_chord * half_chord/2.0) + ): + target_key = f"block_right_{panel_ct:.0f}_{module_id:.0f}" + surface_located_or_not = True + break + + if ( + np.isclose(com[0], -params.pv_array.block_chord_div_by_panel_chord * half_chord) + and modeling_torque_tube and np.isclose(com[1], module_distances[module_id] + params.pv_array.block_chord_div_by_panel_chord * half_chord/2.0) + ): + target_key = f"block_left_{panel_ct:.0f}_{module_id:.0f}" + surface_located_or_not = True + break + + if ( + np.isclose(com[2], params.pv_array.torque_tube_separation/2.0) and modeling_torque_tube and np.isclose(com[1], module_distances[module_id]) + ): + target_key = f"block_front_{panel_ct:.0f}_{module_id:.0f}" + surface_located_or_not = True + break + + if ( + np.isclose(com[2], params.pv_array.torque_tube_separation/2.0) and modeling_torque_tube and np.isclose(com[1], module_distances[module_id] + params.pv_array.block_chord_div_by_panel_chord * half_chord) + ): + target_key = f"block_back_{panel_ct:.0f}_{module_id:.0f}" + surface_located_or_not = True + break + + if ( + np.isclose(com[2], 0.0) and modeling_torque_tube and np.isclose(com[1], module_distances[module_id] + params.pv_array.block_chord_div_by_panel_chord * half_chord/2.0) + ): + target_key = f"block_bottom_{panel_ct:.0f}_{module_id:.0f}" + surface_located_or_not = True + break + + if not surface_located_or_not: + target_key = f"trash_{panel_ct:.0f}" + + if target_key is not None: + if target_key in this_panel_transformed_com: + this_panel_transformed_com[target_key].append(com) + else: + this_panel_transformed_com[target_key] = [com] + + for key, val in this_panel_transformed_com.items(): # all facets for each panel row. + for row_num, com in enumerate(val): + com_array = np.array(com) + + com_array = np.dot(com_array, Ry(tracker_angle_rad).T) + + com_array[0] += xx + com_array[1] += yy + com_array[2] += params.pv_array.elevation + + com_array[0] -= x_center_of_mass + com_array[1] -= y_center_of_mass + + com_array = np.dot(com_array, Rz(array_rotation_rad).T) + + com_array[0] += x_center_of_mass + com_array[1] += y_center_of_mass + + if key in transformed_com: + transformed_com[key].append(com_array) + else: + transformed_com[key] = [com_array] # including all rows + + # actually this is panel array count + panel_ct += 1 - # Store the result of fragmentation, it holds all the small surfaces we need to tag - panel_frags = self.gmsh_model.occ.fragment( - [panel_tag], embedded_lines_tag_list + # Rotate the panel by its tracking angle along the y-axis + # (currently centered at (0.0, 0.0, 0.0)) + self.gmsh_model.occ.rotate( + this_panel_tag_list, + 0.0, + 0.0, + 0.0, + 0, + 1, + 0, + tracker_angle_rad, ) - # extract just the first entry, and remove the 3d entry in position 0 - panel_surfs = panel_frags[0] - panel_surfs.pop(0) - panel_surfs = [k[1] for k in panel_surfs] - - # TODO: USE THESE UNAMBIGUOUS NAMES IN A FUTURE REFACTOR - # self._add_to_domain_markers(f"x_min_{panel_ct:.0f}", [panel_surfs[0]], "facet") - # self._add_to_domain_markers(f"x_max_{panel_ct:.0f}", [panel_surfs[1]], "facet") - # self._add_to_domain_markers(f"y_min_{panel_ct:.0f}", [panel_surfs[2]], "facet") - # self._add_to_domain_markers(f"y_max_{panel_ct:.0f}", [panel_surfs[3]], "facet") - # self._add_to_domain_markers(f"z_min_{panel_ct:.0f}", panel_surfs[4:-1], "facet") - # self._add_to_domain_markers(f"z_max_{panel_ct:.0f}", [panel_surfs[-1]], "facet") - # Translate the panel by (x_center, y_center, elev) self.gmsh_model.occ.translate( - [panel_tag], + this_panel_tag_list, xx, yy, params.pv_array.elevation, ) - # Get the bounding box for this panel - panel_bounding_box = self.gmsh_model.occ.getBoundingBox( - panel_tag[0], panel_tag[1] - ) - - # Store the x_min, y_min, x_max, y_max points - # corresponding to the un-rotated (tracking_angle = 0) configuration - x_min_panel = panel_bounding_box[0] - y_min_panel = panel_bounding_box[1] - z_min_panel = panel_bounding_box[2] - x_max_panel = panel_bounding_box[3] - y_max_panel = panel_bounding_box[4] - z_max_panel = panel_bounding_box[5] - - self.gmsh_model.occ.synchronize() - - # Get the list of 1D surfaces (edges) that make up this panel - surf_tags_for_this_panel = self.gmsh_model.getBoundary( - [panel_tag], oriented=False - ) - - for surf_tag in surf_tags_for_this_panel: - surf_dim = surf_tag[0] - surf_id = surf_tag[1] - com = self.gmsh_model.occ.getCenterOfMass(surf_dim, surf_id) - # sturctures tagging - if np.isclose(com[0], x_min_panel): - self._add_to_domain_markers( - f"left_{panel_ct:.0f}", [surf_id], "facet" - ) - - elif np.isclose(com[0], x_max_panel): - self._add_to_domain_markers( - f"right_{panel_ct:.0f}", [surf_id], "facet" - ) - - elif np.isclose(com[1], y_min_panel): - self._add_to_domain_markers( - f"front_{panel_ct:.0f}", [surf_id], "facet" - ) - - elif np.isclose(com[1], y_max_panel): - self._add_to_domain_markers( - f"back_{panel_ct:.0f}", [surf_id], "facet" - ) - - elif np.isclose(com[2], z_min_panel): - self._add_to_domain_markers( - f"bottom_{panel_ct:.0f}", [surf_id], "facet" - ) - - elif np.isclose(com[2], z_max_panel): - self._add_to_domain_markers( - f"top_{panel_ct:.0f}", [surf_id], "facet" - ) - - panel_ct += 1 - - # Rotate the panel by its tracking angle along the y-axis - # (currently centered at (xx, yy, params.pv_array.elevation)) + # Rotate the panel about the center of the full array as a proxy for changing wind direction (x_center, y_center, 0) self.gmsh_model.occ.rotate( - [panel_tag], - xx, - yy, - params.pv_array.elevation, + this_panel_tag_list, + x_center_of_mass, + y_center_of_mass, 0, - 1, 0, - tracker_angle_rad, + 0, + 1, + array_rotation_rad, ) - # Note that the numpy pt panel array has NOT been shifted, i.e., - # it still represents a panel centered at (0, 0, 0), so the rotation - # can be applied directly about the point (0, 0, 0) without any - # shifting like that which needs to be done above + # Now, apply the same transformations to the numpy representation + # Rotate the panel by its tracking angle along the y-axis + # (currently centered at (0.0, 0.0, 0.0)) numpy_pt_panel_array = np.array(numpy_pt_list) numpy_pt_panel_array = np.reshape(numpy_pt_panel_array, (-1, self.ndim)) @@ -356,27 +656,12 @@ def Rz(theta): numpy_pt_panel_array, Ry(tracker_angle_rad).T ) - # if not hasattr(self, "numpy_pt_array"): - # numpy_pt_array = np.array(numpy_pt_list) - # else: - # numpy_pt_array = np.vcat(numpy_pt_array, np.array(numpy_pt_list)) - + # Translate the panel by (x_center, y_center, elev) numpy_pt_panel_array[:, 0] += xx numpy_pt_panel_array[:, 1] += yy numpy_pt_panel_array[:, 2] += params.pv_array.elevation # Rotate the panel about the center of the full array as a proxy for changing wind direction (x_center, y_center, 0) - self.gmsh_model.occ.rotate( - [panel_tag], - x_center_of_mass, - y_center_of_mass, - 0, - 0, - 0, - 1, - array_rotation_rad, - ) - numpy_pt_panel_array[:, 0] -= x_center_of_mass numpy_pt_panel_array[:, 1] -= y_center_of_mass @@ -394,26 +679,6 @@ def Rz(theta): else: self.numpy_pt_total_array = np.copy(numpy_pt_panel_array) - # Check that this panel still exists in the confines of the domain - bbox = self.gmsh_model.occ.get_bounding_box(panel_tag[0], panel_tag[1]) - - if bbox[0] < params.domain.x_min: - raise ValueError( - f"Panel with location (x, y) = ({xx}, {yy}) extends past x_min wall." - ) - if bbox[1] < params.domain.y_min: - raise ValueError( - f"Panel with location (x, y) = ({xx}, {yy}) extends past y_min wall." - ) - if bbox[3] > params.domain.x_max: - raise ValueError( - f"Panel with location (x, y) = ({xx}, {yy}) extends past x_max wall." - ) - if bbox[4] > params.domain.y_max: - raise ValueError( - f"Panel with location (x, y) = ({xx}, {yy}) extends past y_max wall." - ) - # Fragment all panels from the overall domain self.gmsh_model.occ.fragment(domain_tag_list, panel_tag_list) @@ -423,20 +688,19 @@ def Rz(theta): self.numpy_pt_total_array, (-1, int(2 * self.ndim)) ) - # import matplotlib.pyplot as plt - # for k in self.numpy_pt_total_array: - # plt.plot([k[0], k[3]], [k[1], k[4]]) - # plt.show() - - # it is not necessary to loop over all surfaces, since the panel - # surfaces have been tagged already, but this ensures any ordering change - # doesn't cause problems in the future + # print(transformed_com) + # exit() + + # for all panels + # Loop over all the finalized surfaces after fragmentation and tag everything all_surf_tag_list = self.gmsh_model.occ.getEntities(self.ndim - 1) for surf_tag in all_surf_tag_list: surf_id = surf_tag[1] com = self.gmsh_model.occ.getCenterOfMass(self.ndim - 1, surf_id) + located_this_surface = False + # sturctures tagging if np.isclose(com[0], params.domain.x_min): self._add_to_domain_markers("x_min", [surf_id], "facet") @@ -456,22 +720,88 @@ def Rz(theta): elif np.allclose(com[2], params.domain.z_max): self._add_to_domain_markers("z_max", [surf_id], "facet") + else: + for key, val in transformed_com.items(): + for target_com in val: + # print(target_com) + if np.allclose(np.array(com), target_com): + located_this_surface = True + if "trash" not in key: + self._add_to_domain_markers(key, [surf_id], "facet") + + if not located_this_surface: + print( + f"Warning: Surface {surf_tag} has not been added to domain markers" + ) + + # Since this is not one of the exterior walls, we should check if it extends + # past the boundaries x_min, x_max, ... + this_surf_bbox = self.gmsh_model.occ.get_bounding_box( + self.ndim - 1, surf_id + ) + + # Test that the rotated point still exists in the box domain + if this_surf_bbox[0] < params.domain.x_min: + raise ValueError(f"A panel extends past the x_min wall.") + if this_surf_bbox[0] > params.domain.x_max: + raise ValueError(f"A panel extends past the x_max wall.") + if this_surf_bbox[1] < params.domain.y_min: + raise ValueError(f"A panel extends past the y_min wall.") + if this_surf_bbox[1] > params.domain.y_max: + raise ValueError(f"A panel extends past the y_max wall.") + if this_surf_bbox[2] < 0.0: + raise ValueError( + f"A panel extends past the z_min wall (ground level = 0.0)." + ) + if this_surf_bbox[2] > params.domain.z_max: + raise ValueError(f"A panel extends past the z_max wall.") + + + # mark the panel and connector volumes + all_vol_tag_list = self.gmsh_model.occ.getEntities(self.ndim) + + for vol_tag in all_vol_tag_list: + vol_id = vol_tag[1] + com = self.gmsh_model.occ.getCenterOfMass(self.ndim, vol_id) + + located_this_volume = False + + for key, val in transformed_com.items(): + for target_com in val: + # print(target_com) + if np.allclose(np.array(com), target_com): + located_this_volume = True + if "trash" not in key: + self._add_to_domain_markers(key, [vol_id], "cell") + + # Mark the fluid volume # Volumes are the entities with dimension equal to the mesh dimension vol_tag_list = self.gmsh_model.occ.getEntities(self.ndim) structure_vol_list = [] fluid_vol_list = [] + num_solids = len(vol_tag_list) + for vol_tag in vol_tag_list: vol_id = vol_tag[1] - if vol_id <= params.pv_array.stream_rows * params.pv_array.span_rows: + if vol_id <= num_solids - 1: # Solid Cell + # Since all panel/table structures were created after the box domain + # the final fragmentation removes the original vol_tag of the box and appends + # it to the end, thus, everything with id < num_solids is structure, and + # vol_tag = num_solids (the last-added, fragmented box) is fluid structure_vol_list.append(vol_id) else: # Fluid Cell fluid_vol_list.append(vol_id) - self._add_to_domain_markers("structure", structure_vol_list, "cell") + structure_panel_only_list = [] + structure_connector_only_list = [] +# + # # self._add_to_domain_markers("structure", structure_vol_list, "cell") + # for individual_vol_id in structure_vol_list: + # self._add_to_domain_markers(f"structure_{individual_vol_id:03.0f}", [individual_vol_id], "cell") self._add_to_domain_markers("fluid", fluid_vol_list, "cell") # Record all the data collected to domain_markers as physics groups with physical names @@ -1144,24 +1474,35 @@ def set_length_scales(self, params, domain_markers): internal_surface_tags = [] for panel_id in range(params.pv_array.stream_rows * params.pv_array.span_rows): - internal_surface_tags.append( - domain_markers[f"bottom_{panel_id}"]["gmsh_tags"][0] + + internal_surface_tags.extend( + domain_markers[f"panel_left_{panel_id}"]["gmsh_tags"] ) - internal_surface_tags.append( - domain_markers[f"top_{panel_id}"]["gmsh_tags"][0] + internal_surface_tags.extend( + domain_markers[f"panel_right_{panel_id}"]["gmsh_tags"] ) - internal_surface_tags.append( - domain_markers[f"left_{panel_id}"]["gmsh_tags"][0] + internal_surface_tags.extend( + domain_markers[f"panel_front_{panel_id}"]["gmsh_tags"] ) - internal_surface_tags.append( - domain_markers[f"right_{panel_id}"]["gmsh_tags"][0] + internal_surface_tags.extend( + domain_markers[f"panel_back_{panel_id}"]["gmsh_tags"] ) - internal_surface_tags.append( - domain_markers[f"front_{panel_id}"]["gmsh_tags"][0] + + for module_id in range(params.pv_array.modules_per_span): + internal_surface_tags.extend( + domain_markers[ + f"panel_bottom_{panel_id:.0f}_{module_id:.0f}" + ]["gmsh_tags"] ) - internal_surface_tags.append( - domain_markers[f"back_{panel_id}"]["gmsh_tags"][0] + internal_surface_tags.extend( + domain_markers[ + f"panel_top_{panel_id:.0f}_{module_id:.0f}" + ]["gmsh_tags"] ) + + # print(internal_surface_tags) + # print(len(internal_surface_tags)) + # exit() min_dist = [] diff --git a/pvade/structure/ElasticityAnalysis.py b/pvade/structure/ElasticityAnalysis.py index 4a84d98..e585054 100644 --- a/pvade/structure/ElasticityAnalysis.py +++ b/pvade/structure/ElasticityAnalysis.py @@ -87,6 +87,285 @@ def build_boundary_conditions(self, domain, params): self.bc = build_structure_boundary_conditions(domain, params, self.V) + def calculate_K_for_Robin_BC(self, domain, flow, params): + + # shear modulus of tube + Gt = params.structure.elasticity_modulus_tube / (2 * (1 + params.structure.poissons_ratio_tube)) + + Ipt = np.pi/2*(params.pv_array.torque_tube_outer_radius**4 - params.pv_array.torque_tube_inner_radius**4) + + Gp = params.structure.elasticity_modulus / (2 * (1 + params.structure.poissons_ratio)) + + Ipp = 1/3.0*params.pv_array.panel_thickness * params.pv_array.panel_chord**3 + + # total number of pv rows + total_num_panels = params.pv_array.stream_rows * params.pv_array.span_rows + + for panel_id in range(total_num_panels): + # construct the matrix, size: modules_per_span+1 by modules_per_span+1 + theo_matrix = np.zeros((params.pv_array.modules_per_span+1, params.pv_array.modules_per_span+1)) + theo_vector = np.zeros(params.pv_array.modules_per_span+1) + + connector_locations = np.linspace( + 0, params.pv_array.panel_span, params.pv_array.modules_per_span + 1 + ) + + # ?? double integral is not available until the second time step. Should make it be zero for the + # first time step then it will be updated once the first step is finished? but the assemble is out of + # the time loop. Should we solve fluid then structure? + + total_torque_on_this_panel_name = f"total_torque_panel_{panel_id:.0f}" + total_torque_on_this_panel = getattr(flow, total_torque_on_this_panel_name).value + + theo_matrix[params.pv_array.modules_per_span, :] = 1/Gp/Ipp + theo_vector[params.pv_array.modules_per_span] = -total_torque_on_this_panel/Gp/Ipp + + # contribution from panel: + + # contribution from C0 to matrix: + + theo_matrix[:params.pv_array.modules_per_span, :params.pv_array.fixed_location] += connector_locations[params.pv_array.fixed_location]/(Gp*Ipp) + theo_matrix[:params.pv_array.modules_per_span, 1:params.pv_array.fixed_location] -= connector_locations[1:params.pv_array.fixed_location]/(Gp*Ipp) + + # contribution from C0 to vector: + T_double_integral_at_fixed_location_name = f"double_integral_total_torque_panel_{panel_id:.0f}_{params.pv_array.fixed_location:.0f}" + T_double_integral_at_fixed_location = getattr(flow, T_double_integral_at_fixed_location_name).value + theo_vector[:params.pv_array.modules_per_span] += -T_double_integral_at_fixed_location/Gp/Ipp + + # contribution of panels to matrix + theo_matrix[0, 0] += -connector_locations[0]/(Gp*Ipp) + theo_matrix[1, 0] += -connector_locations[1]/(Gp*Ipp) + + for i in range(2, params.pv_array.fixed_location): + theo_matrix[i, :i] += -connector_locations[i]/(Gp*Ipp) + theo_matrix[i, 1:i] += connector_locations[1:i]/(Gp*Ipp) + + for i in range(params.pv_array.fixed_location, params.pv_array.modules_per_span): + theo_matrix[i, :i+1] += -connector_locations[i+1]/(Gp*Ipp) + theo_matrix[i, 1:i+1] += connector_locations[1:i+1]/(Gp*Ipp) + + # contribution of panels to vector + T_double_integral_array = [] + for i in range(params.pv_array.modules_per_span+1): + name = f"double_integral_total_torque_panel_{panel_id:.0f}_{i:.0f}" + T_double_integral_array.append(getattr(flow, name).value) + T_double_integral_array = np.array(T_double_integral_array) + theo_vector[:params.pv_array.modules_per_span] += np.delete(T_double_integral_array, params.pv_array.fixed_location)/Gp/Ipp + + # contribution from tube: + + # contribution from C0 to matrix: + + theo_matrix[:params.pv_array.modules_per_span, params.pv_array.fixed_location:params.pv_array.modules_per_span+1] += -connector_locations[params.pv_array.fixed_location]/(Gt*Ipt) + theo_matrix[:params.pv_array.modules_per_span, 1:params.pv_array.fixed_location] += -connector_locations[1:params.pv_array.fixed_location]/(Gt*Ipt) + + # contribution from C0 to vector: + theo_vector[:params.pv_array.modules_per_span] += total_torque_on_this_panel*connector_locations[params.pv_array.fixed_location]/Gt/Ipt + + # contribution of tube to matrix + theo_matrix[0, 1:params.pv_array.modules_per_span+1] += connector_locations[0]/(Gt*Ipt) + theo_matrix[1, 1:params.pv_array.modules_per_span+1] += connector_locations[1]/(Gt*Ipt) + + for i in range(2, params.pv_array.fixed_location): + theo_matrix[i, i:params.pv_array.modules_per_span+1] += connector_locations[i]/(Gt*Ipt) + theo_matrix[i, 1:i] += connector_locations[1:i]/(Gt*Ipt) + + for i in range(params.pv_array.fixed_location, params.pv_array.modules_per_span): + theo_matrix[i, i+1:params.pv_array.modules_per_span+1] += connector_locations[i+1]/(Gt*Ipt) + theo_matrix[i, 1:i+1] += connector_locations[1:i+1]/(Gt*Ipt) + + # contribution of tube to vector + theo_vector[:params.pv_array.fixed_location] += -total_torque_on_this_panel/Gt/Ipt*connector_locations[:params.pv_array.fixed_location] + theo_vector[params.pv_array.fixed_location:params.pv_array.modules_per_span] += -total_torque_on_this_panel/Gt/Ipt*connector_locations[params.pv_array.fixed_location] + + R_reaction_torque = np.dot(np.linalg.inv(theo_matrix), theo_vector) # this is the torque applied to tube + + # rotation of tube at each connector location + tube_rotate_matrix = np.zeros((params.pv_array.modules_per_span, params.pv_array.modules_per_span+1)) + tube_rotate_vector = np.zeros(params.pv_array.modules_per_span) + + # contribution from panel: + + # contribution from C0 to matrix: + + tube_rotate_matrix[:params.pv_array.modules_per_span, :params.pv_array.fixed_location] += connector_locations[params.pv_array.fixed_location]/(Gp*Ipp) + tube_rotate_matrix[:params.pv_array.modules_per_span, 1:params.pv_array.fixed_location] -= connector_locations[1:params.pv_array.fixed_location]/(Gp*Ipp) + + # contribution from C0 to vector: + tube_rotate_vector[:params.pv_array.modules_per_span] += T_double_integral_array[params.pv_array.fixed_location]/Gp/Ipp + + + # contribution of panels to matrix + tube_rotate_matrix[0, 0] += -connector_locations[0]/(Gp*Ipp) + tube_rotate_matrix[1, 0] += -connector_locations[1]/(Gp*Ipp) + + for i in range(2, params.pv_array.fixed_location): + tube_rotate_matrix[i, :i] += -connector_locations[i]/(Gp*Ipp) + tube_rotate_matrix[i, 1:i] += connector_locations[1:i]/(Gp*Ipp) + + for i in range(params.pv_array.fixed_location, params.pv_array.modules_per_span): + tube_rotate_matrix[i, :i+1] += -connector_locations[i+1]/(Gp*Ipp) + tube_rotate_matrix[i, 1:i+1] += connector_locations[1:i+1]/(Gp*Ipp) + + # contribution of panels to vector + tube_rotate_vector[:params.pv_array.modules_per_span] += -np.delete(T_double_integral_array, params.pv_array.fixed_location)/Gp/Ipp + + + phi = np.dot(tube_rotate_matrix, R_reaction_torque) + tube_rotate_vector + + if isinstance(params.pv_array.tracker_angle, list): + if panel_id == 0: + assert ( + len(params.pv_array.tracker_angle) + == params.pv_array.stream_rows * params.pv_array.span_rows + ), f"Length of tracker angle list ({len(params.pv_array.tracker_angle)}) not equal to total number of PV tables ({params.pv_array.stream_rows * params.pv_array.span_rows})." + + tracker_angle_rad = np.radians( + params.pv_array.tracker_angle[panel_id] + ) + else: + tracker_angle_rad = np.radians(params.pv_array.tracker_angle) + + if np.linalg.norm(phi) == 0: + K = np.zeros(params.pv_array.modules_per_span) + else: + K = np.abs(12*(np.delete(R_reaction_torque, params.pv_array.fixed_location))/((params.pv_array.block_chord_div_by_panel_chord * params.pv_array.panel_chord)**3)/np.cos(tracker_angle_rad+phi)/(np.sin(tracker_angle_rad+phi)-np.sin(tracker_angle_rad))/(params.pv_array.block_chord_div_by_panel_chord * params.pv_array.panel_chord/2)) + + # assume K is 0 at the fixed connector + K = np.flip(np.insert(K, params.pv_array.fixed_location, 0)) + + # K[0] is the stiffness at the most back connector (highest y) + # K[10] is the stiffness at the most front connector (lowest y) + # for block_bot_surface, it is numbered from lowest y to highest y, so flip K + + # if there are 10 modules per array, there are 11 connectors, K has shape of 10, the K at the fixed connector is not included. + for i in range(params.pv_array.modules_per_span+1): + + name_K = f"spring_stiffness_{panel_id:.0f}_{i:.0f}" + setattr(self, name_K, K[i]) + + # num_panel_right_fixed = params.pv_array.modules_per_span // 2 + # num_panel_left_fixed = params.pv_array.modules_per_span - num_panel_right_fixed + + # # for right fixed part: + + # for panel_id in range(total_num_panels): + # total_torque_right_fixed = 0 + # total_torque_left_fixed = 0 + + # for i in range(num_panel_left_fixed, params.pv_array.modules_per_span): + # name = f"total_torque_panel_{panel_id:.0f}_{i:.0f}" + # total_torque_right_fixed += getattr(flow, name) + + # for i in range(num_panel_left_fixed): + # name = f"total_torque_panel_{panel_id:.0f}_{i:.0f}" + # total_torque_left_fixed += getattr(flow, name) + + # theo_matrix_right_fixed = np.zeros((num_panel_right_fixed+1, num_panel_right_fixed+1)) + # theo_vector_right_fixed = np.zeros(num_panel_right_fixed+1) + + # theo_matrix_right_fixed[0, :] = 1.0 + # theo_vector_right_fixed[0] = total_torque_right_fixed + # theo_matrix_right_fixed[1:, :-1] = params.pv_array.panel_span/params.pv_array.modules_per_span/Gt/Ipt + + # # as the double integral of torque along the span from front to back, left fixed part to right fixed part, it need to be flipped + # T_double_integral_right_fixed = [] + # T_right_fixed = [] + # for i in range(params.pv_array.modules_per_span): + # name = f"double_integral_total_torque_panel_{panel_id:.0f}_{params.pv_array.modules_per_span-1-i:.0f}" + # T_double_integral_right_fixed.append(getattr(flow, name)) + # name = f"total_torque_panel_{panel_id:.0f}_{params.pv_array.modules_per_span-1-i:.0f}" + # T_right_fixed.append(getattr(flow, name)) + # for i in range(num_panel_right_fixed): + # for j in range(i+1): + # theo_vector_right_fixed[i+1] += T_double_integral_right_fixed[-1-j]/(Gp*Ipp) + # theo_vector_right_fixed[i+1] -= (i+1-j)*T_right_fixed[-1-j]/(Gp*Ipp)*params.pv_array.panel_span/params.pv_array.modules_per_span + # for j in range(i): + # theo_matrix_right_fixed[i+1, :-1-j-1] += params.pv_array.panel_span/params.pv_array.modules_per_span/(Gt*Ipt) + + + # for i in range(num_panel_right_fixed): + # for j in range(i+1): + # theo_matrix_right_fixed[i+1, num_panel_right_fixed-j:] -= params.pv_array.panel_span/params.pv_array.modules_per_span/Gp/Ipp + + # R_reaction_torque = np.dot(np.linalg.inv(theo_matrix_right_fixed), theo_vector_right_fixed) + + # tube_rotate_matrix = np.ones((num_panel_right_fixed, num_panel_right_fixed))*params.pv_array.panel_span/params.pv_array.modules_per_span*(1.0/Gt/Ipt) + + # for i in range(num_panel_right_fixed): + # for j in range(i): + # tube_rotate_matrix[i, :num_panel_right_fixed-1-j] += params.pv_array.panel_span/params.pv_array.modules_per_span*(1.0/Gt/Ipt) + + # # check the standalone code, why it need to be flipped. + # phi = np.flip(np.dot(tube_rotate_matrix, R_reaction_torque[:-1])) + + # # rotation of connector + # x_panel = np.arange(num_panel_right_fixed)*(params.pv_array.panel_span/params.pv_array.modules_per_span) # location of panels + + # block_length = params.pv_array.block_chord_div_by_panel_chord * params.pv_array.panel_chord + # block_width = params.pv_array.block_chord_div_by_panel_chord*params.pv_array.panel_span/params.pv_array.modules_per_span/2 + + # # To Do: check the angle, is this correct? + # array_rotation = (params.fluid.wind_direction + 90.0) % 360.0 + # array_rotation_rad = np.radians(array_rotation) + + # for i in range(num_panel_right_fixed): + + # phi_i = phi[i] + + # K_i = 12*(R_reaction_torque[i])/((block_length)**3)/np.cos(array_rotation_rad+phi_i)/(np.sin(array_rotation_rad+phi_i)-np.sin(array_rotation_rad))/block_width + + # name_K = f"spring_stiffness_{panel_id:.0f}_{params.pv_array.modules_per_span+1-i:.0f}" + # setattr(self, name_K, K_i) + + # # for left fixed part: + # theo_matrix_left_fixed = np.zeros((num_panel_left_fixed+1, num_panel_left_fixed+1)) + # theo_vector_left_fixed = np.zeros(num_panel_left_fixed+1) + # theo_matrix_left_fixed[0, :] = 1.0 + # theo_vector_left_fixed[0] = total_torque_left_fixed + + # T_double_integral_left_fixed = [] + # T_left_fixed = [] + + # for i in range(num_panel_left_fixed): + # name = f"double_integral_total_torque_panel_{panel_id:.0f}_{num_panel_left_fixed-1-i:.0f}" + # T_double_integral_left_fixed.append(getattr(flow, name)) + # name = f"total_torque_panel_{panel_id:.0f}_{num_panel_left_fixed-1-i:.0f}" + # T_left_fixed.append(getattr(flow, name)) + + # for i in range(num_panel_left_fixed): + # for j in range(i+1): + # theo_matrix_left_fixed[i+1, j] = (i+1-j)*params.pv_array.panel_span/params.pv_array.modules_per_span/(Gp*Ipp) + # theo_vector_left_fixed[i+1] += T_double_integral_left_fixed[j]/(Gp*Ipp) + # for j in range(i): + # theo_vector_left_fixed[i+1] += (i-j)*T_left_fixed[j]/(Gp*Ipp)*params.pv_array.panel_span/params.pv_array.modules_per_span + + # for i in range(num_panel_left_fixed): + # theo_matrix_left_fixed[i+1:, i+1:] -= params.pv_array.panel_span/params.pv_array.modules_per_span/Gt/Ipt + + # R_reaction_torque = np.dot(np.linalg.inv(theo_matrix_left_fixed), theo_vector_left_fixed) # this is the torque applied to tube + + + # tube_rotate_matrix = np.zeros((num_panel_left_fixed, num_panel_left_fixed)) + # for i in range(num_panel_left_fixed): + # tube_rotate_matrix[i:, i:] += params.pv_array.panel_span/params.pv_array.modules_per_span/Gt/Ipt + + # phi = np.dot(tube_rotate_matrix, R_reaction_torque[1:]) + + # # rotation of connector + # x_panel = np.arange(num_panel_left_fixed)*(params.pv_array.panel_span/params.pv_array.modules_per_span)+params.pv_array.panel_span/params.pv_array.modules_per_span # location of panles + + + # for i in range(num_panel_left_fixed): + + # phi_i = phi[i] + + # K_i = 12*(R_reaction_torque[i+1])/((block_length)**3)/np.cos(array_rotation_rad+phi_i)/(np.sin(array_rotation_rad+phi_i)-np.sin(array_rotation_rad))/block_width + + # name_K = f"spring_stiffness_{panel_id:.0f}_{num_panel_left_fixed-1-i:.0f}" + + # setattr(self, name_K, K_i) + def update_a(self, u, u_old, v_old, a_old, dt, beta, ufl=True): # Update formula for acceleration # a = 1/(2*beta)*((u - u0 - v0*dt)/(0.5*dt*dt) - (1-2*beta)*a0) @@ -126,7 +405,7 @@ def update_fields(self, u, u_old, v_old, a_old, dt, beta, gamma): def avg(self, x_old, x_new, alpha): return alpha * x_old + (1 - alpha) * x_new - def build_forms(self, domain, params, structure): + def build_forms(self, domain, params, structure, flow): """Builds all variational statements This method creates all the functions, expressions, and variational @@ -212,6 +491,9 @@ def c(u, u_): def k_nominal(u, u_): return ufl.inner(P_(u), ufl.grad(u_)) + + def k_nominal_connector(u, u_): + return ufl.inner(P_connector(u), ufl.grad(u_)) # The deformation gradient, F = I + dy/dX def F_(u): @@ -241,6 +523,17 @@ def S_(u): S_svk = structure.lame_lambda * ufl.tr(E) * I + 2.0 * structure.lame_mu * E return S_svk + + # The second Piola–Kirchhoff stress, S + def S_connector(u): + E = E_(u) + I = ufl.Identity(len(u)) + + # return lamda * ufl.tr(E) * I + 2.0 * mu * (E - ufl.tr(E) * I / 3.0) + # TODO: Why does the above form give a better result and where does it come from? + + S_svk = structure.lame_lambda_connector * ufl.tr(E) * I + 2.0 * structure.lame_mu_connector * E + return S_svk # The first Piola–Kirchhoff stress tensor, P = F*S def P_(u): @@ -248,6 +541,12 @@ def P_(u): S = S_(u) # return ufl.inv(F) * S return F * S + + def P_connector(u): + F = F_(u) + S = S_connector(u) + # return ufl.inv(F) * S + return F * S # self.uh_exp = dolfinx.fem.Function(self.V, name="Deformation") @@ -302,15 +601,35 @@ def P_(u): F = ufl.grad(self.u) + ufl.Identity(len(self.u)) J = ufl.det(F) + + self.z_unit_vector = dolfinx.fem.Constant(domain.structure.msh, [0.0,0.0,1.0]) # surface traction, N/m^2 + + self.calculate_K_for_Robin_BC(domain, flow, params) + + dx_structure = ufl.Measure( + "dx", domain=domain.structure.msh, subdomain_data=domain.structure.cell_tags + ) + + + # To Do: how to differentiate the connector part and the panel part, dx_connector and dx_panel self.res = ( - m(self.avg(self.a_old, a_new, self.alpha_m), self.u_) * ufl.dx - + c(self.avg(self.v_old, v_new, self.alpha_f), self.u_) * ufl.dx - + k_nominal(self.avg(self.u_old, self.u, self.alpha_f), self.u_) * ufl.dx - - structure.rho * ufl.inner(self.f, self.u_) * ufl.dx + m(self.avg(self.a_old, a_new, self.alpha_m), self.u_) * dx_structure + + c(self.avg(self.v_old, v_new, self.alpha_f), self.u_) * dx_structure + + k_nominal(self.avg(self.u_old, self.u, self.alpha_f), self.u_) * dx_structure(domain.domain_markers["modules"]["idx"]) + + k_nominal_connector(self.avg(self.u_old, self.u, self.alpha_f), self.u_) * dx_structure(domain.domain_markers["connectors"]["idx"]) + - structure.rho * ufl.inner(self.f, self.u_) * dx_structure(domain.domain_markers["modules"]["idx"]) + - structure.rho_connector * ufl.inner(self.f, self.u_) * dx_structure(domain.domain_markers["connectors"]["idx"]) - ufl.dot(ufl.dot(self.stress_predicted * J * ufl.inv(F.T), n), self.u_) * self.ds ) # - Wext(self.u) + # Robin boundary condition terms + for panel_id in range(params.pv_array.stream_rows * params.pv_array.span_rows): + for i in range(params.pv_array.modules_per_span+1): + name_K = f"spring_stiffness_{panel_id:.0f}_{i:.0f}" + K_springs = dolfinx.fem.Constant(domain.structure.msh, float(getattr(self, name_K))) + self.res -= ufl.dot(K_springs * self.u_, self.z_unit_vector)*self.ds(domain.domain_markers[f"block_bottom_{panel_id:.0f}_{i:.0f}"]["idx"]) + # self.a = dolfinx.fem.form(ufl.lhs(res)) # self.L = dolfinx.fem.form(ufl.rhs(res)) @@ -484,26 +803,52 @@ def solve(self, params, dataIO, structure): try: idx = structure.north_east_corner_dofs[0] - # idx = self.north_east_corner_dofs[0] - nw_corner_accel = self.u.x.array[ + north_east_corner_deformation = self.u.x.array[ structure.ndim * idx : structure.ndim * idx + structure.ndim ].astype(np.float64) - print(nw_corner_accel) + print(f"Deformation: {north_east_corner_deformation}") + except: - nw_corner_accel = np.zeros(structure.ndim, dtype=np.float64) + north_east_corner_deformation = np.zeros(structure.ndim, dtype=np.float64) - nw_corner_accel_global = np.zeros( + try: + idx = structure.north_east_corner_dofs[0] + north_east_corner_acceleration = self.a_old.x.array[ + structure.ndim * idx : structure.ndim * idx + structure.ndim + ].astype(np.float64) + print(f"Acceleration: {north_east_corner_acceleration}") + + except: + north_east_corner_acceleration = np.zeros(structure.ndim, dtype=np.float64) + + # Initialize a buffer to collect everything into + north_east_corner_deformation_global = np.zeros( + (self.num_procs, structure.ndim), dtype=np.float64 + ) + + north_east_corner_acceleration_global = np.zeros( (self.num_procs, structure.ndim), dtype=np.float64 ) - self.comm.Gather(nw_corner_accel, nw_corner_accel_global, root=0) + # Gather all points (many of which are zeros) to rank 0 + self.comm.Gather( + north_east_corner_deformation, north_east_corner_deformation_global, root=0 + ) - # print(f"Acceleration at North West corner = {nw_corner_accel}") + self.comm.Gather( + north_east_corner_acceleration, + north_east_corner_acceleration_global, + root=0, + ) if self.rank == 0: - norm2 = np.sum(nw_corner_accel_global**2, axis=1) + norm2 = np.sum(north_east_corner_deformation_global**2, axis=1) + max_norm2_idx = np.argmax(norm2) + np_deformation = north_east_corner_deformation_global[max_norm2_idx, :] + + norm2 = np.sum(north_east_corner_acceleration_global**2, axis=1) max_norm2_idx = np.argmax(norm2) - np_accel = nw_corner_accel_global[max_norm2_idx, :] + np_acceleration = north_east_corner_acceleration_global[max_norm2_idx, :] accel_pos_filename = os.path.join( params.general.output_dir_sol, "accel_pos.csv" @@ -512,18 +857,35 @@ def solve(self, params, dataIO, structure): if self.first_call_to_solver: with open(accel_pos_filename, "w") as fp: - fp.write("#x-pos,y-pos,z-pos\n") if structure.ndim == 3: - fp.write(f"{np_accel[0]},{np_accel[1]},{np_accel[2]}\n") + fp.write( + "#x-deformation,y-deformation,z-deformation,x-acceleration,y-acceleration,z-acceleration\n" + ) + fp.write( + f"{np_deformation[0]},{np_deformation[1]},{np_deformation[2]}," + ) + fp.write( + f"{np_acceleration[0]},{np_acceleration[1]},{np_acceleration[2]}\n" + ) elif structure.ndim == 2: - fp.write(f"{np_accel[0]},{np_accel[1]}\n") + fp.write( + "#x-deformation,y-deformation,x-acceleration,y-acceleration\n" + ) + fp.write(f"{np_deformation[0]},{np_deformation[1]},") + fp.write(f"{np_acceleration[0]},{np_acceleration[1]}\n") else: with open(accel_pos_filename, "a") as fp: if structure.ndim == 3: - fp.write(f"{np_accel[0]},{np_accel[1]},{np_accel[2]}\n") + fp.write( + f"{np_deformation[0]},{np_deformation[1]},{np_deformation[2]}," + ) + fp.write( + f"{np_acceleration[0]},{np_acceleration[1]},{np_acceleration[2]}\n" + ) elif structure.ndim == 2: - fp.write(f"{np_accel[0]},{np_accel[1]}\n") + fp.write(f"{np_deformation[0]},{np_deformation[1]},") + fp.write(f"{np_acceleration[0]},{np_acceleration[1]}\n") if self.first_call_to_solver: self.first_call_to_solver = False diff --git a/pvade/structure/StructureMain.py b/pvade/structure/StructureMain.py index 6e55265..8f81ea2 100644 --- a/pvade/structure/StructureMain.py +++ b/pvade/structure/StructureMain.py @@ -80,6 +80,10 @@ def __init__(self, domain, params): domain.structure.msh, params.structure.rho ) # Constant(0.) + self.rho_connector = dolfinx.fem.Constant( + domain.structure.msh, params.structure.rho_tube + ) + # Define structural properties self.E = params.structure.elasticity_modulus # 1.0e9 self.poissons_ratio = params.structure.poissons_ratio # 0.3 @@ -89,6 +93,16 @@ def __init__(self, domain, params): * self.poissons_ratio / ((1.0 + self.poissons_ratio) * (1.0 - 2.0 * self.poissons_ratio)) ) + + # we are assuming the connectors has same mechanical properties as the tubes + self.E_connector = params.structure.elasticity_modulus_tube # 1.0e9 + self.poissons_ratio_connector = params.structure.poissons_ratio_tube # 0.3 + self.lame_mu_connector = self.E_connector / (2.0 * (1.0 + self.poissons_ratio_connector)) + self.lame_lambda_connector = ( + self.E_connector + * self.poissons_ratio_connector + / ((1.0 + self.poissons_ratio_connector) * (1.0 - 2.0 * self.poissons_ratio_connector)) + ) if self.rank == 0: print( @@ -110,9 +124,21 @@ def _north_east_corner(x): else: tracker_angle_rad = np.radians(params.pv_array.tracker_angle) - x1 = 0.5 * params.pv_array.panel_chord * np.cos(tracker_angle_rad) - x2 = 0.5 * params.pv_array.panel_thickness * np.sin(tracker_angle_rad) - corner = [x1 - x2, 0.5 * params.pv_array.panel_span] + """ + y + + ^ + | + |---o NE corner (measured on bottom of panel) + | | + | | + |-----> x + """ + + corner = [ + 0.5 * params.pv_array.panel_chord * np.cos(tracker_angle_rad), + 0.5 * params.pv_array.panel_span, + ] east_edge = np.logical_and(corner[0] - eps < x[0], x[0] < corner[0] + eps) north_edge = np.logical_and(corner[1] - eps < x[1], x[1] < corner[1] + eps) @@ -121,12 +147,12 @@ def _north_east_corner(x): return north_east_corner - north_east_corner_facets = dolfinx.mesh.locate_entities_boundary( + north_east_corner_vertices = dolfinx.mesh.locate_entities_boundary( domain.structure.msh, 0, _north_east_corner ) self.north_east_corner_dofs = dolfinx.fem.locate_dofs_topological( - self.elasticity.V, 0, north_east_corner_facets + self.elasticity.V, 0, north_east_corner_vertices ) def build_boundary_conditions(self, domain, params): @@ -189,7 +215,7 @@ def update_fields(self, u, u_old, v_old, a_old, dt, beta, gamma): def avg(self, x_old, x_new, alpha): return alpha * x_old + (1 - alpha) * x_new - def build_forms(self, domain, params): + def build_forms(self, domain, params, flow): """Builds all variational statements This method creates all the functions, expressions, and variational @@ -205,7 +231,7 @@ def build_forms(self, domain, params): params (:obj:`pvade.Parameters.SimParams`): A SimParams object """ - self.elasticity.build_forms(domain, params, self) + self.elasticity.build_forms(domain, params, self, flow) def _assemble_system(self, params): """Pre-assemble all LHS matrices and RHS vectors diff --git a/pvade/structure/boundary_conditions.py b/pvade/structure/boundary_conditions.py index 0aa0ea0..956b857 100644 --- a/pvade/structure/boundary_conditions.py +++ b/pvade/structure/boundary_conditions.py @@ -272,7 +272,7 @@ def build_structure_boundary_conditions(domain, params, functionspace): total_num_panels = params.pv_array.stream_rows * params.pv_array.span_rows for num_panel in range(total_num_panels): - for location in params.structure.bc_list: + for location in params.structure.bc_list: # it is empty location_panel = f"{location}_{num_panel}" # f"front_{num_panel}" , f"back_{num_panel}": # for location in [f"left_{num_panel}"]:# , f"right_{num_panel}": @@ -344,20 +344,12 @@ def connection_point_up_helper(nodes_to_pin_between): return fn_handle - # Start pinning along the lines expressed byt numpy_pt_total_array - # First, determine the total number of rows (number of lines to pin) - num_nodes = np.shape(domain.numpy_pt_total_array)[0] - - # Determine how many pinning lines exist per each panel (e.g., 24 lines distributed on 8 panels means 3 lines per panel) - nodes_per_panel = int(num_nodes / total_num_panels) - - # The torque tube entry (oriented spanwise along the middle, divides panel into upstream and downstream rectangular halves) - # is always the first entry, e.g., [0, ..., ..., 3, ..., ..., 6, ..., ...], [0, 3, 6] are the torque tubes - tube_nodes_idx = np.arange(0, num_nodes, nodes_per_panel, dtype=np.int64) + # # Start pinning along the lines expressed byt numpy_pt_total_array + # The center line of connectors bottom surface is fixed to remove rigid body motion if params.structure.tube_connection == True: # If making torque tube connections, pass only those pinning lines to the BC identification function - tube_nodes = domain.numpy_pt_total_array[tube_nodes_idx, :] + tube_nodes = domain.numpy_pt_total_array[:, :] facet_uppoint = dolfinx.mesh.locate_entities( domain.structure.msh, 1, connection_point_up_helper(tube_nodes) @@ -369,18 +361,18 @@ def connection_point_up_helper(nodes_to_pin_between): bc.append(dolfinx.fem.dirichletbc(zero_vec, dofs_disp, functionspace)) if params.structure.motor_connection == True: - # If making motor mount connections, pass only those pinning lines to the BC identification function - # this is done by making a copy of the numpy_pt_total_array with the torque tube lines *deleted* - # not done in place, so numpy_pt_total_array remains unaltered. - motor_nodes = np.delete(domain.numpy_pt_total_array, tube_nodes_idx, axis=0) + # The bottom surface of the center connector is fixed to represent the motor mount + motor_location = params.pv_array.fixed_location # fixed location along the span, if fixed_location=5, it means the left boundary of the 6th connector is fixed + for panel_id in range(total_num_panels): - facet_uppoint = dolfinx.mesh.locate_entities( - domain.structure.msh, 1, connection_point_up_helper(motor_nodes) - ) - dofs_disp = dolfinx.fem.locate_dofs_topological( - functionspace, 1, [facet_uppoint] + mount_facet = domain.structure.facet_tags.find( + domain.domain_markers[f"block_left_{panel_id:.0f}_{motor_location:.0f}"]["idx"] ) + + dofs_disp = dolfinx.fem.locate_dofs_topological( + functionspace, 2, [mount_facet] + ) - bc.append(dolfinx.fem.dirichletbc(zero_vec, dofs_disp, functionspace)) + bc.append(dolfinx.fem.dirichletbc(zero_vec, dofs_disp, functionspace)) return bc diff --git a/pvade/tests/input/yaml/embedded_box.yaml b/pvade/tests/input/yaml/embedded_box.yaml index 0bd12d5..ad13c66 100644 --- a/pvade/tests/input/yaml/embedded_box.yaml +++ b/pvade/tests/input/yaml/embedded_box.yaml @@ -9,15 +9,15 @@ domain: x_max: 1.0 y_min: -1.0 y_max: 1.0 - z_min: -1.0 - z_max: 1.0 + z_min: 0.0 + z_max: 2.0 l_char: 0.05 pv_array: stream_rows: 1 span_rows: 1 stream_spacing: 1.0 span_spacing: 1.0 - elevation: 0.0 + elevation: 0.5 # 0.0 panel_chord: 1.2 panel_span: 1.1 panel_thickness: 1.0 diff --git a/pvade/tests/input/yaml/sim_params.yaml b/pvade/tests/input/yaml/sim_params.yaml index 7b6e4a6..2e63349 100644 --- a/pvade/tests/input/yaml/sim_params.yaml +++ b/pvade/tests/input/yaml/sim_params.yaml @@ -16,11 +16,11 @@ domain: pv_array: stream_rows: 7 span_rows: 1 - elevation: 1.5 + elevation: 1.45 # 1.5 - 0.5*0.1 = 1.4 stream_spacing: 7.0 panel_chord: 2.0 panel_span: 7.0 - panel_thickness: 0.1 + panel_thickness: 0.1 tracker_angle: -30.0 solver: dt: 0.005 diff --git a/pvade/tests/test_fsi_mesh.py b/pvade/tests/test_fsi_mesh.py index 2d4f475..eba906e 100644 --- a/pvade/tests/test_fsi_mesh.py +++ b/pvade/tests/test_fsi_mesh.py @@ -148,10 +148,10 @@ def test_meshing_3dpanels_rotations(wind_direction, num_stream_rows, num_span_ro # Create the 4 corners of this table corresponding to the *top* surface (+0.5*thickness) top_surface_corners = np.array( [ - [-0.5 * chord, -0.5 * span, 0.5 * thickness], - [0.5 * chord, -0.5 * span, 0.5 * thickness], - [0.5 * chord, 0.5 * span, 0.5 * thickness], - [-0.5 * chord, 0.5 * span, 0.5 * thickness], + [-0.5 * chord, -0.5 * span, thickness], + [0.5 * chord, -0.5 * span, thickness], + [0.5 * chord, 0.5 * span, thickness], + [-0.5 * chord, 0.5 * span, thickness], ] ) diff --git a/pvade/tests/test_mesh_movement.py b/pvade/tests/test_mesh_movement.py index 4de08c8..ec0798d 100644 --- a/pvade/tests/test_mesh_movement.py +++ b/pvade/tests/test_mesh_movement.py @@ -49,7 +49,7 @@ def test_calc_distance_to_panel_surface(): dx = params.domain.x_max - 0.5 * params.pv_array.panel_chord dy = params.domain.y_max - 0.5 * params.pv_array.panel_span - dz = params.domain.z_max - 0.5 * params.pv_array.panel_thickness + dz = params.pv_array.elevation truth_max_dist = np.sqrt(dx * dx + dy * dy + dz * dz) assert np.isclose(max_dist, truth_max_dist) @@ -113,9 +113,7 @@ def __init__(self, domain, x_shift, y_shift, z_shift): assert np.isclose( np.amin(structure_coords_before[:, 1]), -0.5 * params.pv_array.panel_span ) - assert np.isclose( - np.amin(structure_coords_before[:, 2]), -0.5 * params.pv_array.panel_thickness - ) + assert np.isclose(np.amin(structure_coords_before[:, 2]), params.pv_array.elevation) assert np.isclose( np.amax(structure_coords_before[:, 0]), 0.5 * params.pv_array.panel_chord @@ -124,7 +122,8 @@ def __init__(self, domain, x_shift, y_shift, z_shift): np.amax(structure_coords_before[:, 1]), 0.5 * params.pv_array.panel_span ) assert np.isclose( - np.amax(structure_coords_before[:, 2]), 0.5 * params.pv_array.panel_thickness + np.amax(structure_coords_before[:, 2]), + params.pv_array.elevation + params.pv_array.panel_thickness, ) # Move the mesh by the amount prescribed in u_delta diff --git a/pvade/tests/test_solve.py b/pvade/tests/test_solve.py index 8734207..77d5457 100644 --- a/pvade/tests/test_solve.py +++ b/pvade/tests/test_solve.py @@ -310,7 +310,7 @@ def test_fsi2(): # print('pos_data = ', pos_data) - assert np.allclose(pos_data, pos_data_truth) + assert np.allclose(pos_data[:, 0:2], pos_data_truth) print(lift_and_drag_data) # assert np.allclose(lift_and_drag_data[:, 0:3], lift_and_drag_data_truth[:, 0:3]) # needs new truth values to pass, mesh has changed diff --git a/pvade_main.py b/pvade_main.py index 8bc322a..5e52d56 100644 --- a/pvade_main.py +++ b/pvade_main.py @@ -16,6 +16,10 @@ import os from mpi4py import MPI +""" +to run: python -u pvade_main.py --input_file=input/duramat_case_study.yaml --domain.l_char=3.0 --pv_array.torque_tube_radius=0.1 --pv_array.torque_tube_separation=0.2 --pv_array.stream_rows=2 --pv_array.tracker_angle=[20,20] --pv_array.modules_per_span=10 +""" + def main(input_file=None): # Get the path to the input file from the command line @@ -40,7 +44,7 @@ def main(input_file=None): domain.read_mesh_files(params.general.input_mesh_dir, params) else: domain.build(params) - + # If we only want to create the mesh, we can stop here if params.general.mesh_only: list_timings(params.comm, [TimingType.wall]) @@ -60,7 +64,7 @@ def main(input_file=None): structure = Structure(domain, params) structure.build_boundary_conditions(domain, params) # # # Build the fluid forms - structure.build_forms(domain, params) + structure.build_forms(domain, params, flow) else: structure = None @@ -75,6 +79,7 @@ def main(input_file=None): for k in range(params.solver.t_steps): current_time = (k + 1) * params.solver.dt + if ( structural_analysis and (k + 1) % solve_structure_interval_n == 0