Dear everyone
I meet a problem when adding the multi-solid regions with diffe…rent solid materials by modifying the fsi3 case.
<img width="1075" alt="两相板几何" src="https://user-images.githubusercontent.com/97109440/222377634-8d4f7dc3-dfad-4900-9993-cf51452391cb.png">
The tag of boundaries and regions set in gmsh are:
```
1 10 "Solidinterface"
1 11 "Inlet"
1 12 "Outlet"
1 13 "Wall"
1 14 "Bar"
1 15 "Barwall"
1 16 "Circle"
2 17 "Fluid" #region
2 18 "Solid1" #region
2 19 "Solid2" #region
```
To do it, I referred to the `turtleFSI-master/turtleFSI/modules/domains.py and solid.py` and changed them into `domain_solid.py and solid_multiphase.py` as follows:
```
#domain_solid.py
from turtleFSI.modules import *
from dolfin import Constant, inner, inv, grad, div
def assign_solid_domain_properties(dx, dx_s_id, rho_s, nu_s, mu_s, lambda_s, solid_properties, domains, **namespace):
"""
Assigns solid properties to each region.
"""
# 1. Create differential for each solid region, and organize into solid_properties list of dicts
dx_s = {}
if isinstance(dx_s_id, list): # If dx_s_id is a list (i.e, if there are multiple solid regions):
for solid_region in range(len(dx_s_id)):
dx_s[solid_region] = dx(dx_s_id[solid_region], subdomain_data=domains) # Create dx_s for each solid region
dx_s_id_list=dx_s_id
else:
dx_s[0] = dx(dx_s_id, subdomain_data=domains)
dx_s_id_list=[dx_s_id]
# Assign material properties to each region
if len(solid_properties) == 0:
if isinstance(dx_s_id, list): # If dx_s_id is a list (i.e, if there are multiple solid regions):
for solid_region in range(len(dx_s_id)):
solid_properties.append(
{"dx_s_id": dx_s_id[solid_region],
"rho_s": rho_s[solid_region],
"nu_s": nu_s[solid_region],
"mu_s": mu_s[solid_region],
"lambda_s": lambda_s[solid_region]})
else:
solid_properties.append(
{"dx_s_id": dx_s_id,
"rho_s": rho_s,
"nu_s": nu_s,
"mu_s": mu_s,
"lambda_s": lambda_s})
elif isinstance(solid_properties, dict):
solid_properties = [solid_properties]
return dict(dx_s=dx_s, dx_s_id_list=dx_s_id_list, solid_properties=solid_properties)
```
and
```
#solid_multiphase.py
from turtleFSI.modules import *
from dolfin import Constant, inner, grad
def solid_setup(d_, v_, phi, psi, dx_s, dx_s_id_list, solid_properties, k, theta,
gravity, **namespace):
delta = 1E7
# Theta scheme constants
theta0 = Constant(theta)
theta1 = Constant(1 - theta)
F_solid_linear = 0
F_solid_nonlinear = 0
for solid_region in range(len(dx_s_id_list)):
rho_s = solid_properties[solid_region]["rho_s"]
nu_s = solid_properties[solid_region]["nu_s"]
mu_s = solid_properties[solid_region]["mu_s"]
lambda_s = solid_properties[solid_region]["lambda_s"]
# Temporal term and convection
F_solid_linear = (rho_s/k * inner(v_["n"] - v_["n-1"], psi)*dx_s[solid_region]
+ delta * rho_s * (1 / k) * inner(d_["n"] - d_["n-1"], phi) * dx_s[solid_region]
- delta * rho_s * inner(theta0 * v_["n"] + theta1 * v_["n-1"], phi) * dx_s[solid_region])
# Gravity
if gravity is not None:
F_solid_linear -= inner(Constant((0, -gravity * rho_s)), psi)*dx_s[solid_region]
# Stress
F_solid_nonlinear = theta0 * inner(Piola1(d_["n"], lambda_s, mu_s), grad(psi)) * dx_s[solid_region]
F_solid_linear += theta1 * inner(Piola1(d_["n-1"], lambda_s, mu_s), grad(psi)) * dx_s[solid_region]
return dict(F_solid_linear=F_solid_linear, F_solid_nonlinear=F_solid_nonlinear)
```
The input parameters are
```
# Temporal variables
T=10, # End time [s]
dt=0.001, # Time step [s]
theta=0.51, # Temporal scheme
# Physical constants ('FSI 3')
Um=2.0, # Max. velocity inlet, CDF3: 2.0 [m/s]
rho_f=1.0e3, # Fluid density [kg/m3]
mu_f=1.0, # Fluid dynamic viscosity [Pa.s]
solid_properties=[{"dx_s_id": 18, "rho_s": 2.0E3, "nu_s": 0.8, "mu_s": 4.0e6, "lambda_s": 8E6}, {"dx_s_id": 19, "rho_s": 1.0E3, "nu_s": 0.4, "mu_s": 2.0e6, "lambda_s": 4e6}],
# Domain settings
dx_f_id=17, # Domain id of the fluid domain
dx_s_id=[18,19], # Domain id of the solid domain
fluid="fluid", # ["fluid", "no_fluid"] Turn off fluid and only solve the solid problem
solid="solid_multiphase", # ["solid", "no_solid"] Turn off solid and only solve the fluid problem
# Problem specific
folder="CylinderFlap_multisolidphase_interface_results", # Name of the results folder
extrapolation="biharmonic", # No displacement to extrapolate
extrapolation_sub_type="constrained_disp", # Biharmonic type
```
But the results are bad:
<img width="586" alt="结果" src="https://user-images.githubusercontent.com/97109440/222380126-64f400b2-145d-4f86-82fb-a54c8431fe4d.png">
It seems the code cannot rightly recognize the physical region, do some friends here know why? Is that because of the interface between two solid regions?
Thanks in advance!
Best regards
ZHANG