From d2787f90df1d9121b5fae7c4263766b3a7e45f52 Mon Sep 17 00:00:00 2001 From: cflin Date: Thu, 15 Jun 2023 21:58:22 +0800 Subject: [PATCH] todo: fix dT_drho --- .../config.json | 51 +++++++++++-------- .../top-thermoelastic-regular-model/main.cpp | 2 +- src/ThermoelasticTop3d.cpp | 49 ++++++++++-------- src/ThermoelasticTop3d.h | 2 +- 4 files changed, 61 insertions(+), 43 deletions(-) diff --git a/examples/top-thermoelastic-regular-model/config.json b/examples/top-thermoelastic-regular-model/config.json index ac1b932..5434b63 100644 --- a/examples/top-thermoelastic-regular-model/config.json +++ b/examples/top-thermoelastic-regular-model/config.json @@ -1,34 +1,34 @@ { "material": { - "E": 1.0, + "E": 2e11, "poisson_ratio": 0.3, - "thermal_conductivity":0.004, "unit": "W/K", - "thermal_expansion_coefficient": 1e-6, "unit": "1/K" + "thermal_conductivity":43, "unit": "W/(m*K)", + "thermal_expansion_coefficient": 1.21e-5, "unit": "1/K" }, "topology": { "max_loop": 100, "volfrac": 0.5, "r_min": 1.5, "T_ref": 295, - "T_limit": 325, - "R_E": 0, - "R_lambda": 0, + "T_limit": 345, + "R_E": 12, + "R_lambda": 8, "R_beta": 0 }, "model": { "regular_model": { - "lx": 2, - "ly": 20, - "lz": 30 + "lx": 1, + "ly": 50, + "lz": 35 }, "visual_model":"2-box-refined.obj" }, "mechanical_boundary_condition":{ "NBC": [ { - "min": [0, 1, 0], - "max": [1, 1, 0], - "val": [0.0, 0.0, -1.0] + "min": [0, 0.5, 0.5], + "max": [1, 0.5, 0.5], + "val": [0.0, 0.0, -1e1] } ], @@ -37,29 +37,38 @@ "min": [0, 0, 0], "max": [1, 0, 1], "dir": [1, 1, 1] + }, + { + "min": [0, 1, 0], + "max": [1, 1, 1], + "dir": [1, 1, 1] } ] }, "thermal_boundary_condition": { "NBC": [ { - "min": [0, 0.25, 0.5], - "max": [1, 0.25, 0.5], - "heat_flux": 1e-15, "unit": "W" + "min": [0, 0.2, 0.5], + "max": [1, 0.2, 0.5], + "heat_flux": 0.335, "unit": "W" }, { - "min": [0, 0.75, 0.5], - "max": [1, 0.75, 0.5], - "heat_flux": 1e-15, "unit": "W" + "min": [0, 0.8, 0.5], + "max": [1, 0.8, 0.5], + "heat_flux": 0.335, "unit": "W" } ], "DBC": [ { "min": [0, 0, 0], - "max": [1, 0, 1], - "temperature": 312, "unit": "K" - + "max": [1, 1, 0], + "temperature": 325, "unit": "K" + }, + { + "min": [0, 0, 1], + "max": [1, 1, 1], + "temperature": 325, "unit": "K" } ] } diff --git a/examples/top-thermoelastic-regular-model/main.cpp b/examples/top-thermoelastic-regular-model/main.cpp index ae6ed05..217d795 100644 --- a/examples/top-thermoelastic-regular-model/main.cpp +++ b/examples/top-thermoelastic-regular-model/main.cpp @@ -133,6 +133,6 @@ int main() { write_tensor3d(output_dir / "txt" / "top_mach_regular_field_matrix.txt", ten_rho, sp_mech_mesh->GetOrigin(), sp_mech_mesh->GetOrigin() + sp_mech_mesh->GetLenBox()); WriteTensorToVtk(output_dir / "vtk" / "top_mach_regular_field_matrix.vtk", ten_rho, sp_mech_mesh); -// WriteUToVtk(output_dir / "vtk" / "top_thermoelastic_regular_U.vtk",ther_top3d.GetU(),sp_mech_mesh); + WriteUToVtk(output_dir / "vtk" / "top_thermoelastic_regular_U.vtk",ther_top3d.GetU(),sp_mech_mesh); } diff --git a/src/ThermoelasticTop3d.cpp b/src/ThermoelasticTop3d.cpp index 5c152a6..d4af1c9 100644 --- a/src/ThermoelasticTop3d.cpp +++ b/src/ThermoelasticTop3d.cpp @@ -25,9 +25,11 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { int loop = 0; double change = 1.0; double E0_m = sp_mech_top3d_->sp_fea_->sp_material_->E; + double E_min = E0_m * sp_mech_top3d_->sp_para_->E_factor; double lambda0 = sp_mech_top3d_->sp_fea_->sp_material_->thermal_conductivity; double lambda_min = lambda0 * sp_mech_top3d_->sp_para_->E_factor; double alpha0 = sp_mech_top3d_->sp_fea_->sp_material_->thermal_expansion_coefficient; + double alpha_min=alpha0 * sp_mech_top3d_->sp_para_->E_factor; // Precompute Eigen::VectorXd dv(sp_mesh_->GetNumEles()); @@ -75,7 +77,7 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { spdlog::warn("using Eigen built-in direct solver!"); #endif // start iteration - while (change > sp_para_->tol_x && loop < sp_para_->max_loop) { + while (change > sp_para_->tol_x*1e-6 && loop < sp_para_->max_loop) { ++loop; // filter xPhys_col = sp_mech_top3d_->H_ * (xPhys_col.array() / sp_mech_top3d_->Hs_.array()).matrix().eval(); @@ -94,10 +96,10 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { return (1 + R) / down.pow(2); }; auto CalE_Vec = [&](const Eigen::VectorXd &vec_rho) -> Eigen::VectorXd { - return CalR_Vec(vec_rho, sp_para_->R_E) * E0_m; + return E_min + CalR_Vec(vec_rho, sp_para_->R_E).array() * (E0_m - E_min); }; auto CalDEDrho_Vec = [&](const Eigen::VectorXd &vec_rho) -> Eigen::VectorXd { - return CalDRDrho_Vec(vec_rho, sp_para_->R_E) * E0_m; + return CalDRDrho_Vec(vec_rho, sp_para_->R_E) * (E0_m - E_min); }; auto CalLambda_Vec = [&](const Eigen::VectorXd &vec_rho) -> Eigen::VectorXd { return lambda_min + CalR_Vec(vec_rho, sp_para_->R_lambda).array() * (lambda0 - lambda_min); @@ -106,14 +108,14 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { return CalDRDrho(rho, sp_para_->R_lambda) * (lambda0 - lambda_min); }; auto CalDlambdaDrho_Vec = [&](const Eigen::VectorXd &vec_rho) -> Eigen::VectorXd { - return CalDRDrho_Vec(vec_rho, sp_para_->R_lambda) * (lambda0 - lambda_min); + return CalDRDrho_Vec(vec_rho, sp_para_->R_lambda) * (lambda0 - lambda_min); }; auto CalBeta = [&](double rho) { - return CalR(rho, sp_para_->R_beta) * E0_m * alpha0; + return E_min*alpha_min+ CalR(rho, sp_para_->R_beta) * (E0_m * alpha0-E_min*alpha_min); }; auto CalDBetaDrho = [&](double rho) { - return CalDRDrho(rho, sp_para_->R_beta) * E0_m * alpha0; + return CalDRDrho(rho, sp_para_->R_beta) * (E0_m * alpha0-E_min*alpha_min); }; // solve T @@ -150,6 +152,7 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { double beta_rho = CalBeta(xPhys_col(i)); F_th(dofs_m) += beta_rho * (Te - sp_mech_top3d_->sp_para_->T_ref) * Inted_; } + spdlog::info("||Fth|| / ||Fm||: {}",F_th.norm()/sp_mech_top3d_->F_.norm()); Eigen::VectorXd F = Eigen::VectorXd(sp_mech_top3d_->F_) + F_th; sp_mech_top3d_->IntroduceFixedDofs(sp_mech_top3d_->K_, F); @@ -191,7 +194,7 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { Eigen::VectorXd ele_dFth_drho = CalDBetaDrho(xPhys_col(i)) * (Te - sp_mech_top3d_->sp_para_->T_ref) * Inted_;// 24x1 assert(ele_dFth_drho.size() == 24); - v_dFth_drho(Eigen::seqN(i * ele_dFth_drho.rows(), ele_dFth_drho.size())) = ele_dFth_drho; + v_dFth_drho(Eigen::seqN(i * ele_dFth_drho.size(), ele_dFth_drho.size())) = ele_dFth_drho; } auto v_dFth_drho_tri = Vec2Triplet(i_dFth_drho_, j_dFth_drho_, v_dFth_drho); dFth_drho.setFromTriplets(v_dFth_drho_tri.begin(), v_dFth_drho_tri.end()); @@ -209,7 +212,7 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { Eigen::MatrixXd ele_dFth_dT = Eigen::VectorXd::Ones(dofs_th.size()) * 1.0 / 8.0 * beta_rho * Inted_.transpose(); assert(ele_dFth_dT.rows() == 8 && ele_dFth_dT.cols() == 24); - v_dFth_dT(Eigen::seqN(i * ele_dFth_dT.rows(), ele_dFth_dT.size())) = ele_dFth_dT.reshaped(); + v_dFth_dT(Eigen::seqN(i * ele_dFth_dT.size(), ele_dFth_dT.size())) = ele_dFth_dT.reshaped(); } auto v_dFth_dT_tri = Vec2Triplet(i_dFth_dT_, j_dFth_dT_, v_dFth_dT); dFth_dT.setFromTriplets(v_dFth_dT_tri.begin(), v_dFth_dT_tri.end()); @@ -236,12 +239,16 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { } Eigen::VectorXd lambda_t_Mul_dKt_drho_Mul_T = CalDlambdaDrho_Vec(xPhys_col).array() * ce_th.array(); // dc_drho - Eigen::VectorXd dc_drho = lambda_t_Mul_dKt_drho_Mul_T + - lambda_m_Mul_dKm_drho_Mul_U + - 2 * Eigen::VectorXd(dF_drho * sp_mech_top3d_->U_); -// Eigen::VectorXd dc_drho = -// lambda_m_Mul_dKm_drho_Mul_U + -// 2 * Eigen::VectorXd(dF_drho * sp_mech_top3d_->U_); +// Eigen::VectorXd dc_drho = lambda_t_Mul_dKt_drho_Mul_T + +// lambda_m_Mul_dKm_drho_Mul_U + +// 2 * Eigen::VectorXd(dF_drho * sp_mech_top3d_->U_); + Eigen::VectorXd dc_drho = + lambda_t_Mul_dKt_drho_Mul_T + + + lambda_m_Mul_dKm_drho_Mul_U + + + 2 * Eigen::VectorXd(dF_drho * sp_mech_top3d_->U_) + ; // dT_drho Eigen::MatrixXd dT_drho = Eigen::MatrixXd::Zero(sp_thermal_top3d_->sp_mesh_->GetNumEles(), sp_thermal_top3d_->set_dofs_to_load.size()); @@ -285,10 +292,11 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { // dc_dx Eigen::VectorXd dc_dx = sp_mech_top3d_->drho_dx_ * dc_drho; - // dT_dx - Eigen::MatrixXd dT_dx = sp_mech_top3d_->drho_dx_ * dT_drho; + // dT_dx + Eigen::MatrixXd dT_dx = sp_mech_top3d_->drho_dx_ * dT_drho; // mma solver +#define SENSITIVITY_SCALE_COEF 100 size_t num_constraints = 1 + dT_dx.cols();// volume and temperature constraints @@ -298,14 +306,15 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { Eigen::VectorXd variables_tmp = flg_chosen ? xPhys_col(chosen_ele_id) : xPhys_col; double f0val = c; Eigen::VectorXd df0dx = flg_chosen - ? dc_dx(chosen_ele_id).eval() / dc_dx(chosen_ele_id).cwiseAbs().maxCoeff() - : dc_dx / dc_dx.cwiseAbs().maxCoeff(); + ? dc_dx(chosen_ele_id).eval() / dc_dx(chosen_ele_id).cwiseAbs().maxCoeff() *SENSITIVITY_SCALE_COEF + : dc_dx / dc_dx.cwiseAbs().maxCoeff() *SENSITIVITY_SCALE_COEF; // double fval = v - num_variables * sp_para_->volfrac; + Eigen::VectorXd fval = (Eigen::VectorXd(num_constraints) << (v / num_variables - sp_para_->volfrac), - T(v_dof).array() / sp_para_->T_limit - 1).finished(); + T(v_dof).array() / sp_para_->T_limit - 1).finished() * SENSITIVITY_SCALE_COEF; // Eigen::VectorXd dfdx = flg_chosen ? dv(chosen_ele_id) : dv; Eigen::MatrixXd dfdx = (Eigen::MatrixXd(num_variables, num_constraints) - << 1.0 / num_variables * dv, 1.0 / sp_para_->T_limit * dT_dx).finished().transpose(); + << 1.0 / num_variables * dv , 1.0 / sp_para_->T_limit * dT_dx).finished().transpose() *SENSITIVITY_SCALE_COEF; static Eigen::VectorXd low_bounds = Eigen::VectorXd::Zero(num_variables); diff --git a/src/ThermoelasticTop3d.h b/src/ThermoelasticTop3d.h index 5d95745..7c51b27 100644 --- a/src/ThermoelasticTop3d.h +++ b/src/ThermoelasticTop3d.h @@ -37,7 +37,7 @@ namespace da::sha::top { i_dFth_drho_ = (Eigen::VectorXi::LinSpaced(sp_mech_top3d_->sp_mesh_->GetNumEles(), 0, sp_mech_top3d_->sp_mesh_->GetNumEles()) * - Eigen::RowVectorXi::Ones(sp_mech_top3d_->sp_mesh_->Get_DOFS_EACH_ELE())).reshaped(); + Eigen::RowVectorXi::Ones(sp_mech_top3d_->sp_mesh_->Get_DOFS_EACH_ELE())).transpose().reshaped(); j_dFth_drho_ = sp_mech_top3d_->sp_mesh_->GetEleId2DofsMap().transpose().reshaped(); D0_ = sp_mech_top3d_->sp_fea_->computeD(1.0);