Three-Dimensional Cochlea Model in ABAQUS

by Liu-Jie Ren

In this note, we present the details in developing 3D finite element cochlea models in commercial software ABAQUS. The flow chart is as follows

In [9]:
%%html
<img src='./cochlea_model_3d_abaqus_flowchart.svg'>

Model 0: Air-conduction test model

Note, since the at first time, MODEL1 encoungered some problems. To debug the model, we use the simplest first model (model0) to simulate the air conduction traveling wave. The technical details are shown in model 1

MODEL0 contains the following parts

  • fluid domain
    • scala vestibula
    • scala tympani
    • vestibule
    • helicotrema
  • solid domain
    • oval window
    • round window
    • bony lamina
    • basilar membrane

Working procedure: model00.hm -> model00.inp --> dictory MODEL00 with nodes, elements and other files --> spiral shape --> generate basilar membrane --> calculate in abaqus

Straight model

In [179]:
system('echo MODEL00>__dic.log');
system('gfortran readAbaqus3D.f90');
system('a');
system('del a.exe');

 ----------------- ABAQUS READER ----------------- 
 | Working dictory >>> MODEL00 
 | Input file name >>> MODEL00.inp 
 ------------------------------------------------- 
A subdirectory or file MODEL00 already exists. 
 | Creating folder >>> MODEL00 
 *NODE                                                                                                                                                                                
 *ELEMENT,TYPE=S4R,ELSET=shell-ow 
  ./MODEL00/model_ele_001.ele 
 *ELEMENT,TYPE=S4R,ELSET=shell-rw 
  ./MODEL00/model_ele_002.ele 
 *ELEMENT,TYPE=S4R,ELSET=shell-lam 
  ./MODEL00/model_ele_003.ele 
 *ELEMENT,TYPE=C3D8R,ELSET=fluid_sv_others 
  ./MODEL00/model_ele_004.ele 
 *ELEMENT,TYPE=C3D8R,ELSET=fluid_sv_bm 
  ./MODEL00/model_ele_005.ele 
 *ELEMENT,TYPE=C3D8R,ELSET=fluid_sv_heli 
  ./MODEL00/model_ele_006.ele 
 *ELEMENT,TYPE=C3D8R,ELSET=fluid_vest 
  ./MODEL00/model_ele_007.ele 
 *ELEMENT,TYPE=C3D8R,ELSET=fluid_st_others 
  ./MODEL00/model_ele_008.ele 
 *ELEMENT,TYPE=C3D8R,ELSET=fluid_st_bm 
  ./MODEL00/model_ele_009.ele 
 *ELEMENT,TYPE=C3D8R,ELSET=fluid_st_heli 
  ./MODEL00/model_ele_010.ele 
 *NSET, NSET=line4, UNSORTED 
  ./MODEL00/model_bc_001.nset 
 *NSET, NSET=line3, UNSORTED 
  ./MODEL00/model_bc_002.nset 
 *NSET, NSET=line2, UNSORTED 
  ./MODEL00/model_bc_003.nset 
 *NSET, NSET=line1, UNSORTED 
  ./MODEL00/model_bc_004.nset 
 *NSET, NSET=rw_fixed, UNSORTED 
  ./MODEL00/model_bc_005.nset 
 *NSET, NSET=ow_fixed, UNSORTED 
  ./MODEL00/model_bc_006.nset 
 *NSET, NSET=ow_center 
  ./MODEL00/model_bc_007.nset 
 *ELSET, ELSET=ow_inner 
  ./MODEL00/model_bc_008.elset 
 *ELSET, ELSET=ow_outer 
  ./MODEL00/model_bc_009.elset 
 *SURFACE, NAME = fluid-surf, TYPE = ELEMENT 
  ./MODEL00/model_surf_001.surf
In [180]:
% tospiral(1) make a spiral shape, tospiral(0) keeps a strainght model
tospiral(0)

Mission Complete
In [182]:
folder = textread('__dic.log','%s'); folder = folder{1};
line1 = 'model_bc_004.nset';
line2 = 'model_bc_003.nset';
line3 = 'model_bc_002.nset';
line4 = 'model_bc_001.nset';
generateBM(folder,line1,line2,line3,line4);

The length of the BM mid-line is about 34.900000 mm
Mission Complete !
In [185]:
%%file MODEL00/model00-ac.inp
*Heading
*preprint, echo=yes, mode=yes, history=no, contact=no
** --------------------------------------------------------
**          M O D E L - D E F I N I T I O N                
** --------------------------------------------------------
*part, name=part1
*end part
*assembly, name=asse1
*instance, name=inst1, part=part1
** node and elements
*node
*include, input=model_node_spiral.cnn
*include, input=input-inst1-node-bm.inp

*ELEMENT, TYPE=AC3D8, ELSET=FLUID
*include, input=model_ele_004.ele
*include, input=model_ele_005.ele
*include, input=model_ele_006.ele
*include, input=model_ele_007.ele
*include, input=model_ele_008.ele
*include, input=model_ele_009.ele
*include, input=model_ele_010.ele

*ELEMENT, TYPE=S4R, ELSET=OW
*include, input=model_ele_001.ele
*ELEMENT, TYPE=S4R, ELSET=RW
*include, input=model_ele_002.ele
*ELEMENT, TYPE=S4R, ELSET=WALL02
*include, input=model_ele_003.ele

*ELEMENT, TYPE=S4R, ELSET=BM
*include, input=input-inst1-elem-bm.inp
** nodal thickness
*nodal thickness
*include, input=input-distribution-thick.inp
** bm-element-mats
** *shell section, elset=bm, material=mat-bm, nodal thickness
*include, input=input-inst1-bm-set.inp
*include, input=input-inst1-bm-sec.inp
** sections
*SHELL SECTION, ELSET=WALL02, MATERIAL=MAT-BONE
1,
*SHELL SECTION, ELSET=OW, MATERIAL=MAT-OW
1,
*SHELL SECTION, ELSET=RW, MATERIAL=MAT-RW
1,
*SOLID SECTION, ELSET=FLUID, MATERIAL=MAT-FLUID
,
** sets and elements and surfaces
*SURFACE, NAME=SURF-FLUID, TYPE=ELEMENT
*include, input=model_surf_001.surf
*SURFACE, NAME=SURF-SOLID, TYPE=ELEMENT
OW,     SNEG
RW,     SNEG
BM,     SPOS
WALL02, SPOS
*SURFACE, NAME=OW-PRESSURE, TYPE=ELEMENT
OW,     SPOS 
*NSET, NSET=OW-CENTER
*include, input=model_bc_007.nset
*ELSET, ELSET=WALL
WALL02
*NSET, NSET=WALLNODE, ELSET=WALL02
*NSET, NSET=BM-MIDLINE
*include, input=input-inst1-bm-mid.inp
*NSET, NSET=RW-FIXED
*include, input=model_bc_005.nset
*RIGID BODY, REF NODE=OW-CENTER, ELSET=OW
*end instance
*tie, name=tie1, adjust=yes, position tolerance=0.01, no thickness
inst1.surf-fluid, inst1.surf-solid
*end assembly
** --------------------------------------------------------
**          M O D E L - S E T T I N G S                    
** --------------------------------------------------------
** MATERIALS
*include, input=input-inst1-bm-mat.inp
*material, name=mat-bone
*Damping, structural=0.05
*density
2.200e-9
*elastic
14100
*material, name=mat-ow
*Damping, structural=0.05
*density
1.200e-9
*elastic
1.0, 0.3
*material, name=mat-rw
*Damping, structural=0.05
*density
1.200e-9
*elastic
1.0, 0.3
*material, name=mat-fluid
*acoustic medium
1E6,
*acoustic medium, volumetric drag
0.0, 500.
*density
1.000e-9
** BOUNDARIES
*boundary
INST1.OW-CENTER, 1 ,1
INST1.OW-CENTER, 2 ,2
INST1.OW-CENTER, 4 ,4
INST1.OW-CENTER, 5 ,5
INST1.OW-CENTER, 6 ,6
INST1.WALLNODE, 1, 1
INST1.WALLNODE, 2, 2
INST1.WALLNODE, 3, 3
INST1.WALLNODE, 4, 4
INST1.WALLNODE, 5, 5
INST1.WALLNODE, 6, 6
INST1.RW-FIXED, 1, 1
INST1.RW-FIXED, 2, 2
INST1.RW-FIXED, 3, 3
INST1.RW-FIXED, 4, 4
INST1.RW-FIXED, 5, 5
INST1.RW-FIXED, 6, 6
** STEP CALC
*step, name=step001, nlgeom=no, perturbation
*Steady State Dynamics, direct, friction damping=NO
2000, 2000, 1, 1.
*BOUNDARY, REAL
INST1.OW-CENTER,3,3,-1.0
*BOUNDARY, imaginary
INST1.OW-CENTER,3,3, 0.0
** LOAD
**Dsload, real
**inst1.press, P, 1.0
*output, field, variable=preselect
*output, history, variable=preselect
*end step

Created file 'C:\Users\RenLi\__matlab\__model_cochlea\3d_model_abaqus\MODEL00\model00-ac.inp'.
In [194]:
%%html
<center><img width='1200px' src='./MODEL00/straight-ac-2000Hz.gif'> straight model</center>
straight model

Spiral model

In [197]:
% tospiral(1) make a spiral shape, tospiral(0) keeps a strainght model
tospiral(1)

Mission Complete
In [198]:
folder = textread('__dic.log','%s'); folder = folder{1};
line1 = 'model_bc_004.nset';
line2 = 'model_bc_003.nset';
line3 = 'model_bc_002.nset';
line4 = 'model_bc_001.nset';
generateBM(folder,line1,line2,line3,line4);

The length of the BM mid-line is about 35.220971 mm
Mission Complete !
In [200]:
!copy MODEL00/model00-ac.inp MODEL00/model00-ac-spiral.inp
In [207]:
%%html
<center><img width='1200px' src='./MODEL00/spiral-ac-2000Hz.gif'> spiral model</center>

<p>Problems:<br/>

1. The spiral shape interact with each at the base<br/>
2. Tip ends has element distortion<br/>

These problems will be corrected later.
</p>
spiral model

Problems:
1. The spiral shape interact with each at the base
2. Tip ends has element distortion
These problems will be corrected later.

Model 1: Rectangular Section

Here we developped the most commonly used cochlea box model, with two sections (the scala vestibula and the scala tympany).

Note:

  • At first, model1 encounted problems, so we develop model 0 first to verify the set-ups from the very beginning.

Step 01: Define the model in HyperMesh

note that the model was in unit: mm in HyperMesh, we should revise to unit: mm-MPa-tons in abaqus.

  1. Firstly we build 1/4 of the cochlear geometry (due to y and z symmetry). The shaded part is half of the basilar membrane. The geometrical parameters are the same with the published article as follows

    Ren L J, Hua C, Ding G H, et al. Three-dimensional finite element hydrodynamical modeling of straight and spiral cochlea[C]//AIP Conference Proceedings. AIP Publishing, 2018, 1965(1): 030003.

  2. Then we emerge the surfaces and mesh the solid fluid part, with the following domains
    • fluid of the vestibule
    • fluid of the scala vestibula (SV), above the BM
    • fluid of the scala vestibula (SV), others
    • fluid of the helicotrema note that using the reflect tool we can create the other half of the vestibula and SV and helicotrema.
  3. The area of the oval window is about 3.2mm$^2$, consider its shape as an ellipse of $a=1.25$ and $b=0.8$
  4. Use the reflect and scale tool to create the scala tympani (ST), note that the scala height ratio SV:ST=5:3; but currently we will set it as 1:1. If needed, we may adjust the nodes locations later.
  5. The edges/boundaries of the fluid domain are shell elements
  6. Define element and node sets
    • element sets are mostly included by the default component
      • for a better view, define fluid sets for the left / right side of the model (fluid_yplus and fluid_yminus)
      • elements of the ow (inner part) ow_inner, (outer part) ow_outer
    • node sets
      • to create the BM, define node sets as follows: line1,line2,line3,line4 
                  line 2
     ---------------------
  3  |                   | line 4   
  3  |                   | line 4
     ---------------------
          line 1
along the x positive and y positive directions
      • center of the oval window ow_center
    • surface data
      • the shell surface can be difined with shell elements
      • the fluid surface is defined fluid_surf_sv, fluid_surf_st
  7. Export to model.inp (save to MODEL01.inp)
In [37]:
%%html
<center><img width='400px' src='./MODEL01/figures/figure0001.png'> The box model geometry</center>
<center><img width='800px' src='./MODEL01/figures/figure0002.png'> The box model mesh</center>
The box model geometry
The box model mesh

Step 02: Process the data with fortran code

  1. use file __dic.log to define the work dictory, in this case, MODEL01
  2. use the fortran code readAbaqus3D.f90 to process the data
In [165]:
%%file ./readAbaqus3D.f90
! 本程序读取 abaqus 的文件
! 2015-07-27 by RenLiujie 
! 2018-08-22 需要读取logfile来进行读取文件夹
! 2018-08-23 添加了surface的读取来进行读取文件夹
! 2018-09-17 添加了读取 c3d6 单元的功能
! 2018-09-18 添加了读取 c3d4 单元的功能
! 2019-04-08 增加了补0的操作
! 2019-05-05 增加了ELSET的读取
program ReadAbaqus3D
implicit none

character(len=120)  :: logfile, folder, inputfile
character(len=180)  :: cmdline
character(len=4)    :: mode
integer             :: endoffile

integer             :: nod, ele, etype, bc, surf
integer             :: n1, n2, n3, n4, n5, n6, n7, n8, itmp, itmp2, i
integer             :: intarray(8)
real*8              :: x1, x2, x3
character(len=3)    :: charetype, charbc, charsurf

! 读取 log 文件, 获得当前工作目录
logfile = '__dic.log'
open(unit = 11, file = trim(logfile), status = 'old', action = 'read')
read(11,*) folder
close(unit=11)
n1 = system('cls')
write(*,*) '----------------- ABAQUS READER -----------------'
write(*,*) '| Working dictory >>> ', trim(folder)
write(*,*) '| Input file name >>> ', trim(folder),'.inp'
write(*,*) '-------------------------------------------------'
n1 = system('mkdir '//trim(folder));
write(*,*) '| Creating folder >>> ', trim(folder)
! 读取 abaqus3D inp 文件
open(unit = 11, file = './'//trim(folder)//'.inp', status = 'old', action = 'read')
endoffile = 0
etype = 0
nod = 0
ele = 0
bc = 0
surf = 0
do while(.true.)
    read(11,'(A180)',iostat = endoffile) cmdline
    if(endoffile.ne.0) exit
    if( cmdline(1:1).eq.'*' ) then
        mode = 'none'
        if(cmdline(2:5).eq.'NODE') then
            mode = 'NODE'
            write(*,*) cmdline
            close(unit = 12)
            write(inputfile,*) './',trim(folder),'/model_node.cnn'
            open(unit = 12, file = inputfile(2:), status = 'replace', action = 'write')
        endif
        if( (cmdline(2:5).eq.'ELEM') .and. (cmdline(15:16).eq.'S3') ) then ! 三角形单元
            mode = 'ELS3'
            etype = etype + 1
            write(*,*) trim(cmdline)
            close(unit = 12)
            write(charetype,'(i3.3)') etype
            write(inputfile,*) './',trim(folder),'/model_ele_',charetype,'.ele'
            write(*,*) trim(inputfile)
            open(unit = 12, file = inputfile(2:), status = 'replace', action = 'write')
        endif
        if( (cmdline(2:5).eq.'ELEM') .and. (cmdline(15:16).eq.'S4') ) then ! 四边形单元
            mode = 'ELS4'
            etype = etype + 1
            write(*,*) trim(cmdline)
            close(unit = 12)
            write(charetype,'(i3.3)') etype
            write(inputfile,*) './',trim(folder),'/model_ele_',charetype,'.ele'
            write(*,*) trim(inputfile)
            open(unit = 12, file = inputfile(2:), status = 'replace', action = 'write')
        endif
        if( (cmdline(2:5).eq.'ELEM') .and. (cmdline(15:18).eq.'C3D8') ) then ! 六面体单元
            mode = 'ELD8'
            etype = etype + 1
            write(*,*) trim(cmdline)
            close(unit = 12)
            write(charetype,'(i3.3)') etype
            write(inputfile,*) './',trim(folder),'/model_ele_',charetype,'.ele'
            write(*,*) trim(inputfile)
            open(unit = 12, file = inputfile(2:), status = 'replace', action = 'write')
        endif
        if( (cmdline(2:5).eq.'ELEM') .and. (cmdline(15:18).eq.'C3D6') ) then ! 六面体单元
            mode = 'ELD6'
            etype = etype + 1
            write(*,*) trim(cmdline)
            close(unit = 12)
            write(charetype,'(i3.3)') etype
            write(inputfile,*) './',trim(folder),'/model_ele_',charetype,'.ele'
            write(*,*) trim(inputfile)
            open(unit = 12, file = inputfile(2:), status = 'replace', action = 'write')
        endif
        if( (cmdline(2:5).eq.'ELEM') .and. (cmdline(15:18).eq.'C3D4') ) then ! 四面体单元
            mode = 'ELD4'
            etype = etype + 1
            write(*,*) trim(cmdline)
            close(unit = 12)
            write(charetype,'(i3.3)') etype
            write(inputfile,*) './',trim(folder),'/model_ele_',charetype,'.ele'
            write(*,*) trim(inputfile)
            open(unit = 12, file = inputfile(2:), status = 'replace', action = 'write')
        endif
        if( cmdline(2:5).eq.'NSET') then
            mode = 'NSET'
            bc = bc + 1
            write(*,*) trim(cmdline)
            close(unit = 12)
            write(charbc,'(i3.3)') bc
            write(inputfile,*) './',trim(folder),'/model_bc_',charbc,'.nset'
            write(*,*) trim(inputfile)
            open(unit = 12, file = inputfile(2:), status = 'replace', action = 'write')
        endif
        if( cmdline(2:5).eq.'ELSE') then
            mode = 'ELSE'
            bc = bc + 1
            write(*,*) trim(cmdline)
            close(unit = 12)
            write(charbc,'(i3.3)') bc
            write(inputfile,*) './',trim(folder),'/model_bc_',charbc,'.elset'
            write(*,*) trim(inputfile)
            open(unit = 12, file = inputfile(2:), status = 'replace', action = 'write')
        endif
        if( cmdline(2:5).eq.'SURF') then
            mode = 'SURF'
            surf = surf + 1
            write(*,*) trim(cmdline)
            close(unit = 12)
            write(charsurf,'(i3.3)') surf
            write(inputfile,*) './',trim(folder),'/model_surf_',charsurf,'.surf'
            write(*,*) inputfile
            open(unit = 12, file = inputfile(2:), status = 'replace', action = 'write')
        endif
    else
        if(mode.eq.'NODE') then
            nod = nod + 1
            read(cmdline,*) itmp, x1, x2, x3
            ! write(12,*) nod, x1, x2
            write(12,*) itmp,',',x1,',',x2,',',x3
        endif
        if(mode.eq.'ELS3') then
            ele = ele + 1
            read(cmdline,*) itmp, n1, n2, n3
            ! write(12,*) ele, n1, n2, n3
            write(12,*) itmp,',',n1,',',n2,',',n3
        endif
        if(mode.eq.'ELS4') then
            ele = ele + 1
            read(cmdline,*) itmp, n1, n2, n3, n4
            ! write(12,*) ele, n1, n2, n3
            write(12,*) itmp,',',n1,',',n2,',',n3,',',n4
        endif
        if(mode.eq.'ELD8') then
            ele = ele + 1
            read(cmdline,*) itmp, n1, n2, n3, n4, n5, n6, n7
            read(11,'(A120)',iostat = endoffile) cmdline
            read(cmdline,*) n8
            ! write(12,*) ele, n1, n2, n3
            write(12,*) itmp,',',n1,',',n2,',',n3,',',n4,',',n5,',',n6,',',n7,',',n8
        endif
        if(mode.eq.'ELD6') then
            ele = ele + 1
            read(cmdline,*) itmp, n1, n2, n3, n4, n5, n6
            ! write(12,*) ele, n1, n2, n3
            write(12,*) itmp,',',n1,',',n2,',',n3,',',n4,',',n5,',',n6
        endif
        if(mode.eq.'ELD4') then
            ele = ele + 1
            read(cmdline,*) itmp, n1, n2, n3, n4
            ! write(12,*) ele, n1, n2, n3
            write(12,*) itmp,',',n1,',',n2,',',n3,',',n4
        endif
        if(mode.eq.'NSET') then
            do itmp = 8,1,-1
                read(cmdline,*,iostat = itmp2) intarray(1:itmp)
                if(itmp2.eq.0) then 
                    do i=1,itmp; write(12,*) intarray(i); enddo
                    exit
                endif
            enddo
        endif
        if(mode.eq.'ELSE') then
            write(12,*) trim(cmdline)
        endif
        if(mode.eq.'SURF') then
            write(12,*) trim(cmdline)
        endif
    endif
enddo

close(unit = 11); close(unit = 12); close(unit = 13)
endprogram ReadAbaqus3D

Created file 'C:\Users\RenLi\__matlab\__model_cochlea\3d_model_abaqus\readAbaqus3D.f90'.
In [227]:
system('echo MODEL01>__dic.log');
system('gfortran readAbaqus3D.f90');
system('a');
system('del a.exe');

 ----------------- ABAQUS READER ----------------- 
 | Working dictory >>> MODEL01 
 | Input file name >>> MODEL01.inp 
 ------------------------------------------------- 
A subdirectory or file MODEL01 already exists. 
 | Creating folder >>> MODEL01 
 *NODE                                                                                                                                                                                
 *ELEMENT,TYPE=S4R,ELSET=shell-ow-inner 
  ./MODEL01/model_ele_001.ele 
 *ELEMENT,TYPE=S4R,ELSET=shell-ow-liga 
  ./MODEL01/model_ele_002.ele 
 *ELEMENT,TYPE=S4R,ELSET=shell-rw 
  ./MODEL01/model_ele_003.ele 
 *ELEMENT,TYPE=S4R,ELSET=shell-lam 
  ./MODEL01/model_ele_004.ele 
 *ELEMENT,TYPE=S4R,ELSET=shell-wall 
  ./MODEL01/model_ele_005.ele 
 *ELEMENT,TYPE=C3D8R,ELSET=fluid_sv_others 
  ./MODEL01/model_ele_006.ele 
 *ELEMENT,TYPE=C3D8R,ELSET=fluid_sv_bm 
  ./MODEL01/model_ele_007.ele 
 *ELEMENT,TYPE=C3D8R,ELSET=fluid_sv_heli 
  ./MODEL01/model_ele_008.ele 
 *ELEMENT,TYPE=C3D8R,ELSET=fluid_vest 
  ./MODEL01/model_ele_009.ele 
 *ELEMENT,TYPE=C3D8R,ELSET=fluid_st_others 
  ./MODEL01/model_ele_010.ele 
 *ELEMENT,TYPE=C3D8R,ELSET=fluid_st_bm 
  ./MODEL01/model_ele_011.ele 
 *ELEMENT,TYPE=C3D8R,ELSET=fluid_st_heli 
  ./MODEL01/model_ele_012.ele 
 *NSET, NSET=line1, UNSORTED 
  ./MODEL01/model_bc_001.nset 
 *NSET, NSET=line2, UNSORTED 
  ./MODEL01/model_bc_002.nset 
 *NSET, NSET=line3, UNSORTED 
  ./MODEL01/model_bc_003.nset 
 *NSET, NSET=line4, UNSORTED 
  ./MODEL01/model_bc_004.nset 
 *NSET, NSET=ow-center 
  ./MODEL01/model_bc_005.nset 
 *ELSET, ELSET=fluid-yp 
  ./MODEL01/model_bc_006.elset 
 *ELSET, ELSET=fluid-ym 
  ./MODEL01/model_bc_007.elset 
 *SURFACE, NAME = fluid-surf, TYPE = ELEMENT 
  ./MODEL01/model_surf_001.surf

Step03: The spiral shape

Use the function tospiral.m to generate the spiral shape of the model

In [228]:
%%file tospiral.m
function tospiral(spiral_or_not)
    folder = textread('__dic.log','%s'); folder = folder{1};
    nod = load([folder,'/model_node.cnn']);
    fid = fopen([folder,'/model_node_spiral.cnn'],'w');
    list = nod(:,1);
    nod = nod(:,2:4);
    locspr = nod;
    L = 35; r1 = L/3/pi;
    theta1 = -pi/2; dtheta=5*pi; theta2=theta1+dtheta;
    temp = [1 theta1; dtheta (theta2^2-theta1^2)/2]\[r1;L];
    a=temp(1); b=temp(2);
    r2=a+b*theta2;
    xc=r1*cos(theta1); yc=r1*sin(theta1);
    for ii=1:size(nod,1)
    x = nod(ii,1);
    if(x>=0)
        y=nod(ii,2);
        theta=(-a+sqrt(a^2+2*b*(a*theta1+b/2*theta1^2+x)))/b;
        r = a+b*theta;
        locspr(ii,1)=r*cos(theta)-y*cos(theta)-xc;
        locspr(ii,2)=r*sin(theta)-y*sin(theta)-yc;
        locspr(ii,3)=locspr(ii,3)+x/L*3.5;
    end
    if spiral_or_not==1
        fprintf(fid,'%d,%f,%f,%f\n',list(ii),locspr(ii,1),locspr(ii,2),locspr(ii,3));
    else
        fprintf(fid,'%d,%f,%f,%f\n',list(ii),nod(ii,1),nod(ii,2),nod(ii,3));
    end
    end
    fclose(fid);
    fprintf(1,'Mission Complete\n');
end

Created file 'C:\Users\RenLi\__matlab\__model_cochlea\3d_model_abaqus\tospiral.m'.
In [229]:
% tospiral(1) make a spiral shape, tospiral(0) keeps a strainght model
tospiral(1)

Mission Complete

Step 04: Create the Basilar Membrane

We use the MATLAB code generateBM.m to generate the BM nodes, mesh, thickness and material properties +

In [218]:
%%file generateBM.m
function generateBM(folder,line1,line2,line3,line4)

start = 1000000;
folder = textread('__dic.log','%s'); folder = folder{1};
nod = load([folder,'/model_node_spiral.cnn']); nod_list = nod(:,1);
line_inn = load([folder,'/',line1]);
line_out = load([folder,'/',line2]);
line_bas = load([folder,'/',line3]);
line_tip = load([folder,'/',line4]);

start = size(nod,1);

Ta = 0.0075; Tb = 0.0025;
Eya = 100; Eyb = 2.5; 
material = 'orth'; %iso
damping = 'stru';
structuraldamping = 0.15;
betadamping = 5e-6;
alphadamping = 200;

%               line-out  ---->
% ***********************************************************************           
% * ^                                                                   *           
% * | line-base                                                         *  line-tip 
% * |                                                                   *           
% *                                                                     *           
% ***********************************************************************           
%               line-inner  ---->
%

if isempty(line_tip)~=1; M = length(line_tip); end
if isempty(line_bas)~=1; M = length(line_bas); end
N = length(line_out); 
L = zeros(N,1);
%% calculate the length of the BM
for jj = 1:N
    n1 = line_inn(jj); n2 = line_out(jj);
    n1_index = find(nod_list==n1);
    n2_index = find(nod_list==n2);
    loc1 = nod(n1_index,2:end);
    loc2 = nod(n2_index,2:end);
    locc = (loc1+loc2)/2;
    if jj == 1;
        L(jj) = 0;
        locc0 = locc;
    else
        dL = norm(locc-locc0);
        L(jj)  = L(jj-1) + dL;
        locc0 = locc;
    end
end
fprintf('The length of the BM mid-line is about %f mm\n', L(end));

%% Generate the BM nodes
bmnodlist = zeros(N,M);
bmnod     = zeros(3,N,M);
fid_nod = fopen([folder,'/input-inst1-node-bm.inp'],'w'); 
fid_thk = fopen([folder,'/input-distribution-thick.inp'],'w');
fid_ele = fopen([folder,'/input-inst1-elem-bm.inp'],'w');
fid_set = fopen([folder,'/input-inst1-bm-set.inp'],'w');
fid_mat = fopen([folder,'/input-inst1-bm-mat.inp'],'w');
fid_sec = fopen([folder,'/input-inst1-bm-sec.inp'],'w');
fid_mid = fopen([folder,'/input-inst1-bm-mid.inp'],'w');

kk = start;
%% node & node thickness
if isempty(line_bas)~=1
    iistart = 2;
    for ii = 1:1
        for jj = 1:M
            n1 = line_bas(jj);
            n1_index = find(nod_list==n1);
            loc1 = nod(n1_index,2:end);
            bmnodlist(ii,jj) = n1;
            bmnod(:,ii,jj) = loc1;
        end
    end
else
    iistart = 1;
end
if isempty(line_tip)~=1
    iiend = N-1;
else
    iiend = N;
end

for ii = iistart:iiend
    n1 = line_inn(ii); n2 = line_out(ii);
    n1_index = find(nod_list==n1);
    n2_index = find(nod_list==n2);
    loc1 = nod(n1_index,2:end);
    loc2 = nod(n2_index,2:end); 
    bmnodlist(ii,1) = n1; bmnodlist(ii,M) = n2;
    bmnod(:,ii,1) = loc1; bmnod(:,ii,M) = loc2;
    for jj = 2:M-1
        locjj = loc1+(loc2-loc1)/(M-1)*(jj-1);
        bmnod(:,ii,jj) = locjj;
        kk = kk + 1;
        bmnodlist(ii,jj) = kk;
        fprintf(fid_nod,'%10d,%12.7e,%12.7e,%12.7e\n',kk,...
            locjj(1),locjj(2),locjj(3));
    end
end

if(iiend==N-1)
    for ii = N:N
        for jj = 1:M
            n1 = line_tip(jj);
            n1_index = find(nod_list==n1);
            loc1 = nod(n1_index,2:end);
            bmnodlist(ii,jj) = n1;
            bmnod(:,ii,jj) = loc1;
        end
    end
end
%% node thickness
for ii = 1:N
thick = Ta+(Tb-Ta)/L(end)*L(ii); 
for jj = 1:M
    n1 = bmnodlist(ii,jj);
    fprintf(fid_thk,'%10d,%12.7e\n',n1,thick);
end
fprintf(fid_mid,'%d\n',bmnodlist(ii,floor((1+M)/2)));
end

%% element
kk = start; Lend = L(end);
for ii = 1:N-1

    density = 1.0e-9;
    nuxy =0.03; nuxz = 0.3; nuyz = 0.3;
    locx = (L(ii)+L(ii+1))/2;
    Ey = Eya+(Eyb-Eya)*locx/Lend;
    Ex = Ey/10; Ez = Ey/10;
    Gxy = Ey/2; Gyz = Gxy; 
    Gxz = Ex/2/(1+nuxz);

    % *Orientation, name=ORI001, DEFINITION=NODES
    % node a, node b, node c
    % normaldirection, angle
    % *shell section, ... orientation=ORI001
    % *Elastic, type=ENGINEERING CONSTANTS
    % 1.,  2.,  3., 0.1, 0.2, 0.3, 10., 20.
    % 30.,   
    nc = bmnodlist(ii,(M+1)/2); 
    nb = bmnodlist(ii,M); 
    na = bmnodlist(ii+1,(M+1)/2);
    fprintf(fid_set,'*elset, elset=bm%010d\n',kk);
    fprintf(fid_sec,'*orientation, name=ori-bm%010d, definition=nodes\n',kk);
    fprintf(fid_sec,'%d,%d,%d\n', na,nb,nc);
    fprintf(fid_sec,'3, 1.0\n');
    fprintf(fid_sec,'*shell section, elset=bm%010d, material=mat-bm%010d, orientation=ori-bm%010d, nodal thickness\n',...
        kk,kk,kk);
    fprintf(fid_mat,'*material, name=mat-bm%010d\n',kk);
    if damping(1)=='s'
        fprintf(fid_mat,'*damping, structural=%10.3e\n',structuraldamping);
    else
        fprintf(fid_mat,'*damping, beta=%10.3e, alpha=%10.3e\n',betadamping, alphadamping);
    end
    fprintf(fid_mat,'*density\n');
    fprintf(fid_mat,'%12.7e\n',density);
    if material(1) == 'o'
        fprintf(fid_mat,'*elastic, type=engineering constants\n');
        fprintf(fid_mat,'%12.7e,%12.7e,%12.7e,%12.7e,%12.7e,%12.7e,%12.7e,%12.7e\n',Ex,Ey,Ez,...
            nuxy,nuxz,nuyz,Gxy,Gxz);
        fprintf(fid_mat,'%12.7e,\n',Gyz);
    else
        fprintf(fid_mat,'*elastic \n');
        fprintf(fid_mat,'%12.7e,%12.7e\n',Ey,nuxz);
    end
    
    for jj=1:M-1
        kk = kk + 1;
        n1 = bmnodlist(ii,jj); n2 = bmnodlist(ii+1,jj);
        n3 = bmnodlist(ii+1,jj+1); n4 = bmnodlist(ii,jj+1);
        fprintf(fid_ele,'%8d,%8d,%8d,%8d,%8d\n',kk,n1,n2,n3,n4);
        fprintf(fid_set,'%10d\n',kk);
    end
end
fclose(fid_nod);
fclose(fid_thk);
fclose(fid_ele);
fclose(fid_set);
fclose(fid_mat);
fclose(fid_sec);
fprintf(1,'Mission Complete !\n');
end

Created file 'C:\Users\RenLi\__matlab\__model_cochlea\3d_model_abaqus\generateBM.m'.
In [230]:
folder = textread('__dic.log','%s'); folder = folder{1};
line1 = 'model_bc_001.nset';
line2 = 'model_bc_002.nset';
line3 = 'model_bc_003.nset';
line4 = 'model_bc_004.nset';
generateBM(folder,line1,line2,line3,line4);

The length of the BM mid-line is about 35.220971 mm
Mission Complete !

Step 5: AC simulation code

generate an abaqus .inp code for calculation!

In [234]:
%%file MODEL01/model01-ac-spiral.inp
*Heading
*preprint, echo=yes, mode=yes, history=no, contact=no
** --------------------------------------------------------
**          M O D E L - D E F I N I T I O N                
** --------------------------------------------------------
*part, name=part1
*end part
*assembly, name=asse1
*instance, name=inst1, part=part1
** node and elements
*node
*include, input=model_node_spiral.cnn
*include, input=input-inst1-node-bm.inp

*ELEMENT, TYPE=AC3D8, ELSET=FLUID
*include, input=model_ele_006.ele
*include, input=model_ele_007.ele
*include, input=model_ele_008.ele
*include, input=model_ele_009.ele
*include, input=model_ele_010.ele
*include, input=model_ele_011.ele
*include, input=model_ele_012.ele
*ELEMENT, TYPE=S4R, ELSET=OW-INNER
*include, input=model_ele_001.ele
*ELEMENT, TYPE=S4R, ELSET=OW-LIGA
*include, input=model_ele_002.ele
*ELEMENT, TYPE=S4R, ELSET=RW
*include, input=model_ele_003.ele
*ELEMENT, TYPE=S4R, ELSET=WALL02
*include, input=model_ele_004.ele
*ELEMENT, TYPE=S4R, ELSET=WALL01
*include, input=model_ele_005.ele

*ELEMENT, TYPE=S4R, ELSET=BM
*include, input=input-inst1-elem-bm.inp
** nodal thickness
*nodal thickness
*include, input=input-distribution-thick.inp
** bm-element-mats
** *shell section, elset=bm, material=mat-bm, nodal thickness
*include, input=input-inst1-bm-set.inp
*include, input=input-inst1-bm-sec.inp
** sections
*SHELL SECTION, ELSET=WALL01, MATERIAL=MAT-BONE
1,
*SHELL SECTION, ELSET=WALL02, MATERIAL=MAT-BONE
1,
*SHELL SECTION, ELSET=OW-INNER, MATERIAL=MAT-OW
1,
*SHELL SECTION, ELSET=OW-LIGA, MATERIAL=MAT-OW
1,
*SHELL SECTION, ELSET=RW, MATERIAL=MAT-RW
1,
*SOLID SECTION, ELSET=FLUID, MATERIAL=MAT-FLUID
,
** sets and elements and surfaces
*SURFACE, NAME=SURF-FLUID, TYPE=ELEMENT
*include, input=model_surf_001.surf
*SURFACE, NAME=SURF-SOLID, TYPE=ELEMENT
OW-INNER, SNEG
OW-LIGA,  SNEG
RW,       SNEG
BM,       SPOS
WALL01,   SNEG
WALL02,   SPOS
*NSET, NSET=OW-CENTER
*include, input=model_bc_005.nset
*ELSET, ELSET=WALL
WALL01
WALL02
*ELSET, ELSET=fluid_yplus
*include, input=model_bc_006.elset
*ELSET, ELSET=fluid_yminus
*include, input=model_bc_007.elset
*NSET, NSET=WALLNODE, ELSET=WALL
*NSET, NSET=BM-MIDLINE
*include, input=input-inst1-bm-mid.inp
*RIGID BODY, REF NODE=OW-CENTER, ELSET=OW-INNER
*end instance
*tie, name=tie1, adjust=yes, position tolerance=0.01, no thickness
inst1.surf-fluid, inst1.surf-solid
*end assembly
** --------------------------------------------------------
**          M O D E L - S E T T I N G S                    
** --------------------------------------------------------
** MATERIALS
*include, input=input-inst1-bm-mat.inp
*material, name=mat-bone
*Damping, structural=0.05
*density
2.200e-9
*elastic
14100
*material, name=mat-ow
*Damping, structural=0.05
*density
1.200e-9
*elastic
1.0, 0.3
*material, name=mat-rw
*Damping, structural=0.05
*density
1.200e-9
*elastic
1.0, 0.3
*material, name=mat-fluid
*acoustic medium
1E6,
*acoustic medium, volumetric drag
0.0, 500.
*density
1.000e-9
** BOUNDARIES
*boundary
INST1.OW-CENTER, 1 ,1
INST1.OW-CENTER, 2 ,2
INST1.OW-CENTER, 4 ,4
INST1.OW-CENTER, 5 ,5
INST1.OW-CENTER, 6 ,6
INST1.WALLNODE, 1, 1
INST1.WALLNODE, 2, 2
INST1.WALLNODE, 3, 3
INST1.WALLNODE, 4, 4
INST1.WALLNODE, 5, 5
INST1.WALLNODE, 6, 6
** STEP CALC
*step, name=step001, nlgeom=no, perturbation
*Steady State Dynamics, direct, friction damping=NO
2000, 2000, 1, 1.
*BOUNDARY, REAL
INST1.OW-CENTER,3,3,-1.0
*BOUNDARY, imaginary
INST1.OW-CENTER,3,3, 0.0
** LOAD
**Dsload, real
**inst1.press, P, 1.0
*output, field, variable=preselect
*output, history, variable=preselect
*end step

Created file 'C:\Users\RenLi\__matlab\__model_cochlea\3d_model_abaqus\MODEL01\model01-ac-spiral.inp'.
In [235]:
%%html
<center><img width='1200px' src='./MODEL01/spiral-ac-2000Hz.gif'> spiral model</center>
spiral model